

thanks a lot! this indeed has helped me a lot! thanks again!  

Thanks,
Mahalakshmi Quoting Barry Rowlingson <b.rowlingson@lancaster.ac.uk>:

> On Tue, Oct 12, 2010 at 9:04 PM,  <govindas@msu.edu> wrote:
>>
>>
>> Thanks a lot for cracking this Barry.... But, I am still not able to
>> completely understand this as I am not so experienced in R! I got an image
>> display with 4 plots in it. If I am not wrong, are these the plots of
>> rainfall Values (or the gridded dataset) for the 1st four days? can you
>> explain this function's overall purpose, so that I could know if I got it
>> right!
>
>  All the metadata comes from the .ctl file, which looks like this:
>
> DSET ^rf0.5_1975.grd
> options template
> TITLE 0.5 degranalyzed normal grids
> UNDEF -999.0
> XDEF  69  LINEAR  66.5 0.5
> YDEF  65  LINEAR  6.5 0.5
> ZDEF   1 linear 1 1
> TDEF 365 LINEAR 1jan1975 1DY
> * FOR LEAP YEARS CHANGE NO. OF RECORDS TO 366 *
> VARS  1
> rf 0 99 GRIDDED RAINFALL
> ENDVARS
>
>  I don't understand all of it, but that's where I get the X and Y size
> of the grid, and also how I know its 365 time grids. The UNDEF line
> tells me that -999.0 is missing data or no data in some sense.
>
>> I have some other doubts as well which i have mentioned below!
>
>  Yes, it wasn't exactly well commented code! Here goes:
>
>  > c = file(f,open="rb")
>  > seek(c,x*y*(tm-1)*size)                        #what is its purpose?
>
>  The file() function opens a 'connection' to the file, and then 'seek'
> skips a number of bytes into the connection so that the next readBin
> starts there. In this case we skip 69*65*(tm-1)*4 bytes, which gets us
> to the tm'th time grid (there are 4 (the size of each number) times 69
> times 65 bytes per time grid).
>
>  > m=matrix(readBin(c,type,x*y,size),x,y)
>
>  This then reads the binary data as 'double' values, 69*65 of them,
> and puts them in a 69x65 matrix.
>
>  > close(c)
>  > m[abs(m-NAvalue)<.Machine$double.eps] = NA             # a datavalue in
>  > matrix is checked for its lowest possible value and assigned "NA" 
> if its so?
>
> Its fairly common to use an extreme low or high value for missing
> data. In this case the value is specified in the .ctl file (see
> above). I use a test against .Machine$double.eps because testing for
> exact equality with floating point numbers doesn't always work as
> expected.
>
>  > par(mfrow=c(2,2))
>  > for(p in 1:4){
>  > image(readTime("rf0.5_1975.grd",p))
>  > }
>
>  The 'par' splits the window into 4 and then the loop does indeed plot
> the data for time=1 to time=4.
>
>  The next thing to do (I really should go to bed now) is assign X and
> Y coordinates. In the .ctl they are:
>
> XDEF  69  LINEAR  66.5 0.5
> YDEF  65  LINEAR  6.5 0.5
>
>  so I think they are:
>
>  x = seq(66.5,len=69, by=0.5)
>  y = seq(6.5,len=65, by=0.5)
>
> You can then create a list of x, y and z values this way:
>
>  r3=list(x=seq(66.5,len=69, by=0.5),y=seq(6.5,len=65,
> by=0.5),z=readTime("rf0.5_1975.grd",3))
>
> and get a properly georeferenced image matrix, which you could overlay
> on a map of india from the 'maps' package:
>
>  library(maps)
>  map("world","india")
>  image(r3,add=TRUE)
>
>  which is nice.
>
>  Info about the grads metadata file is here:
>
> http://www.iges.org/grads/gadoc/descriptorfile.html
>
> Barry
>
	[[alternative HTML version deleted]]

