Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add uk_os_OSGM15_GB.tif geoid file for ODN height, UK #25

Merged
merged 5 commits into from
Jun 18, 2020

Conversation

and-marsh
Copy link
Contributor

@and-marsh and-marsh commented Jun 2, 2020

Hi,
I decided to add lacked geoid file for EPSG(5701)ODN height, UK

I have read OSGeo/PROJ#1315 and OSGeo/PROJ#2173 issues and tried to build uk_os_OSGM15_GB.tif geoid file from source file OSTN15_OSGM15_DataFile.txt obtained from Ordnance Survey developer pack.

I attached shell script to generate output file. There I used custom projected CRS ETRS89 / British National Grid since the points of the source grid are configured in this system in accordane to the OSGM15_Notice_of_release_for_developers.pdf attached to the developer pack

1 The record number is a sequential number starting at 1 for the origin point (0,0) and finishing at 876 951 for the north-east corner (700 000, 1 250 000).
2 ETRS89, National Grid projection, grid intersection easting coordinate in metres.
3 ETRS89, National Grid projection, grid intersection northing coordinate in metres.
4 The shift in eastings, at the intersection, between ETRS89 and OSGB36 National Grid, that is: ETRS89 east + OSTN15 east shift = OSGB36 National Grid easting.
5 The shift in northings, at the intersection, between ETRS89 and OSGB36 National Grid, that is: ETRS89 north + OSTN15 north shift = OSGB36 National Grid northing.
6 The height of the Geoid above the ETRS89 ellipsoid, in metres, at the intersection, that is: ETRS89 height – OSGM15 Geoid height = orthometric height above mean sea level.
7 The Geoid datum flag is a number representing the local height datum or area of applicability of the transformation. See the table below for details of the datum flag references.

So, I also conducted tests on data obtained from the developer pack for transformation EPSG:4258 -> EPSG:27700+5701:

» tail -n +2 OSTN15_OSGM15_TestInput_ETRStoOSGB.txt | awk 'BEGIN {FS=","; OFS=" "} {print $3, $2, $4}' > OSTN15_OSGM15_TestInput_ETRStoOSGB_prepared.txt
» tail -n +2 OSTN15_OSGM15_TestOutput_ETRStoOSGB.txt | awk 'BEGIN {FS=","; OFS=" "} {print $2, $3, $4}' > OSTN15_OSGM15_TestOutput_ETRStoOSGB_prepared.txt
» cat OSTN15_OSGM15_TestInput_ETRStoOSGB_prepared.txt | ./gdaltransform -s_srs EPSG:4258 -t_srs EPSG:27700+5701 > OSTN15_OSGM15_TestOutput_ETRStoOSGB_actual.txt
» paste OSTN15_OSGM15_TestOutput_ETRStoOSGB_prepared.txt OSTN15_OSGM15_TestOutput_ETRStoOSGB_actual.txt | awk '{print $1 - $4, $2 - $5, $3 - $6}' > OSTN15_OSGM15_TestOutput_ETRStoOSGB_delta.txt

Results of tests - difference of coordinates of points obtained after transformation and provided by the developer pack, all values in meters:

dX dY dZ
-0.000318582 -0.000551564 -53.481
0.000633833 -0.000919592 -0.000389089
7.2703e-05 -0.00739838 -0.00018952
-0.00138396 -0.000875295 0.000456284
0.000221762 -0.00131693 -0.000206267
-0.000542716 -0.000539408 0.00110022
0.001324 -0.000243979 -0.000238674
-0.00233302 0.00188051 0.00013368
0.000737046 -0.000864641 0.000803671
-0.00105092 0.00201758 0.00307603
-0.000426731 0.000731854 0.00101824
-0.0027537 -0.00136221 -0.000218854
0.000512831 -0.00324241 0.000102695
-0.00134748 -0.0014907 0.000287787
0.00327066 -0.00707292 0.000363898
-0.00180256 0.000853918 -0.000280263
-0.000257473 -0.000671562 -0.00128095
-0.000260841 -0.000666039 -0.00132467
0.0019246 -0.00129102 0.0003028
-0.0008565 -0.00101684 6.87085e-05
-0.000778786 -0.00105426 0.000261991
-0.000254567 -0.0017251 0.00155019
-0.000627372 -0.000768617 -0.000869601
-0.000847119 -0.00156796 0.000940164
0.00222652 -0.00138332 -0.000440808
-0.00085512 -0.001071 0.000150201
0.00375297 -0.00123378 -0.000203778
-0.000434079 -0.00125638 -0.000157201
-0.000387423 -0.000895149 0.000140499
-0.00100261 -0.000757105 -0.00150954
0.00156496 -0.00262059 -57.988
0.000466969 -0.00224644 -56.672
0.000189872 -0.000596542 -0.000403209
-0.000686879 -0.000627394 -0.000776666
-0.000251475 -0.00116233 -52.044
-0.000101695 -0.00143951 -53.555
-0.000386164 -0.00091319 -55.367
-0.00118104 -0.00050829 -48.951
-0.000338373 -0.00061752 -48.901
-0.000162162 -0.00131532 -50.701

