search

OSGeo / PROJ-data

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

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

Closed and-marsh closed 4 years ago

and-marsh commented 4 years ago

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

I have read https://github.com/OSGeo/PROJ/issues/1315 and https://github.com/OSGeo/PROJ/issues/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: https://github.com/OSGeo/PROJ/pull/2250

rouault commented 4 years ago

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

2) 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

3) 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

4) You need to update the copyright_and_licenses.csv file as well

and-marsh commented 4 years ago

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

and-marsh commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

Pixel size in meters: 1000x1000 Number of pixels:

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:

rouault commented 4 years ago

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: https://github.com/OSGeo/gdal/commit/309046522ca3ed728b61f00db6948880a073ce2e

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.

and-marsh commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

@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 commented 4 years ago

@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 commented 4 years ago

Index regenerated

SilverSloth commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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.