Incorrect conversion from ETRS89 to OSGB36+ODN with OSGM15 geoid

[hernando]

I'm trying to get accurate conversions from ERTS89 (EPSG:4937) to OSG36B+ODN (EPSG:7405) using a geoid model grid that I derived myself from publicly available data (the grid is only 1.3 MB so I could share attach it if necessary). I will give details on how I created the grid later, but basically I followed the instructions of the the grid_tools from PROJ-data, the grid doesn't seem to be the problem anyway.

I'm using the test coordinates provided in the OS developer pack

Example of problem

I modified the WKT string for 7405 to include PROJ4_GRIDS in the vertical datum, then I used this string as input to cs2cs to test a few transformations.

$ EPSG_7405=$(projinfo -q -o WKT1:GDAL EPSG:7405)
$ EPSG_7405_ext=$(echo $EPSG_7405 | sed -e 's/Newlyn",2005/Newlyn",2005, EXTENSION["PROJ4_GRIDS", "uk_os_etrs89_osgm15.tif"]/')
$ export PROJ_LIB=paths to proj.db and the grid files

$ echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 "$EPSG_7405" -f "%.6f"
91492.146319    11318.804552 100.000000

$ echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 "$EPSG_7405_ext" -f "%.6f"
91495.585557    11326.276036 46.517545

Problem description

