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

Fri Nov 18 21:29:05 CET 2011

Thanks for the references Roger,

The "shift bug" had escaped me! Now the overlay between polygons and grids fits better visually.

On the other hand, GDALinfo("mygrib.grb", returnRAT=TRUE) returns empty names for all bands and a warning for each of them: "statistics not supported by this driver".

I found a workaround strategy converting a reference GRIB file to netCDF format (using the ncl_convert2nc software) and then reading it with the "ncdf" package. This way I could obtain the metada (abbreviations, long names, units) and compare the summary statistics with those of the unnamed bands from RGDAL. The summary statistics match, so I was able to obtain the names this way....So I'm guessing that it is GrADS who's doing some extrapolations/processing, at least in my case, explaining the discrepancies. 

It is worth noting that readGDAL transforms temperatures in Kelvin to Celsius when reading GRIB1 files...all other units in my case remain untouched.

One other thing is that the "rgdal" and "ncdf" packages in R, as well as GrADS for that matter, return bands in different orders. So you really have to go through this matching process to guess variable names/identities when reading GRIB1 files with many bands using RGDAL (version 1.8.0, released 2011/01/12)

Thanks again,

Ariel OB

________________________________________
From: Roger Bivand [via R-sig-geo] [ml-node+s2731867n7007977h13 at n2.nabble.com]
Sent: Friday, November 18, 2011 6:56 AM
To: Ariel Ortiz Bobea
Subject: Re: inconsistencies between readGDAL and other software (GrADS) when reading GRIB1 files

On Thu, 17 Nov 2011, Ariel Ortiz-Bobea wrote:

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

The readGDAL() function in the rgdal package uses GDAL and its drivers, so
see:

http://www.gdal.org/frmt_grib.html

for details. Note also:

https://stat.ethz.ch/pipermail/r-sig-geo/2011-January/010760.html
http://grass.osgeo.org/wiki/GRIB

for a bug in the GDAL driver and a workaround, and other GRIB threads on
this list.

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

As far as I can see, band names are not supported in GDAL in:

http://www.gdal.org/classGDALDataset.html

or

http://www.gdal.org/classGDALRasterBand.html

but may be in http://www.gdal.org/classGDALRasterAttributeTable.html.

Could you please try running:

GDALinfo("mygrib.grb", returnRAT=TRUE)

to see if you get the band names?

Roger

> 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.
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]</user/SendEmail.jtp?type=node&node=7007977&i=0>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: [hidden email]</user/SendEmail.jtp?type=node&node=7007977&i=1>

_______________________________________________
R-sig-Geo mailing list
[hidden email]</user/SendEmail.jtp?type=node&node=7007977&i=2>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Economic Geography Section
Department of Economics
Norwegian School of Economics and Business Administration
Helleveien 30
N-5045 Bergen, Norway

________________________________
If you reply to this email, your message will be added to the discussion below:
http://r-sig-geo.2731867.n2.nabble.com/inconsistencies-between-readGDAL-and-other-software-GrADS-when-reading-GRIB1-files-tp7006160p7007977.html
To unsubscribe from inconsistencies between readGDAL and other software (GrADS) when reading GRIB1 files, click here<http://r-sig-geo.2731867.n2.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=7006160&code=YW9ydGl6Ym9iZWFAYXJlYy51bWQuZWR1fDcwMDYxNjB8LTYxMzExNzA2NQ==>.
NAML<http://r-sig-geo.2731867.n2.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.InstantMailNamespace&breadcrumbs=instant+emails%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>

-----
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-tp7006160p7009592.html
Sent from the R-sig-geo mailing list archive at Nabble.com.