Height of some points weren't transformed since these points are outside ODN height area of use - UK - Great Britain mainland onshore. Accuracy of height transformation is up to 1.5 mm.

Related PR: OSGeo/PROJ#2250

@rouault
Copy link
Member

rouault commented Jun 2, 2020

  1. The quality checks run by the CI raise the following warnings:
The following warnings were found (the file will probably be read correctly by PROJ, but corrective action may be required):
 - This file uses a RasterTypeGeoKey = PixelIsArea convention. While this will work properly with PROJ, a georeferencing using RasterTypeGeoKey = PixelIsPoint is more common.
 - Missing TYPE metadata item
 - GDAL area_of_use metadata item is missing. Typically used to capture plain text information about where the grid applies
 - TIFF tag ImageDescription is missing. Typically used to capture plain text information about what the grid does
 - TIFF tag Copyright is missing. Typically used to capture information on the source and licensing terms
 - TIFF tag DateTime is missing. Typically used to capture when the grid has been generated/converted
 - Image is uncompressed. Could potentially benefit from compression
 - Given image dimension, tiling could be appropriate

You could solve most of them by running the grid_tools/vertoffset_grid_to_gtiff.py script of this repo as a post-processing step on your warped TIFF file. For the first point regarding PixelIsPoint, adding -mo AREA_OR_POINT=Point to gdalwarp should fix it

  1. Regarding the grid resolution, gdalwarp automatically selects square pixels. But given the latitude of UK and a target geographic CRS, this is probably not appropriate. Specifying an explicit -tr where the resolution in longitude is roughly the resolution in latitude * cos(lat) would probably be a bit more appropriate. Actually I would start by determinining the res_latitude ( = res_northing / 111319 ), and then compute res_longitude from it

  2. Due to warp changing the shape of the grid, you'll likely get values at 0 at the edges of the reprojected raster. Two alternatives: either subset the raster afterwards to remove them, or add a -dstnodata -32768 so that those values are at a tagged invalid value

  3. You need to update the copyright_and_licenses.csv file as well

@and-marsh
Copy link
Contributor Author

Thank you! I'll try to follow your remarks.

@and-marsh and-marsh marked this pull request as draft June 2, 2020 09:45
@and-marsh
Copy link
Contributor Author

@rouault, I have some questions to you and hope you will not mind helping me
Here is gdalinfo of new grid with almost all fixes applied:

Driver: GTiff/GeoTIFF
Files: test-2-processed.tif
Size is 1443, 1266
Coordinate System is:
GEOGCRS["ETRS89",
    DATUM["European Terrestrial Reference System 1989",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4258]]
