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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
