Add NGA EGM 2008 1x1 grid_alternative entry

[msmitherdc]

In order to account for the missing Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz file in the CDN, we need to add a grid_alternative entry to point to the 1' geotiff file.

('Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz','us_nga_egm2008_1.tif','egm08_1.gtx','GTiff','geoid_like',0,NULL,'https://cdn.proj.org/us_nga_egm2008_1.tif',1,1,NULL)

1' file can be found at https://grid-partner-share.s3.amazonaws.com/egm2008/us_nga_egm2008_1.tif

[rouault]

Hi @msmitherdc Who created / what was the creation process for that file ?

This file is 427 MB large. I don't think we want to include that directly in this github repo, nor in the PROJ-data tarball. We can upload it to cdn.proj.org, but should probably have a hint somewhere to keep track of such files that are CDN only. Probably a cdn_only_files.geojson stored in this repository that can be ingested by https://github.com/OSGeo/PROJ-data/blob/master/regenerate_index_html.py to include it in the generated https://github.com/OSGeo/PROJ-data/blob/master/files.geojson .
Otherwise people using projsync wouldn't get it.
And a mention of it in https://github.com/OSGeo/PROJ-data/blob/master/us_nga/us_nga_README.txt

CC @hobu @kbevers @jjimenezshaw

[rouault]

Looking at the metadata of the 1' file, I see:

  AREA_OR_POINT=Point
  target_crs_epsg_code=3855
  TIFFTAG_COPYRIGHT=Derived from work by USA NGA
  TIFFTAG_IMAGEDESCRIPTION=WGS 84 (EPSG::4979) to EGM2008 height (EPSG::3855)
  TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL

compared to the existing us_nga_egm08_25.tif:

  area_of_use=World
  AREA_OR_POINT=Point
  target_crs_epsg_code=3855
  TIFFTAG_COPYRIGHT=Derived from work by NGA. Public Domain
  TIFFTAG_DATETIME=2019:12:27 00:00:00
  TIFFTAG_IMAGEDESCRIPTION=WGS 84 (EPSG:4979) to EGM2008 height (EPSG:3855). Converted from egm08_25.gtx (last modified at 2018/10/08)
  TYPE=VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL

So the optional "area_of_use=World" is missing. And the mention of the license (Public Domain?) in the TIFFTAG_COPYRIGHT too

[msmitherdc]

I believe it was created from the .gtx file and yes, I just meant it to go into the cdn. Its referenced here https://earth-info.nga.mil/php/download.php?file=egm-08interpolation but not included. This is a file we've had for a while and not sure if it came direct from NGA or not. I can work with NGA to see if I can get something official

[msmitherdc]

Hmm looks like the file came from https://s3-eu-west-1.amazonaws.com/download.agisoft.com/gtg/us_nga_egm2008_1.tif https://www.agisoft.com/downloads/geoids/ and it says Converted from [USA NGA](http://earth-info.nga.mil/) data under Public Domain license.

[rouault]

ok, we could potentially recreate the file from the data and hints at https://geographiclib.sourceforge.io/C++/doc/geoid.html#geoidinst since @cffk is a trusted source for that topic!

At first glance, agisoft file seems OK:

 gdallocationinfo -wgs84 us_nga_egm2008_1.tif 2 49
Report:
  Location: (10920P,2460L)
  Band 1:
    Value: 44.6771507263184

compared to GeographicLib PGM dataset which has a scale of 0.003 and offset of -108:

 gdallocationinfo -wgs84 geoids/egm2008-1.pgm 2 49
Report:
  Location: (120P,2460L)
  Band 1:
    Value: 50892

$ python -c "print 50892 * 0.003 - 108"
44.676000000000016

[msmitherdc]

I’ve reached out to NGA to get any updates direct from them.

[rouault]

Procedure to generate a GeoTIFF file ranging from ~ -180 to ~ 180 deg longitude from the PGM file which ranges from ~ 0 to ~ 360:

