My suggestion is to convert your .grd files to self-descriptive ones, for
example netcdf. Considering that you have some experience with GrADS, a helpful
tutorial is available
on http://cookbooks.opengrads.org/index.php?title=Recipe-017:_Subsetting_and_reformating_MERRA_data_with_Lats4d#Examples,
but the same results could be reached by using NCO or CDO.
This procedure would make easier the process of importing the data into R,
using for example the packages "ncdf" or "RNetCDF".
Hope this helps,
Thiago.
--- On Tue, 12/10/10, govin...@msu.edu <govin...@msu.edu> wrote:
From: govin...@msu.edu <govin...@msu.edu>
Subject: Re: [R-sig-Geo] reading .grd file in R?
To: "Barry Rowlingson" <b.rowling...@lancaster.ac.uk>
Cc: r-sig-geo@stat.math.ethz.ch
Date: Tuesday, 12 October, 2010, 20:35
thanks a lot! this indeed has helped me a lot! thanks again!
Thanks,
Mahalakshmi Quoting Barry Rowlingson <b.rowling...@lancaster.ac.uk>:
> On Tue, Oct 12, 2010 at 9:04 PM, <govin...@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]]
-----Inline Attachment Follows-----
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