The reference output for this transformation is 91492.146 11318.804 46.519.
The first transformation gets the horizontal part right, but the vertical component is not transformed. This is expected, as the file referred by proj.db for OSGM15 is not found (the file is OSTN15_OSGM15_GB.txt and I haven't found it anywhere, where does it come from?).
In the second transformation the vertical transformation is almost perfect, but the horizontal transformation is ignoring the horizontal grid shift for ETRS89 to OSGB36.

The grid file I made was obtained by downloading a GeoTIFF from Agisoft which contains the same data as the OSTN15_OSGM15_DataFile.txt file from the OS Developers Pack but in a more convenient format. I used QGIS to resample the data and change the source CRS from the projected OSGB36 to ETRS89 (so it's geographic) and then I used vertoffset_grid_to_gtiff.py from the Docker image osgeo/gdal:alpine-normal-latest to convert it to the standard GeoTIFF expected by Proj.

Expected Output

$ echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 "$EPSG_7405" -f "%.6f"
91492.146319    11318.804552 100.000000
$ echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 "$EPSG_7405_ext" -f "%.6f"
91495.146319   11326.804552 46.517545

Environment Information

• PROJ master branch commit 5e3aba3
• Operation System Information: Ubuntu 18.04

Installation method

• Compiled from source

[hernando]

I've also tried modifying the SQL table grid_alternatives adding this to customizations.sql:

INSERT INTO grid_alternatives(original_grid_name,
                              proj_grid_name,
                              old_proj_grid_name,
                              proj_grid_format,
                              proj_method,
                              inverse_direction,
                              package_name,
                              url, direct_download, open_license, directory)
VALUES ('OSTN15_OSGM15_GB.txt','uk_os_ETRS89toOSGM15.tif','OSTN15_OSGM15_GB.txt','GTiff','geoid_like',0,NULL,NULL,0,0,NULL);

After doing this, the output of projinfo -s EP SG:4937 -t EPSG:7405 -o PROJ shows this:

Candidate operations found: 5                                                                                           
Note: using '--spatial-test intersects' would bring more results (11)                                                   
-------------------------------------                                                                                   
Operation No. 1:                                                                                                        
                                                                                                                        
unknown id, ETRS89 to ODN height (2) + Inverse of OSGB 1936 to ETRS89 (2) + British National Grid, 0.038 m, UK - Great Britain mainland onshore                                                                                                 
                                                                                                                        
PROJ string:                                                                                                            
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=uk_os_ETRS89toOSGM15.tif +multiplier=1 +step +inv +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif +step +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy

But

echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 EPSG:7405 -f "%.6f"

still outputs 91492.146319 11318.804552100.000000.

I'm trying to feed that PROJ string to proj like this: proj $(projinfo -s EPSG:4937 -t EPSG:7405 -o PROJ -q), but I get:

Rel. 7.1.0, August 1st, 2020
<proj>: 
can't initialize operations that take non-angular input coordinates
program abnormally terminated

[rouault]

Instead of hacking the EPSG:7405 definition, you probably just need to add an entry in data/sql/grid_alternatives.sql to map the EPSG "official" name 'OSTN15_OSGM15_GB.txt' to your grid. And if it works and the grid can be licensed under the ones we can accept, it would be more than welcome!

I used QGIS to resample the data and change the source CRS from the projected OSGB36 to ETRS89 (so it's geographic)

Change from projected to geographic + datum shifts will probably cause some "artefacts" that can explain you will not get exactly the reference output, but that's probably acceptable. Care should be done to having the horizontal shift OSGB36->ETRS89 right when reprojecting the grid.
It would be nice if you contribute the grid that you provide full instructions to reproduce what you did. My own preference would go to something based on GDAL : command line use of gdalwarp and/or some GDAL Python script

[rouault]

After doing this, the output of projinfo -s EP SG:4937 -t EPSG:7405 -o PROJ shows this:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=uk_os_ETRS89toOSGM15.tif +multiplier=1 +step +inv +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif +step +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy

That looks very good !

But echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 EPSG:7405 -f "%.6f"
still outputs 91492.146319 11318.804552100.000000.

Hum, strange. Try setting PROJ_DEBUG=3 to have verbose details

I'm trying to feed that PROJ string to proj like this

Pipelines can't work with the proj ancient utility. Use cct instead

echo 49.92226393730 -6.29977752014 100.000 | cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=uk_os_ETRS89toOSGM15.tif +multiplier=1 +step +inv +proj=hgridshift +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif +step +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy

[hernando]

Nice, it worked with cct!
Now the question is why cs2cs doesn't work. I've already tried PROJ_DEBUG, but the output is so long that I get lost. I've attached a file with it.

Regarding the grid, I will try to produce it using GDAL. Assuming that I'm resampling the original data from a projected CRS into a geographic one and that the interpolation is bilinear, having mm errors is totally acceptable in my opinion. Strictly speaking, the question is how much does my grid differ from the original data at the benchmark positions of the ODN network, but I have neither the knowledge nor the time to do that analysis. On the other hand, making sure that the resampling is using the correct horizontal grid is a very good point, because that could introduce some additional mm of error that is totally avoidable and I don't know what QGIS did in that regard.

[hernando]

I naively thought I could attach a file... Here you are: https://gist.github.com/hernando/af03ca08de428baaab1399ca73a248fd#file-cs2c2_debug_output

[rouault]

The log is interesting since at the end it mentions "Using coordinate operation Inverse of Transformation from ODN height to ETRS89 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of OSGB 1936 to ETRS89 (2) + British National Grid", so it ignores the vertical grid, but earlier in it, it successfully opens the vertical grid. I'd suspect maybe something related to the area of use of the transformation.
You could check with "projinfo -s EP SG:4937 -t EPSG:7405 --spatial-test intersects" if the declared BBOX[] of the transformation that uses both the horizontal & vertical grids includes the point you want to transform

Assuming that I'm resampling the original data from a projected CRS into a geographic one and that the interpolation is bilinear,

"gdalwarp -r bilinear -s_srs EPSG:27700 -t_srs EPSG:4258 ..." should do the job (assuming the projected grid is in EPSG:27700), and use PROJ underneath to apply the horizontal shift grid.

[hernando]

You nailed it! The coordinates 49.92226393730 -6.29977752014 are in fact a bit west to the box reported for that transformation: BBOX[49.93,-7.06,58.71,1.8]. That point is on-land in the Isles of Scilly (never heard of before), and they are not very far from the west coast of Wales (https://www.google.com/maps/place/49%C2%B055'20.2%22N+6%C2%B017'59.2%22W). I surprised me that they could be out of the bounding box, but I did a bit of research and discovered that they have their own area of use (https://epsg.io/2802-area) and vertical datum (despite the geoid model for the transformation is the same) and grid_transformations.sql contains a vertical transformation for them.

So, if I do echo 49.92226393730 -6.29977752014 100.000 | cs2cs EPSG:4937 EPSG:27700+5749, the transformation is correct.

I guess you can close this ticket, thanks for the help. I will submit the PRs for the grid if I can.

[rouault]

Note: pending pull request to add that grid in OSGeo/PROJ-data#25