gdal_translate geoids/egm2008-1.pgm west.vrt -srcwin 10800 0 10800 10801 -a_ullr -180.00833333333333333  90.008333333333333333  -0.0083333333333333  -90.008333333333333
gdal_translate geoids/egm2008-1.pgm east.vrt -srcwin 0 0 10800 10801
gdalbuildvrt out.vrt west.vrt east.vrt
gdal_translate out.vrt us_nga_egm08_1min.tif -co compress=deflate -co tiled=yes -co num_threads=12 -a_offset -108 -a_scale 0.003 -co predictor=2 -mo AREA_OR_POINT=Point -mo target_crs_epsg_code=3855 -a_srs EPSG:4979 # some missing metadata to add

echo 2 49 0 |  ~/proj/proj/build_cmake/bin/cct +proj=vgridshift +grids=us_nga_egm08_1min.tif +multiplier=1
  2.0000000000   49.0000000000      44.6760           inf

The resulting file is actually 78 MB large due to using integer storage !

[jjimenezshaw]

The resulting file is actually 78 MB large due to using integer storage !

What?! the file us_nga_egm08_25.tif (with 2.5 min step) is 80 MB!

I remember that 3 years ago the file was available in the NGA webpage, and I downloaded it. Unfortunately I deleted it when I changed my laptop. I did compile the fortran code to produce an ascii file from the binary provided.

If the final output is really only 78 MB we could include it. I was expecting something more than 6 times the 2.5' version, that would be around 500 MB!
If it is finally ~500 MB, definitely we should do something special for it. (BTW, is it really needed at that 1 minute accuracy?)

If we can compress it that much, should we do the same to existing files? What are the drawbacks?

Having it in the list of alternatives but not in the standard package could be a problem with the option ONLY_BEST?

[rouault]

the file us_nga_egm08_25.tif (with 2.5 min step) is 80 MB!

yes, because if uses Float32 encoding, whereas the PGM file used for the 1' file uses scaled UInt16 data type. And the PREDICTOR=3 mode that applies to floating-point numbers is less efficient than PREDICTOR=2 on UInt16 data. The 2.5'' could also be reduced in a similar manner. The scaled encoding "only" allows a 0.003 m precision on the data, but the GeographicLib page mentions "However, the resulting quantization error (up to 1.5 mm) is typically smaller than the linear interpolation errors."

Having it in the list of alternatives but not in the standard package could be a problem with the option ONLY_BEST?

Well, if the 1' grid is listed in grid_alternatives but not available locally or PROJ_NETWORK is not enabled, then yes when enabling ONLY_BEST, it would logically complain

[msmitherdc]

Here is the official file from NGA s3://grid-partner-share/egm2008/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz
and a presentation about EGM s3://grid-partner-share/egm2008/NPavlis&al_EGU2008.ppt

And I'm sorry that the file is 800Mb.

[jjimenezshaw]

Thanks @msmitherdc
How can I download that s3 url?

[rouault]

https://grid-partner-share.s3.amazonaws.com/egm2008/Und_min1x1_egm2008_isw%3D82_WGS84_TideFree.gz but there's indeed a permission issue. Anyway I don't think that file is directly usable as it must contain some "harmonics coefficients" or things like that, and that's GeographicLib that has done the heavy work of converting that to a grid

[msmitherdc]

Permissions should be fixed

[jjimenezshaw]

If I understood the format correctly, this script should convert to gtx (actually the data is pretty similar)

import struct

def main(src, dst):
    cols = 21600
    rows = 10801
    with open(dst, 'wb') as gtx, open(src, 'rb') as infile:
        gtx.write(struct.pack('>ddddii', -90, 0, 1.0/60, 1.0/60, rows, cols))

        rs = []
        for count in range(rows):
            infile.read(4)
            line = infile.read(cols*4)
            if not line:
                return 1
            rs.append(line)
            infile.read(4)

        for row in reversed(rs):
            gtx.write(row)

    return 0

if __name__ == "__main__":
    src = 'Und_min1x1_egm2008_isw=82_WGS84_TideFree'
    dst = 'egm2008.gtx'
    exit(main(src, dst))

The file does not contain the harmonics, but the geoid undulations directly. Before and after each row it has 4 bites with the integer value 86400.

[rouault]

If I understood the format correctly, this script should convert to gtx (actually the data is pretty similar)

cool! I assumed it wasn't a grid... and do you get the same values as https://grid-partner-share.s3.amazonaws.com/egm2008/us_nga_egm2008_1.tif / https://s3-eu-west-1.amazonaws.com/download.agisoft.com/gtg/us_nga_egm2008_1.tif ?