OSGeo / PROJ-data

Repository for proj datum grids (for use by PROJ 7 or later)
Other
70 stars 32 forks source link

Add vertical grid file for Belgian geoid model hBG18 #77

Closed AlexBass05 closed 2 years ago

AlexBass05 commented 2 years ago

The grid file allows converting from ellipsoidal height in ETRS89 (EPSG:4937) to orthometric height in Ostend height (EPSG:5710) using be_ign_hBG18.tif

The GeoTIFF file was created with the attached build_hGB18.sh script.

Corresponding PR in PROJ: https://github.com/OSGeo/PROJ/pull/3044

Created as a draft PR as it requires an EPSG database update to include this change request: https://epsg.org/closed-change-request/browse/id/2022.007

rouault commented 2 years ago

The validation script displays

The following warnings were found (the file will probably be read correctly by PROJ, but corrective action may be required):
 - Geotransform with a positive pixel height, that is a south-up image, is supported, but a unusual formulation

This warning could be avoided by first doing a "gdalwarp hBG18.dat temp.tif" that would produce a "north-up" raster, and then using vertoffset_grid_to_gtiff on it to add the tags.

AlexBass05 commented 2 years ago

The validation script displays

The following warnings were found (the file will probably be read correctly by PROJ, but corrective action may be required):
 - Geotransform with a positive pixel height, that is a south-up image, is supported, but a unusual formulation

This warning could be avoided by first doing a "gdalwarp hBG18.dat temp.tif" that would produce a "north-up" raster, and then using vertoffset_grid_to_gtiff on it to add the tags.

Thanks for pointing this out. I regenerated the GeoTIFF (applying your suggestion) and ran regenerate_index_html.py again.

rouault commented 2 years ago

Thanks for pointing this out.

ok, actually I missed that the original raster didn't have square pixels, and that gdalwarp produces square pixels by default, so you have to invoke it with more arguments, something

like "gdalwarp hBG18.dat temp.tif -overwrite -tr 0.015 0.010 -te 0.9925 52.505 7.0075 48.495" you should check that the output of gdalinfo on the original and target grid results in the same reported extent, resolution and dimensions. and with the gdalsrsinfo utility on a few sample points on the original and resulting grid that you get the same values like "gdallocationinfo -geoloc be_ign_hBG18.tif 4 50"

AlexBass05 commented 2 years ago

ok, actually I missed that the original raster didn't have square pixels, and that gdalwarp produces square pixels by default, so you have to invoke it with more arguments,

Once more: Thanks for the advice.

The tiff has been updated and I added some testing output to the generation script. The index has been regenerated as well.