[R-sig-Geo] reading HDF5 files
Roger Bivand
Roger.Bivand at nhh.no
Tue Aug 14 11:48:38 CEST 2007
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
More information about the R-sig-Geo
mailing list