[R-sig-Geo] inconsistencies between readGDAL and other software (GrADS) when reading GRIB1 files

[Ariel Ortiz-Bobea]

Hi,

I've been importing GRIB1 files into R through the readGDAL command in the
GDAL package.

However, the GRIB1 files are multi-band and during the import process I
"lose" the variable names and get generic names: band1, band2, etc. 

I have tried to circumvent this by using other software (GrADS) to obtain
summary statistics for each variable (of which I know the names) and then
compare these to summary statistics of each variable/band in R. Matching
summary statistics would allow me to get back the real names. 

However, there are discrepancies between summary statistics. More precisely,
there are differences between minimums, means and maximums between both
approaches AND these differences are not of the same magnitude (ruling out a
homogenous shift in values).

Moreover, GrADS gives temperatures in Kelvin and readGDAL seems to be
transforming them to Celsius BUT the difference between minimum, maximum and
means is not consistently 273.15 as one would expect. For instance, for soil
temperature the difference is 271.441 (for the mean), 270.514 (for the max)
and 273.149 (for the min), but for Dew point it is 269.436 (for the mean),
272.58 (for the max) and 273.349 (for the min). So it is not even consistent
across different temperature variables.

I wonder if this has to do with how GrADS or readGDAL read these GRIB1 files
(or the way I'm importing this)... would there be any hidden processing
going on that would explain this? Would using netCDF format files rather
than GRIB1 help in any way? Any suggestions or comments on how to import
this data in a manner in a proper way?

Thanks a lot in advance for any guidance on this,

Ariel

Here is my example file:
http://dl.dropbox.com/u/45311184/US48_merged_AWIP32.1985120100.RS.flx (in
GRIB1 format)

R code:
-------
# import
   flx<-  readGDAL("1985/US48_merged_AWIP32.1985120100.RS.flx")
# fix missing values (imported as 9999)
   for (n in 1:10 ){
      is.na(flx[[n]]) <- flx[[n]] > 9000
   }
# print summary statistics for each variable/band
   for (n in 1:10 ){
      print(summary(flx[[n]])[c(4,6,1)])
   }

GrADS:
--------
code: http://dl.dropbox.com/u/45311184/gradscode.txt
files to be included in the same folder with the .flx file 
http://dl.dropbox.com/u/45311184/US48_merged_AWIP32.1985120100.RS.flx.ctl
http://dl.dropbox.com/u/45311184/US48_merged_AWIP32.1985120100.RS.flx.idx


-----
Ariel Ortiz-Bobea
PhD Candidate in Agricultural & Resource Economics
University of Maryland - College Park
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/inconsistencies-between-readGDAL-and-other-software-GrADS-when-reading-GRIB1-files-tp7006160p7006160.html
Sent from the R-sig-geo mailing list archive at Nabble.com.