Data axis to CRS axis mapping: 2,1
Origin = (-9.400175143506939,61.135879203514222)
Pixel Size = (0.008982829987676,-0.008983192446932)
Metadata:
  area_of_use=Great Britain mainland onshore
  AREA_OR_POINT=Point
  target_crs_epsg_code=5701
  TIFFTAG_COPYRIGHT=Derived from work by Ordnance Survey. The 2-Clause BSD License https://opensource.org/licenses/bsd-license.php
  TIFFTAG_DATETIME=2020:06:02 00:00:00
  TIFFTAG_IMAGEDESCRIPTION=ETRS89 (EPSG:4937) to ODN height (EPSG:5701). Converted from OSTN15_OSGM15_DataFile.txt
  TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -9.4001751,  61.1358792) (  9d24' 0.63"W, 61d 8' 9.17"N)
Lower Left  (  -9.4001751,  49.7631576) (  9d24' 0.63"W, 49d45'47.37"N)
Upper Right (   3.5620485,  61.1358792) (  3d33'43.37"E, 61d 8' 9.17"N)
Lower Right (   3.5620485,  49.7631576) (  3d33'43.37"E, 49d45'47.37"N)
Center      (  -2.9190633,  55.4495184) (  2d55' 8.63"W, 55d26'58.27"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = geoid_undulation
  NoData Value=-32768
  Unit Type: metre

So, I have added all necessary metadata and tried to change pixel size by setting -tr flag. I don't quite understand how exactly did you get 111319 in resolution latitude formulа. See my calculations:

>>> res_lat=1000/111319
>>> res_long=res_lat*math.cos(res_lat)
>>> res_lat
0.008983192446931791
>>> res_long
0.008982829987675677

and then

./gdalwarp -r bilinear -dstnodata "-32768" -tr 0.008982829987675677 0.008983192446931791 -s_srs "${ETRS89_BRITISH_NAT_GRID}" -t_srs EPSG:4258 osgm15_temp osgm15_temp.tif
./gdal_translate -mo "AREA_OR_POINT=Point" osgm15_temp.tif ${TARGET_FILE}

I have no enough experience and knowledge to estimate if new transformation is more correct, so, I'm very appreciate your help! Pleas, take a look at the new results of tests:

-0.000318582 -0.000551564 -53.481
0.000633833 -0.000919592 -0.000405646
7.2703e-05 -0.00739838 -5.59963e-05
-0.00138396 -0.000875295 0.000383306
0.000221762 -0.00131693 0.000678684
-0.000542716 -0.000539408 0.00100392
0.001324 -0.000243979 -6.31774e-05
-0.00233302 0.00188051 0.000131025
0.000737046 -0.000864641 0.000998632
-0.00105092 0.00201758 0.0019535
-0.000426731 0.000731854 0.000908441
-0.0027537 -0.00136221 -5.7632e-05
0.000512831 -0.00324241 3.45652e-05
-0.00134748 -0.0014907 0.000128022
0.00327066 -0.00707292 0.000843169
-0.00180256 0.000853918 2.98292e-05
-0.000257473 -0.000671562 -0.000982741
-0.000260841 -0.000666039 -0.00102905
0.0019246 -0.00129102 0.000273473
-0.0008565 -0.00101684 5.88329e-05
-0.000778786 -0.00105426 0.000332238
-0.000254567 -0.0017251 0.000929893
-0.000627372 -0.000768617 -0.000683017
-0.000847119 -0.00156796 0.000448668
0.00222652 -0.00138332 -0.000416919
-0.00085512 -0.001071 -0.000686977
0.00375297 -0.00123378 -0.00183045
-0.000434079 -0.00125638 0.000175505
-0.000387423 -0.000895149 0.000290185
-0.00100261 -0.000757105 -0.00118098
0.00156496 -0.00262059 -57.988
0.000466969 -0.00224644 -56.672
0.000189872 -0.000596542 -7.97021e-05
-0.000686879 -0.000627394 -0.000487888
-0.000251475 -0.00116233 -52.044
-0.000101695 -0.00143951 -53.555
-0.000386164 -0.00091319 -55.367
-0.00118104 -0.00050829 -48.951
-0.000338373 -0.00061752 -48.901
-0.000162162 -0.00131532 -50.701

If it's acceptable I will update PR and copyright_and_licenses.csv file.

Thank you!

@rouault
Copy link
Member

rouault commented Jun 2, 2020

how exactly did you get 111319 in resolution latitude formulа

the length in metres of one degree of longitude at the equator of the GRS80 ellipsoid = pi * 6378137 / 180

res_long=res_lat*math.cos(res_lat)

Should rather be res_long=res_lat*math.cos(min_latitude_of_grid_expressed_in_radian)

I have no enough experience and knowledge to estimate if new transformation is more correct

you could possibly compute the mean square root error on a few sample points and compare

@and-marsh
Copy link
Contributor Author

and-marsh commented Jun 3, 2020

Hi, again!
I took into account your comments and that's what i got:

>>> res_lat=1000/111319
>>> res_lat
0.008983192446931791

>>> min_lat=49.7631576
>>> res_long=res_lat*math.cos(math.radians(min_lat))
>>> res_long
0.00580268140457895

rmse=0.0007641565216395795

>>> mid_lat=55.4495184
>>> res_long=res_lat*math.cos(math.radians(mid_lat))
>>> res_long
0 .005094657066217041

rmse=0.000762624030848561

>>> max_lat=61.1358792
>>> res_long=res_lat*math.cos(math.radians(max_lat))
res_long=0.004336492995274525

rmse=0.0007727981480493811

RMSE for the very first geoid were:
rmse=0.0008830190950406347

So, i choose variant with middle latitude. New gdalinfo:

Driver: GTiff/GeoTIFF
Files: ./PROJ-data/uk_os/uk_os_OSGM15_GB.tif
Size is 2544, 1266
Coordinate System is:
GEOGCRS["ETRS89",
    DATUM["European Terrestrial Reference System 1989",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4258]]
Data axis to CRS axis mapping: 2,1
Origin = (-9.400175143506939,61.135879203514222)
Pixel Size = (0.005094657066217,-0.008983192446932)
Metadata:
  area_of_use=Great Britain mainland onshore
  AREA_OR_POINT=Point
  target_crs_epsg_code=5701
  TIFFTAG_COPYRIGHT=Derived from work by Ordnance Survey. The 2-Clause BSD License https://opensource.org/licenses/bsd-license.php
  TIFFTAG_DATETIME=2020:06:03 00:00:00
  TIFFTAG_IMAGEDESCRIPTION=ETRS89 (EPSG:4937) to ODN height (EPSG:5701). Converted from OSTN15_OSGM15_DataFile.txt
  TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -9.4001751,  61.1358792) (  9d24' 0.63"W, 61d 8' 9.17"N)
Lower Left  (  -9.4001751,  49.7631576) (  9d24' 0.63"W, 49d45'47.37"N)
Upper Right (   3.5606324,  61.1358792) (  3d33'38.28"E, 61d 8' 9.17"N)
Lower Right (   3.5606324,  49.7631576) (  3d33'38.28"E, 49d45'47.37"N)
Center      (  -2.9197714,  55.4495184) (  2d55'11.18"W, 55d26'58.27"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = geoid_undulation
  NoData Value=-32768
  Unit Type: metre

Also I added comment with adding metadata to the script.

@and-marsh and-marsh marked this pull request as ready for review June 3, 2020 08:53
@and-marsh
Copy link
Contributor Author

For RMSE calculation I used test OSTN15_OSGM15_TestInput_ETRStoOSGB.txt dataset without points that are outside area of use. RMSE was calculated only for Z coordinates with python:

mse = sklearn.metrics.mean_squared_error(actual_filtered, expected_filtered)
rmse = math.sqrt(mse)

@rouault
Copy link
Member

rouault commented Jun 3, 2020

What's the size (in pixels) and resolution of the input grid ?
Perhaps look at the maximum absolute error too ?
For consistency with other geoid grids, EPSG:4937 (ETRS89 3D) should also be used for the GeoTIFF information (what is reported in "Coordinate System is" by gdalinfo)

@and-marsh
Copy link
Contributor Author

and-marsh commented Jun 3, 2020

Pixel size in meters: 1000x1000
Number of pixels:

  • Northing - 1250
  • Easting - 700

There are used EPSG:4258 for Belfast and Malin grids, so I suppose it would be better to leave EPSG:4258 for this grid too.

Driver: GTiff/GeoTIFF
Files: ./PROJ-data/uk_os/uk_os_OSGM15_Belfast.tif
Size is 326, 376
Coordinate System is:
GEOGCRS["ETRS89",
    DATUM["European Terrestrial Reference System 1989",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4258]]

Moreover, vertoffset_grid_to_gtiff.py sets EPSG:4258 even if EPSG:4937 is used for input file

Absolute errors:

  • For minimal latitude: max_abs_err=0.00191889845089932
  • For middle latitude: max_abs_err=0.0019401348965004672
  • For maximal latitude: max_abs_err=0.0019515429794978445

@rouault
Copy link
Member

rouault commented Jun 3, 2020

There are used EPSG:4258 for Belfast and Malin grids

ah, not for me. It reports EPSG:4937. There's a related fix in GDAL 3.1: OSGeo/gdal@3090465

Number of pixels: Northing - 1250 Easting - 700

OK, so we have a consistent number of samples in the resampled product with 1266. But I obviously got the maths wrong for the longitude resolution: it should be res_lat / cos(lat) . We should aim to go closer to the original number of samples in the longitude direction: I'd say greater than 700 but not significantly larger than 1000 or so.

Source data is `OSTN15_OSGM15_DataFile.txt` file obtained from
OrdnanceSurvey developer pack (see README)
@and-marsh
Copy link
Contributor Author

and-marsh commented Jun 3, 2020

it should be res_lat / cos(lat)

Recalculated:

# RMSE = 0.0007570669375634179m
# MAX_ERROR = 0.0019630308326981094m

We should aim to go closer to the original number of samples in the longitude direction

Now I understand the logic! Check it out

Size is 818, 1266

ah, not for me. It reports EPSG:4937.

I changed -t_srs for gdalwarp to EPSG:4937, but can't check it. Pleas make sure of this.

@rouault
Copy link
Member

rouault commented Jun 3, 2020

Looking at EPSG:7711 metadata (ETRS89 to ODN height (2)), it mentions an accuracy of 8 mm. So here our RMSE is one order of magnitude below it (7.6e-4 m), and our max error is 2mm. @kbevers any opinion if given the above using a grid referenced to geographic coordinates compared to the original one referenced to projected coordinates is a reasonable compromise ?

@rouault
Copy link
Member

rouault commented Jun 16, 2020

@and-marsh Could you add a note in the README about the slight added error due to the reprojection (with mensions of the RMSE and max error) ? with that, I guess we should be good to go

@and-marsh
Copy link
Contributor Author

@rouault I have added some notes to the README. Can I ask you regenerate index file, since for now I don't have GDAL with Python bindings? Thanks a lot!

@rouault rouault added this to the 1.1.0 milestone Jun 18, 2020
@rouault rouault merged commit 0bbb5d1 into OSGeo:master Jun 18, 2020
rouault added a commit to OSGeo/PROJ that referenced this pull request Jun 18, 2020
@rouault
Copy link
Member

rouault commented Jun 18, 2020

Index regenerated

@and-marsh and-marsh deleted the pr-add-main-uk-geoid branch June 18, 2020 12:51
@SilverSloth
Copy link

Referring to the “Transformations_and_OSGM15_UserGuide.pdf“ contained in the developers pack, page 10, Table and note 6. I believe the developer pack provides “Geoid height”. Note 6 further explains that: ETRS height - OSGM15 Geoid height = orthometric height above mean sea level. So to obtain the vertical shift from ETRS89 (GRS80) to orthometric, the Geoid height must be multiplied by minus 1. The grid currently has positive instead of negative shifts.

@rouault
Copy link
Member

rouault commented Jul 1, 2020

The grid currently has positive instead of negative shifts.

yes, that's correct. Geoid models as used in PROJ (and in most other formats) must be the value of the geoid undulation, that is the ellipsoidal height for a orthometric altitude of 0.
The value in the uk_os_OSGM15_GB.tif file are consistent with the EGM96 geoid model:

$ gdallocationinfo -wgs84 uk_os_OSGM15_GB.tif 0 50
Report:
  Location: (593P,1239L)
  Band 1:
    Value: 44.4401626586914

$ gdallocationinfo -wgs84 us_nga_egm96_15.tif 0 50
Report:
  Location: (720P,160L)
  Band 1:
    Value: 45.0383186340332

@SilverSloth
Copy link

Oh, it’s a OSGM15 Geoid height model relative to ETRS89/GRS80 ellipsoid, not a vertical shift from ETRS89/GRS80 to OSGB36/ODN. Thanks for helping the penny drop. Not all vertical shifts are Geoid models, so would be helpful to know this one is.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants