readGDAL() and HDF5 files
On Sat, 24 Jan 2009 12:56:13 +0100 (CET),
Roger Bivand <Roger.Bivand at nhh.no> wrote:
On Fri, 23 Jan 2009, Sebastian P. Luque wrote:
On Sat, 24 Jan 2009 08:11:15 +1100, Michael Sumner <mdsumner at utas.edu.au> wrote:
Hi Sebastian, When you say the geographic coordinates are not a
regular grid - is it that the actual grid is Mercator but the NetCDF file stores an X and Y vector separately for each unique longitude and latitude? I've seen this many times, but never with enough metadata to determine that without guessing. I've seen some documents that refer to the source grid as being in Mercator, but never any that mention explicitly the method used to generate the NetCDF file from those.
If this is the case for your data, I've had success by figuring out an offset/scale value that work sufficiently by assuming a Mercator grid
Specifically this one but sometimes with an extra X offset to overcome hemisphere shift.
CRSargs(CRS("+proj=merc"))
I don't have the details available today, but I can dig them up if that sounds promising. Also, I'd be interested to hear any details you have about the grids you are using, whether they use this Mercator-kludge or not.
Thanks for the feedback.
You're right in that the geographic coordinates are not a regular grid, but separate data subsets in the NetCDF file, with the actual grid being in a polar stereographic projection (the one in my first post). The actual coordinates of the grid are not available in any file (NetCDF nor HDF5), and neither these nor the geographical coordinates in the HDF5 files, but that doesn't matter in my case because I need to subset a smaller area, grid at an appropriate resolution using the geographical coordinates, and reproject to a different projection for my area. After all, it's simple enough to access the geographical coordinates from any of the NetCDF files to map the grid onto such a coordinate system.
Therefore, for working with these grids in R, I'd be happy if I could just load the grid correctly using the arbitrary grid coordinates that readGDAL() builds from the HDF5 file. I haven't managed to read the NetCDF files into R (R segfaults), so things look better with HDF5 for working in R with these data.
Please do provide full version data for R, rgdal, GDAL, and the GDAL drivers, including the HDF5 libraries. I did - please do the same: all tests on Windows, OSGeo4W GDAL 1.5.4 with released drivers:
I did provide most of this information in my first post, although I
forgot to include info on non-R software, so here it is again:
R> sessionInfo()
R version 2.8.1 (2008-12-22)
x86_64-pc-linux-gnu
locale:
LC_CTYPE=en_CA.UTF-8;LC_NUMERIC=C;LC_TIME=en_CA.UTF-8;LC_COLLATE=en_CA.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_CA.UTF-8;LC_PAPER=en_CA.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_CA.UTF-8;LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rgdal_0.6-5 sp_0.9-29 slmisc_0.7.0 lattice_0.17-20
loaded via a namespace (and not attached):
[1] grid_2.8.1
and Debian sid libgdal1-1.5.2-3, which lags behind slightly the upstream
version. I'm also using Debian sid's libhdf5-serial-1.6.6-0 for the
HDF5 libraries. OSGeo4W is only for MS Windows IIUC. So we may not be
able to properly compare our results.
For what it's worth, this is what I tried with one of the NetCDF files
(using Debian sid's libnetcdf4 1:3.6.2-3.1, where ncdump-hdf is
equivalent to the upstream ncdump facility):
---<--------------------cut here---------------start------------------->---
$ ncdump-hdf -h ice_conc_nh_200901151200.nc
netcdf ice_conc_nh_200901151200 {
dimensions:
time = 1 ;
ni = 760 ;
nj = 1120 ;
variables:
long time(time) ;
time:long_name = "reference time of sea ice file" ;
time:units = "seconds since 1981-01-01 00:00:00" ;
float lat(nj, ni) ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
float lon(nj, ni) ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
char polar_stereographic ;
polar_stereographic:grid_mapping_name = "polar_stereographic" ;
polar_stereographic:straight_vertical_longitude_from_pole = -45.f ;
polar_stereographic:latitude_of_projection_origin = 90.f ;
polar_stereographic:standard_parallel = 70.f ;
polar_stereographic:false_easting = -3850.f ;
polar_stereographic:false_northing = 5850.f ;
polar_stereographic:proj4_string = "+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45" ;
short ice_concentration(time, nj, ni) ;
ice_concentration:long_name = "sea ice concentration" ;
ice_concentration:standard_name = "ice_concentration" ;
ice_concentration:units = "%" ;
ice_concentration:coordinates = "lon lat" ;
ice_concentration:grid_mapping = "polar_stereographic" ;
ice_concentration:source = "EUMETSAT OSI SAF" ;
ice_concentration:missing_value = -32767s ;
ice_concentration:_FillValue = -32767s ;
ice_concentration:valid_min = 0.f ;
ice_concentration:valid_max = 100.f ;
ice_concentration:scale_factor = 0.0099999998f ;
ice_concentration:add_offset = 0.f ;
byte confidence_level(time, nj, ni) ;
confidence_level:long_name = "confidence level" ;
confidence_level:units = "1" ;
confidence_level:coordinates = "lon lat" ;
confidence_level:grid_mapping = "polar_stereographic" ;
confidence_level:missing_value = '\0' ;
confidence_level:_FillValue = '\0' ;
confidence_level:valid_min = '\1' ;
confidence_level:valid_max = '\5' ;
confidence_level:comment = "0 Unprocessed; 1 Erroneous; 2 Unreliable; 3 Acceptable; 4 Good; 5 Excellent" ;
short quality_index(time, nj, ni) ;
quality_index:long_name = "quality index" ;
quality_index:units = "1" ;
quality_index:coordinates = "lon lat" ;
quality_index:grid_mapping = "polar_stereographic" ;
quality_index:missing_value = 0s ;
quality_index:_FillValue = 0s ;
quality_index:comment = "Contains bitmap with quality flags, see Users Manual for details" ;
// global attributes:
:title = "Total Sea Ice Concentration from EUMETSAT OSI SAF" ;
:conventions = "CF-1.0" ;
:netcdf_version_id = "3.5.1 of Feb 7 2008 11:53:33 $" ;
:creation_date = "2009-01-16" ;
:product_version = "2.1" ;
:software_version = "3.2" ;
:reference = "OSI SAF Sea Ice Product Manual, Andersen S., Breivik L.A., Eastwood S., God?y ?., Lind M., Porcires M., Schyberg H., v3.5, January 2007" ;
:comment = "For more information about the OSI SAF Sea Ice products, see http://saf.met.no" ;
:sensor = "Multi sensor" ;
:spatial_resolution = "10.0km" ;
:area = "Northern Hemipshere" ;
:southernmost_latitude = 31.08831f ;
:northernmost_latitude = 89.934723f ;
:westernmost_longitude = -180.f ;
:easternmost_longitude = 179.92558f ;
:start_date = "2009-01-15 UTC" ;
:start_time = "00:00:00 UTC" ;
:stop_date = "2009-01-15 UTC" ;
:stop_time = "23:59:59 UTC" ;
:field_type = "daily" ;
:institution = "EUMETSAT OSI SAF" ;
:institution_references = "http://www.osi-saf.org" ;
:contact = "osisaf-manager at met.no" ;
:operational_status = "operational" ;
}
$ gdalinfo ice_conc_nh_200901151200.nc
Segmentation fault
---<--------------------cut here---------------end--------------------->---
As can be seen, these files have more datasets. Presumably the segfault
I'm getting in R stems from the same problem with gdal above:
---<--------------------cut here---------------start------------------->---
R> ncdata <- "ice_conc_nh_200901151200.nc"
R> p4s <- NA
R> ice <- readGDAL(ncdata, p4s=p4s)
*** caught segfault ***
address 0x1050, cause 'memory not mapped'
Traceback:
1: .Call("RGDAL_OpenDataset", as.character(filename), TRUE, PACKAGE = "rgdal")
2: .local(.Object, ...)
3: initialize(value, ...)
4: initialize(value, ...)
5: new("GDALReadOnlyDataset", filename)
6: GDAL.open(fname)
7: readGDAL(ncdata, p4s = p4s)
Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
---<--------------------cut here---------------end--------------------->---
although note that I don't know how to tell readGDAL to access the
dataset that contains the actual grid I want (ice_concentration) for
this file format.
Anyway, based on your analysis below, it does seem as if the file header
is not providing metadata that GDAL can actually understand. Otherwise
gdalinfo should not segfault like that. I also get the same messed up
image I get in R if I use gdaltranslate to see a TIFF image of the file.
Thanks for looking into this. I can live using your data "molesting"
approach for now! :-)
Cheers,
Seb
and unreleased Windows binary rgdal_0.6-6 from:
Then at least we are comparing like with like. Nothing should segfault under any circumstances (though the GDAL HDF5 driver segfaults for me when asking for non-existent data[01]). A fuller report on the NetCDF case would be helpful.
My feeling at the moment is that the file is not self-documenting itself in a portable way, because GDAL does read the data in an order that isn't what the file header claims.
The example I showed (and that Roger reproduced) builds a grid with the wrong dimensions.
The dimensions are what the file declares - I suspect that software within the originator institution knows that they are reversed, and does the same, but this isn't portable. I further suspect that the 1120x1 block size is not helpful, and probably should be 760x1. See below for analysis.
This is what the header of the whole file (this is a new file, with the same structure as the one I showed) says (using HDF5 tools, http://hdf.ncsa.uiuc.edu/HDF5):
---<--------------------cut
here---------------start------------------->--- $ h5dump -H
ice_conc_nh_200901011200.hdf HDF5 "ice_conc_nh_200901011200.hdf" {
GROUP "/" { GROUP "Data" { DATASET "data[00]" { DATATYPE
H5T_IEEE_F32LE DATASPACE SIMPLE { ( 760, 1120 ) / ( 760, 1120 ) }
ATTRIBUTE "Description" { DATATYPE H5T_STRING { STRSIZE 16; STRPAD
H5T_STR_NULLTERM; CSET H5T_CSET_ASCII; CTYPE H5T_C_S1;
}
DATASPACE SIMPLE { ( 1 ) / ( 1 ) }
}
}
}
DATASET "Header" { DATATYPE H5T_COMPOUND { H5T_STRING { STRSIZE 32;
STRPAD H5T_STR_NULLTERM; CSET H5T_CSET_ASCII; CTYPE H5T_C_S1; }
"source"; H5T_STRING { STRSIZE 16; STRPAD H5T_STR_NULLTERM; CSET
H5T_CSET_ASCII; CTYPE H5T_C_S1; } "product"; H5T_STRING { STRSIZE 16;
STRPAD H5T_STR_NULLTERM; CSET H5T_CSET_ASCII; CTYPE H5T_C_S1; }
"area"; H5T_STRING { STRSIZE 128; STRPAD H5T_STR_NULLTERM; CSET
H5T_CSET_ASCII; CTYPE H5T_C_S1; } "projstr"; H5T_STD_U32LE "iw";
H5T_STD_U32LE "ih"; H5T_STD_U16LE "z"; H5T_IEEE_F32LE "Ax";
H5T_IEEE_F32LE "Ay"; H5T_IEEE_F32LE "Bx"; H5T_IEEE_F32LE "By";
H5T_STD_U32LE "year"; H5T_STD_U16LE "month"; H5T_STD_U16LE "day";
H5T_STD_U16LE "hour"; H5T_STD_U16LE "minute";
}
DATASPACE SIMPLE { ( 1 ) / ( 1 ) }
}
}
}
---<--------------------cut
here---------------end--------------------->---
and the grid should have 1120 rows and 760 columns, as displayed by gdalinfo on the actual data:
No, see below, you are reversing axes in gdalinfo:
---<--------------------cut here---------------start------------------->--- $ gdalinfo HDF5:"ice_conc_nh_200901011200.hdf"://Data/data[00] Driver: HDF5Image/HDF5 Dataset Files: none associated Size is 1120, 760 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], TOWGS84[0,0,0,0,0,0,0], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9108"]], AXIS["Lat",NORTH], AXIS["Long",EAST], AUTHORITY["EPSG","4326"]] Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 760.0) Upper Right ( 1120.0, 0.0) Lower Right ( 1120.0, 760.0) Center ( 560.0, 380.0) Band 1 Block=1120x1 Type=Float32, ColorInterp=Undefined Metadata: data[00]:Description=Ice Conc ---<--------------------cut here---------------end--------------------->---
but readGDAL() apparently read the dimensions in the opposite order:
No, exactly as the GDAL driver does, 760 rows and 1120 columns, left is 0, right is 1120, upper is 0, lower is 760.
upper left (0,0) upper right (1120, 0)
lower left (0,760) lower right (1120, 760)
So either the GDAL driver is broken, or the data in the file is not in sync with its metadata.
So far I get your figure by molesting the input metadata:
ice <- readGDAL("HDF5:\"conc_200901011200.hdf\"://Data/data[00]",
p4s=as.character(NA)) is.na(ice$band1) <- ice$band1 < 0 vice <-
ice$band1 mice <- matrix(vice, ncol=760, byrow=TRUE)
image(t(mice[1120:1,]), asp=1)
As I said before, the input HDF5 file has a completely wrong description of its own projection, so my conclusion would be that it is so badly configured as not to be usable in a portable way, since both its geometry and projection are declared in error.
Hope this helps,
Roger
---<--------------------cut here---------------start------------------->---
R> summary(ice)
Object of class SpatialGridDataFrame Coordinates: min max x 0 1120 y 0 760 Is projected: NA proj4string : [NA] Number of points: 2 Grid attributes: cellcentre.offset cellsize cells.dim x 0.5 1 1120 y 0.5 1 760 Data attributes: Min. 1st Qu. Median Mean 3rd Qu. Max. -32800 -99 -99 -117 0 100 ---<--------------------cut here---------------end--------------------->---
Does this make sense? and could this be the problem? The data that are read into the `ice' object do look ok via summary(ice), but the way they're mapped onto the grid does not.
Cheers,
-- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no
Cheers,
Seb