[R-sig-Geo] reading HDF5 files

Thu Aug 16 09:19:21 CEST 2007

On Tue, 14 Aug 2007, Tim Keitt wrote:

> On 8/14/07, Paul Hiemstra <p.hiemstra at geo.uu.nl> 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
>>> 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
>
> This seems a case where failing gracefully would save a lot of
> trouble.

I beg to disagree. The conversion from HDF5 to GTiff in gdal_transform (so 
outside R and rgdal) erroneously turned a NULL coordinate reference system 
in the HDF5 file to a WGS84 longlat. Setting the correct coordinate 
reference system using -a_srs to the one cited in the HDF5 *metadata* by 
hand corrected the problem, so that the dataset is now known as 
Stereographic, as it should be. The only check that failed was that data 
specifically declared as being in geographical coordinates had coordinate 
values outside the feasible range.

Why the HDF5 file used the metadata rather than the coordinate reference 
system to store the CRS is beyond me, as is the behaviour of 
gdal_translate in erroneously imputing geographical coordinates, but the 
work-around does work.

> Perhaps a warning and a NULL CRS? Then things can be patched
> up in R. I am not a fan of the defensive programming style prevalent
> in R. (Although I fully recognize that it is a matter of taste and
> many people would disagree with me.)

There are ways round here too - open as a GDALDataset, and build the grid 
description and coordinate reference system string by hand, which would be 
equivalent to using gdal_translate to do it. I don't think R is very 
defensive as such, but just a warning in this case would be unhelpful for 
the user.

Roger

>
> THK
>
>
>> 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
>>>>>>
>>>>>>
>>>>>>
>>>>
>>>>
>>>
>>>
>>
>>
>> --
>> Drs. Paul Hiemstra
>> Department of Physical Geography
>> Faculty of Geosciences
>> University of Utrecht
>> Heidelberglaan 2
>> P.O. Box 80.115
>> 3508 TC Utrecht
>> Phone:  +31302535773
>> Fax:    +31302531145
>> http://intamap.geo.uu.nl/~paul
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> 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. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no