reading HDF5 files
On Tue, 14 Aug 2007, Paul Hiemstra wrote:
Mike,
I tried, and failed, to directly load the HDF5 files into R using rgdal.
The next thing I tried was to convert the .h5 file to a GTiff using
FWTools. This seemed to work, but when I tried to load the file into R
through rgdal I got the following error message:
test = readGDAL("regenradar_hdf5/RAD_NL21_PCP_NA_200707200000.tif")
test = readGDAL("RAD_NL21_PCP_NA_200707200000.tif")
RAD_NL21_PCP_NA_200707200000.tif has GDAL driver GTiff
and has 256 rows and 256 columns
Closing GDAL dataset handle 0x0226a220... destroyed ... done.
Error in `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
Geographical CRS given to non-conformant data
Please state the gdalinfo of the input HDF5, and the output gTiff. Your use of gdal_translate seems to be mangling the (IIRC) projected coordinate reference system of the HDF5 to a geographical coordinate reference system, with the result that the coordinate values are "out of bounds". I think you need to show the full gdal_translate command, and to look at the help page: http://www.gdal.org/gdal_translate.html and try out the -a_srs option. Maybe also look at the -sds option if you are not sure that the subset is being chosen correctly. Once gdalinfo on the output gTiff shows a geographical coordinate reference system, rgdal will work (unless this is near-polar and you have a raster cell boundary > |90|, which does happen sometimes - then you need to adjust the input data). Roger
traceback()
8: stop("Geographical CRS given to non-conformant data")
7: `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">)
6: SpatialGrid(grid, proj4string)
5: initialize(value, ...)
4: initialize(value, ...)
3: new("SpatialGridDataFrame", SpatialGrid(grid, proj4string), data = data)
2: SpatialGridDataFrame(grid = grid, data = data.frame(df), proj4string = CRS(p4s))
1: readGDAL("RAD_NL21_PCP_NA_200707200000.tif")
sessionInfo()
R version 2.5.1 (2007-06-27) i386-pc-mingw32 locale: LC_COLLATE=Dutch_Netherlands.1252;LC_CTYPE=Dutch_Netherlands.1252;LC_MONETARY=Dutch_Netherlands.1252;LC_NUMERIC=C;LC_TIME=Dutch_Netherlands.1252 attached base packages: [1] "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base" other attached packages: rgdal sp "0.5-14" "0.9-14" Do you have any idea what caused this error and how I can circumvent it? Many thanks in advance, Paul Michael Sumner schreef:
Note that readGDAL can read the subdataset directly (escaping the quotes
appropriately):
d <- readGDAL("HDF5:\"RAD_NL21_PCP_NA_200707200000.h5\"://image1/image_data")
It also seems that your gdal_translate to TIFF output has the extension
.img, but then you attempt to read from a tif.
Just a thought.
Cheers, Mike.
==============Original message text===============
On Fri, 10 Aug 2007 19:09:00 +1000 Roger Bivand wrote:
Paul,
I have added a CC to the R-sig-geo list, since I myself have no experiece
of using this format in practice. If you are not a member, please join at:
https://stat.ethz.ch/mailman/listinfo/r-sig-geo so that you will see any replies. Did you try gdalinfo on the output file? The error message suggests that the output file has a geographical coordinate reference system, with coordinate values outside the feasible range. From the HDF5 description, however, it ought to be projected, right? Do you need to make more use of the arguments taken by gdalinfo to give the output file the correct projected CRS? Roger On Fri, 10 Aug 2007, Paul Hiemstra wrote: Dear Roger, I tried the route through a GTiff before I sent you the e-mail but somehow it did not work. When I use gdalinfo on one of the .h5 files I get this: PS E:\LSO data\Regenbui 20-7\regenradar_rasters\regenradar_hdf5> gdalinfo RAD_NL21_PCP_NA_200707200000.h5 Driver: HDF5/Hierarchical Data Format Release 5 Files: RAD_NL21_PCP_NA_200707200000.h5 Size is 512, 512 Coordinate System is `' Metadata: map_projection:projection_indication=Y map_projection:projection_name=STEREOGRAPHIC map_projection:projection_proj4_params=+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.388 +b=6356.912 +x_0=0 +y_0=0 geographic: geo_number_rows=256 geographic: geo_number_columns=256 geographic: geo_pixel_size_x=2.500000 geographic: geo_pixel_size_y=-2.500000 geographic: geo_par_pixel=X,Y geographic: geo_dim_pixel=KM,KM geographic: geo_column_offset=0.000000 geographic: geo_row_offset=1490.906006 geographic: geo_pixel_def=LU geographic: geo_product_corners=0.000000 49.769001 0.000000 55.296001 9.743000 54.818001 8.337000 49.373001 calibration: calibration_flag=Y calibration: calibration_formulas=GEO = 0.5 * PV + -31.5 calibration: calibration_missing_data=255 calibration: calibration_out_of_image=255 statistics: stat_min_value=-1.500000 statistics: stat_max_value=45.000000 image1: image_product_name=RAD_NL21_PCP_H0.8_NA image1: image_geo_parameter=REFLECTIVITY_[DBZ] image1: image_size=65536d image1: image_bytes_per_pixel=1d overview: product_group_name=RAD_NL21_PCP_NA overview: products_missing=N/A overview: hdftag_version_number=3.5 overview: number_image_groups=1 overview: number_satellite_groups=0d overview: number_radar_groups=2d overview: number_station_groups=0 overview: product_datetime_start=20-JUL-2007;00:00:11.000 overview: product_datetime_end=20-JUL-2007;00:00:11.000 overview: number_lightning_groups=0d overview: number_classification_groups=0d overview: number_grid_groups=0d overview: number_point_groups=0d overview: number_vector_groups=0d radar1: radar_name=De_Bilt radar1: radar_location=5.179000 52.103001 radar2: radar_name=Den_Helder radar2: radar_location=4.790000 52.955002 Subdatasets: SUBDATASET_0_NAME=HDF5:"RAD_NL21_PCP_NA_200707200000.h5"://image1/image_data SUBDATASET_0_DESC=[256x256] //image1/image_data (8-bit unsigned character) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 512.0) Upper Right ( 512.0, 0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0) Then I use gdal_translate: PS E:\LSO data\Regenbui 20-7\regenradar_rasters\regenradar_hdf5> gdal_translate -of GTiff HDF5:"RAD_NL21_PCP_NA_200707200000.h5"://image1/image_data RAD_NL21_PCP_NA_200707200000.img Then I try to read the file using rgdal produces the following error: test = readGDAL("regenradar_hdf5/RAD_NL21_PCP_NA_200707200000.tif") regenradar_hdf5/RAD_NL21_PCP_NA_200707200000.tif has GDAL driver GTiff and has 256 rows and 256 columns Closing GDAL dataset handle 0x0130b898... destroyed ... done. Error in `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) : Geographical CRS given to non-conformant data Do you have an idea how to solve this? kind regards, Paul Roger Bivand schreef: Dear Paul, There is information in the README.windows file in the package, this should display it: file.show(system.file("README.windows", package="rgdal")) If you do not want to build against the FWTools binary, then just install FWTools, and use the gdal_translate command line tool to convert your data to a format that the standard rgdal windows binary can read, such as gTiff. You'll need to play around a bit to find out what the data are called. On Linux, this is easier, because you can build GDAL against HDF5 before installing rgdal, and there it "just works". Hope this helps, Roger On Thu, 9 Aug 2007, Paul Hiemstra wrote: Dear mr. Bivand My name is Paul Hiemstra, a PhD student of Edzer. I am trying to read HDF5 files through rgdal. If I use GDALinfo() on one of the files it produces the following error: Error in .local(.Object, ...) : GDAL Error 4: `RAD_NL21_PCP_NA_200707200000.h5' not recognised as a supported file format. On the mailing list (http://finzi.psych.upenn.edu/R/Rhelp02a/archive/92270.html) you >> > metioned using the dll files from FWTools to get support for HDF5 into rgdal. How can this be done? I also tried the hdf5 package on CRAN, but that does not seem to work, returning an empty object. many thanks, Paul
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