Quantcast

[gdal-dev] Convert precipitation accumulations from HDF5 to GeoTiff

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[gdal-dev] Convert precipitation accumulations from HDF5 to GeoTiff

Paul Meems
I'm trying to use the precipitation accumulations from our national weather agency: https://data.knmi.nl/datasets/radar_tar_corr_accum_03h/2.0
The data is in HDF5 with two subsets.
I only need one subset and I need the data in WGS84.

I'm using GDAL v2.2 from GisInternals.com
I tried using GDALWarp but then I get an error about missing the transformation matrix.
Using Google I figured out I need to do GDALTranslate first to set the georeferenced bounds.

This is the result of GDALInfo:
Driver: HDF5/Hierarchical Data Format Release 5
Files: RAD_NL25_RAC_03H_201612041200.h5
Size is 512, 512
Coordinate System is `'
Metadata:
  geographic_geo_column_offset=0 
  geographic_geo_dim_pixel=KM,KM
  geographic_geo_number_columns=700 
  geographic_geo_number_rows=765 
  geographic_geo_par_pixel=X,Y
  geographic_geo_pixel_def=LU
  geographic_geo_pixel_size_x=1.0000013 
  geographic_geo_pixel_size_y=-1.0000055 
  geographic_geo_product_corners=0 49.362053 0 55.973602 10.856429 55.388977 9.0092793 48.895298 
  geographic_geo_row_offset=3649.9792 
  geographic_map_projection_projection_indication=Y
  geographic_map_projection_projection_name=STEREOGRAPHIC
  geographic_map_projection_projection_proj4_params=+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0
  image1_calibration_calibration_flag=Y
  image1_calibration_calibration_formulas=GEO = 0.010000 * PV + 0.000000
  image1_calibration_calibration_missing_data=0 
  image1_calibration_calibration_out_of_image=65535 
  image1_image_bytes_per_pixel=2 
  image1_image_data_CLASS=IMAGE
  image1_image_data_VERSION=1.2
  image1_image_geo_parameter=PRECIP_[MM]
  image1_image_number_accumulated_products=37 
  image1_image_product_name=RAD_NL25_RAC_H1.5_03H
  image1_image_size=535500 
  image1_statistics_stat_max_value=9.96 
  image1_statistics_stat_min_value=0 
  image2_calibration_calibration_flag=Y
  image2_calibration_calibration_formulas=GEO=0.001*PV+-32.768
  image2_calibration_calibration_missing_data=0 
  image2_calibration_calibration_out_of_image=65535 
  image2_image_bytes_per_pixel=2 
  image2_image_data_CLASS=IMAGE
  image2_image_data_VERSION=1.2
  image2_image_geo_parameter=ADJUSTMENT_FACTOR_[DB]
  image2_image_number_accumulated_products=37 
  image2_image_product_name=RAD_NL25_RAC_H1.5_03H
  image2_image_size=535500 
  image2_statistics_stat_max_value=0 
  image2_statistics_stat_min_value=0 
  overview_hdftag_version_number=3.6
  overview_number_image_groups=2 
  overview_number_radar_groups=1 
  overview_number_satellite_groups=0 
  overview_number_station_groups=0 
  overview_product_datetime_end=04-DEC-2016;12:00:00.000
  overview_product_datetime_start=04-DEC-2016;09:00:00.000
  overview_product_group_name=RAD_NL25_RAC_03H
  overview_products_missing=NA
  radar1_radar_adjustment=F=-0.00dB
  radar1_radar_location=5.17834 52.101681 
  radar1_radar_name=DeBilt
Subdatasets:
  SUBDATASET_1_NAME=HDF5:"RAD_NL25_RAC_03H_201612041200.h5"://image1/image_data
  SUBDATASET_1_DESC=[765x700] //image1/image_data (16-bit unsigned integer)
  SUBDATASET_2_NAME=HDF5:"RAD_NL25_RAC_03H_201612041200.h5"://image2/image_data
  SUBDATASET_2_DESC=[765x700] //image2/image_data (16-bit unsigned integer)
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)
                                                                                                                                             
This is my GDALTranslate command, the projectionstring and a_ullr come from the meta data:
gdal_translate -of GTiff -a_nodata 65535 -a_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" -a_ullr 0 55.973602 9.0092793 48.895298 "HDF5:\"RAD_NL25_RAC_03H_201612041200.h5\"://image1/image_data" "test-translate.tif"

This is my GDALWarp command:
gdalwarp -overwrite -s_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0" -t_srs EPSG:3857 -dstnodata 65535 -of GTiff "test-translate.tif" "test-warp.tif"

The warped tif is not where I expect it to be.
I expect a origin around 52 en 5 degrees (The Netherlands)
I'm getting: Origin = (169.559973821725920,89.530840368222059)

I've checked with the data provider and the projectionstring is correct.
What am I doing wrong?

Thanks,

Paul

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: [gdal-dev] Convert precipitation accumulations from HDF5 to GeoTiff

Rutger
Hey Paul,

If you assign the RD projection with 'a_srs', you should also assign the 'ullr' coordinates in that same projection, not in lat/lon degrees like you do now.


Regards,
Rutger
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Convert precipitation accumulations from HDF5 to GeoTiff

Paul Meems
Thanks Rutger,

You are completely right. I assumed the coordinates were in the same projection as the project string.

Menno van Scheers replied me directly and suggested I should use 0.0, -3650.000, 800.000, -4415.000 as ullr coordinates.
I tried it and now the GeoTiff is where I expect it to be.




Paul

Paul Meems 
Release manager, configuration manager
and forum moderator of MapWindow GIS.
www.mapwindow.org

Owner of MapWindow.nl - Support for
Dutch speaking users.
www.mapwindow.nl


The MapWindow GIS project has moved to GitHub!


Download the latest MapWinGIS mapping engine.

Download the latest MapWindow 5 open source desktop application.


2016-12-13 11:29 GMT+01:00 Rutger <[hidden email]>:
Hey Paul,

If you assign the RD projection with 'a_srs', you should also assign the
'ullr' coordinates in that same projection, not in lat/lon degrees like you
do now.


Regards,
Rutger



--
View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-Convert-precipitation-accumulations-from-HDF5-to-GeoTiff-tp5299661p5299666.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Convert precipitation accumulations from HDF5 to GeoTiff

Rutger
Hey,

Based on the metadata from the file you posted i would guess the x-coordinate of the right edge (lrx) should be about ~700km or ~700000m.

The column offset is 0, and there are 700 columns in the file which are slightly over 1km. How would you end up with 800?

You can calculate the other corner coordinates in a similar fashion.


Regards,
Rutger
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Convert precipitation accumulations from HDF5 to GeoTiff

Paul Meems
Thanks again Rutger for your observation.

The bounding box was sent to me and I was already very happy the tiff is located in The Netherlands ;)

I did a closer look at the meta data and with your explanation I come up with this 'formula':

The ul is 0.0, -1 * geographic_geo_row_offset
lr is geographic_geo_number_columns, -1 * (geographic_geo_number_rows + geographic_geo_row_offset)
Resulting in 0.0 3650 700, -4415

Thanks,

Paul

 
2016-12-13 14:17 GMT+01:00 Rutger <[hidden email]>:
Hey,

Based on the metadata from the file you posted i would guess the
x-coordinate of the right edge (lrx) should be about ~700km or ~700000m.

The column offset is 0, and there are 700 columns in the file which are
slightly over 1km. How would you end up with 800?

You can calculate the other corner coordinates in a similar fashion.


Regards,
Rutger




--
View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-Convert-precipitation-accumulations-from-HDF5-to-GeoTiff-tp5299661p5299678.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev



_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Convert precipitation accumulations from HDF5 to GeoTiff

Rutger
Hey,

Yes, that looks good. Although instead of using the -1, i think its meant to be the actual y-resolution as provided in the metadata, but i'm never 100% sure on those things. Its tempting to ignore it in this case, since its so close to 1, but imagine what will happen if they increase the resolution after an upgrade. If you read the number of rows from the data, and use a hard coded resolution, the geographical extent would change instead of the resolution.

# geographic_geo_pixel_size_x
yres = -1.0000052 # km
xres = 1.0000013

# geographic_geo_column_offset
xoffset = 0
yoffset = 3649.9802

# geographic_geo_number_columns
ncol = 700
nrow = 765

ulx = xoffset * xres
uly = yoffset * yres
lrx = ulx + ncol * xres
lry = uly + nrow * yres

Which gives:
(0.0, -3649.9991798970395, 700.0009100000001, -4415.003157897039)

You can also calculate them by transforming the given lat/lon corner coordinates to the proj4 projection string given in the metadata. If i do that as well, i get:
(0.0, -3649.9991119177544, 700.000903671186, -4415.003881997636)

Since both should be similar (there could be some rounding in the provided corner coordinates), its a good way to check if you did it right.

The last thing to think about it what the coordinates are refering to, the outer edge of the corner pixel, or the center of the corner pixel. Given "geographic_geo_pixel_def=LU" in the metadata, i would guess it the actual outer edge (which is GDAL-style as well).

Here is a Jupyter Notebook which i used (but nbviewer is down right now, check it later):
http://nbviewer.jupyter.org/gist/RutgerK/70fe9d76d72b0fbfdeebda0568753828

And i have to apologize to Even for not using OSR's coordinate conversion in that example. ;)


Regards,
Rutger


 
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Convert precipitation accumulations from HDF5 to GeoTiff

Rutger
Sorry, what i said makes no sense at all. You shouldn't multiply the offset with the y-resolution, although you should take the sign, like you did. Its a bit strange that they provide the offset as a positive value.


Regards,
Rutger
Loading...