Add NGA EGM 2008 1x1 grid_alternative entry #111
Comments
|
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 . |
|
Looking at the metadata of the 1' file, I see: compared to the existing us_nga_egm08_25.tif: So the optional "area_of_use=World" is missing. And the mention of the license (Public Domain?) in the TIFFTAG_COPYRIGHT too |
|
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 |
|
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 |
|
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: compared to GeographicLib PGM dataset which has a scale of 0.003 and offset of -108: |
|
I’ve reached out to NGA to get any updates direct from them. |
|
Procedure to generate a GeoTIFF file ranging from ~ -180 to ~ 180 deg longitude from the PGM file which ranges from ~ 0 to ~ 360: The resulting file is actually 78 MB large due to using integer storage ! |
What?! the file 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 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 |
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."
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 |
|
Here is the official file from NGA s3://grid-partner-share/egm2008/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz And I'm sorry that the file is 800Mb. |
|
Thanks @msmitherdc |
|
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 |
|
Permissions should be fixed |
|
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. |
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 ? |
|
Those tif files that you mention @rouault have 21601 columns, while the original has 21600. In the file
And some stats:
Apart from that, comparing with QGIS they have the same values (in the common columns). (I'm always unsure about center or corner of the pixel. Also in the python script I did to convert to gtx. Anybody with more knowledge should check that, given the description I pasted above) |
This data is "pixel is point". I.e., the data give the heights at the corners of the cells (at lat + lon = integer multiples of 1' in this case). |
|
Is there some agreement that the 1.5 mm max difference due to quantization if using a UInt16 data type + offset/scale is acceptable ? |
|
For me 1.5 mm is fine. I asked a fried from the university and agrees. The model is worse than that. |
Taking a few "random" points, when doing bilinear interpolation, the difference between both versions is about within +/- 1.5 mm quantization error... It might be that I'm just out of luck, and it might be possible that in some areas, there is a significant difference... @msmitherdc Did you find that the 1' model was leading to more accurate results than the 2.5' one ? |
|
actually, my issue was that reading a EGM08 vertical datum file via PDAL was returning this message |
ah, then it is something else... This error is a bit surprising as PROJ shouldn't try to fetch from the CDN a grid it doesn't know in grid_alternatives. I can't reproduce with: It tries to access a local Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz file, but doesn't try to get it from the CDN |
|
the geotiff keys on the file that triggers the lookup look like |
This should be reported as a PDAL issue with precise instructions. It might be a PROJ issue, or a GDAL issue or a PDAL issue, but that's unclear at that point. |
|
My hunch is it is a libgeotiff issue, since we're just handing it all off to that in this case. |
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
The text was updated successfully, but these errors were encountered: