reading HDF5 files
On Tue, 14 Aug 2007, Paul Hiemstra wrote:
Mike,
This is the result of using gdalinfo (FWTools) on the .tif file,
E:\LSO data\Regenbui 20-7\regenradar_rasters\regenradar_tif>gdalinfo
RAD_NL21_PCP_NA_200707200000.tif
Driver: GTiff/GeoTIFF
Files: RAD_NL21_PCP_NA_200707200000.tif
RAD_NL21_PCP_NA_200707200000.aux
Size is 256, 256
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.2572235630016,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Metadata:
AREA_OR_POINT=Area
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_product_corners=0.000000 49.769001 0.000000 55.296001
9.743000 54.818001 8.337000 49.373001
calibration=calibration_out_of_image=255
statistics=stat_max_value=45.000000
image1=image_bytes_per_pixel=1d
overview=number_vector_groups=0d
radar1=radar_location=5.179000 52.103001
radar2=radar_location=4.790000 52.955002
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 256.0)
Upper Right ( 256.0, 0.0)
Lower Right ( 256.0, 256.0)
Center ( 128.0, 128.0)
Band 1 Block=256x32 Type=Byte, ColorInterp=Gray
Metadata:
image_data=DISPLAY_ORIGIN=UL
The metadata records the correct srs, but the file has the wrong one - it is staring you in the face, isn't it? Try adding -a_srs "+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" (all in one line, no line breaks) to the gdal_translate call (untried). Roger
cheers, 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