[R-sig-Geo] reading .grd file in R?

[Barry Rowlingson]

On Tue, Oct 12, 2010 at 9:04 PM,  <govindas at 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