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.

Sign up for GitHub

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

Jump to bottom

Plans to support Geodetic data Grid eXchange Format (GGXF) #4370

Open
gmoney1729 opened this issue Jan 7, 2025 · 10 comments
Open

Plans to support Geodetic data Grid eXchange Format (GGXF) #4370

gmoney1729 opened this issue Jan 7, 2025 · 10 comments
Labels
feature request

Comments

@gmoney1729
Copy link

Are there plans to support the newer OGC standard GGXF for geoid models? NGS recently published one in this format under their NAPGD2022 geopotential datum release.

Thanks.
@rouault
Copy link
Member

rouault commented Jan 7, 2025

Are there plans to support the newer OGC standard GGXF for geoid models?

Is there funding for that?

In the mean time they still do produce under their good old .b format which is handled by GDAL and could be used to produced a GeoTIFF file:

$ gdalinfo /vsicurl/https://alpha.ngs.noaa.gov/NAPGD2022/data/geoid2022/SGEOID2022.NA.N.v1.a.b
Driver: NOAA_B/NOAA GEOCON/NADCON5 .b format
Files: /vsicurl/https://alpha.ngs.noaa.gov/NAPGD2022/data/geoid2022/SGEOID2022.NA.N.v1.a.b
Size is 10801, 5401
Coordinate System is:
GEOGCRS["Unspecified geographic CRS",
    DATUM["Unspecified datum based on GRS80 ellipsoid",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    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]]]
Data axis to CRS axis mapping: 2,1
Origin = (169.991666666666674,90.008333333333340)
Pixel Size = (0.016666666666667,-0.016666666666667)
Corner Coordinates:
Upper Left  ( 169.9916667,  90.0083333) (169d59'30.00"E, 90d 0'30.00"N)
Lower Left  ( 169.9916667,  -0.0083333) (169d59'30.00"E,  0d 0'30.00"S)
Upper Right (     350.008,      90.008) (350d 0'30.00"E, 90d 0'30.00"N)
Lower Right (     350.008,      -0.008) (350d 0'30.00"E,  0d 0'30.00"S)
Center      (     260.000,      45.000) (260d 0' 0.00"E, 45d 0' 0.00"N)
Band 1 Block=10801x1 Type=Float32, ColorInterp=Undefined

@gmoney1729
Copy link
Author

Not that I'm aware of at the moment but will keep my eyes peeled. I noticed that I could read the .b and .bin files for the 3 regional geoids except the NA geoid in .bin format where gdal throws an error 4. Not an issue, just an observation. Thanks.

@gmoney1729 gmoney1729 closed this as not planned Won't fix, can't repro, duplicate, stale Jan 7, 2025
@rouault
Copy link
Member

rouault commented Jan 7, 2025

keeping open, as this will (unfortunately) be something someone will have to address at some point

@rouault rouault reopened this Jan 7, 2025
@busstoptaktik
Copy link
Member

But note that the X in GGXF is for "eXchange". GGXF is by definition NOT for use in actual systems, but for exchange between system specific formats. Hence, PROJ should never directly support the use of the GGXF format, but might consider supporting the functional model behind it

@rouault
Copy link
Member

rouault commented Jan 8, 2025

but might consider supporting the functional model behind it

To clarify what I in mind: If/when we will "support" it, we'll convert to our internal format (GeoTIFF). So a set of scripts that would live in https://github.com/OSGeo/PROJ-data/tree/master/grid_tools

Regarding GEOID2022 specificity, the work is not just handling the format, but potentially extending our proj=gridshift operation to support the time-based velocity component that is added on top of the static geoid.

@jjimenezshaw
Copy link
Contributor

jjimenezshaw commented Jan 8, 2025 •
edited
Loading

But note that the X in GGXF is for "eXchange". GGXF is by definition NOT for use in actual systems, but for exchange between system specific formats. Hence, PROJ should never directly support the use of the GGXF format, but might consider supporting the functional model behind it

How is then NGS publishing GEOID2022.v1.a.ggxf?
Should GDAL read GGXF files and convert to GTG?

edit: this message had a race condition with @rouault ;)

@rouault
Copy link
Member

rouault commented Jan 8, 2025

Should GDAL read GGXF files and convert to GTG?

GDAL can somehow read them using its netCDF driver. What is missing is the georeferencing, since bright minds in the CRS SWG have decided, for good reasons I don't want to understand, and against my recommendation, not to use the netCDF CF conventions for georeferencing. At the end of the day, a GDAL driver is probably not needed. Some custom python script using h5py or the GDAL Python API will probably do the job.

@jjimenezshaw
Copy link
Contributor

GDAL can somehow read them using its netCDF driver

How could I do that? If I try with GDAL 3.4.1 I get this

gdal_translate -if netcdf GEOID2022.v1.a.ggxf -of gtiff geoid2022.tif
ERROR 4: `GEOID2022.v1.a.ggxf' not recognized as a supported file format

@rouault
Copy link
Member

rouault commented Jan 19, 2025

If I try with GDAL 3.4.1 I get this

maybe your libnetcdf lacks HDF5 support ?
Anyway don't expect gdal_translate to produce a GeoTIFF. There needs to be custom scripting around to translate the custom georeferencing that was picked up by GGXF.

@rouault
Copy link
Member

rouault commented Jan 20, 2025 

$ docker run --rm -it -v $HOME:$HOME ghcr.io/osgeo/gdal gdalinfo -if netcdf /home/even/Téléchargements/GEOID2022.v1.a.ggxf
Driver: netCDF/Network Common Data Format
Files: /home/even/Téléchargements/GEOID2022.v1.a.ggxf
Size is 512, 512
Metadata:
  /DGEOID2022 Uncertainty/dgeoid2022-uncertainty/NC_GLOBAL#affineCoeffs={-90,0.25,0,0,0,0.25}
  /DGEOID2022 Uncertainty/dgeoid2022-uncertainty/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /DGEOID2022 Uncertainty/NC_GLOBAL#gridParameters=geoidVelocityUncertainty
  /DGEOID2022 Uncertainty/NC_GLOBAL#interpolationMethod=biquadratic
  /DGEOID2022/dgeoid2022/NC_GLOBAL#affineCoeffs={-90,0.25,0,0,0,0.25}
  /DGEOID2022/dgeoid2022/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /DGEOID2022/NC_GLOBAL#gridParameters=geoidVelocity
  /DGEOID2022/NC_GLOBAL#interpolationMethod=biquadratic
  /SGEOID2022 Uncertainty/NC_GLOBAL#gridParameters=geoidHeightUncertainty
  /SGEOID2022 Uncertainty/NC_GLOBAL#interpolationMethod=biquadratic
  /SGEOID2022 Uncertainty/sgeoid2022-uncertainty/NC_GLOBAL#affineCoeffs={-90,0.08333333333333,0,0,0,0.08333333333333}
  /SGEOID2022 Uncertainty/sgeoid2022-uncertainty/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/American Samoa/NC_GLOBAL#affineCoeffs={-22,0.016666666666667,0,172,0,0.016666666666667}
  /SGEOID2022/American Samoa/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/American Samoa/NC_GLOBAL#gridPriority=2
  /SGEOID2022/Guam and Northern Mariana Islands/NC_GLOBAL#affineCoeffs={5,0.016666666666667,0,130,0,0.016666666666667}
  /SGEOID2022/Guam and Northern Mariana Islands/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/Guam and Northern Mariana Islands/NC_GLOBAL#gridPriority=2
  /SGEOID2022/NC_GLOBAL#gridParameters=geoidHeight
  /SGEOID2022/NC_GLOBAL#interpolationMethod=biquadratic
  /SGEOID2022/North America/NC_GLOBAL#affineCoeffs={0,0.016666666666667,0,170,0,0.016666666666667}
  /SGEOID2022/North America/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/North America/NC_GLOBAL#gridPriority=1
  NC_GLOBAL#city=Silver Spring
  NC_GLOBAL#content=producerDefinedContent
  NC_GLOBAL#Conventions=GGXF-1.0, ACDD-1.3
  NC_GLOBAL#country=United States of America
  NC_GLOBAL#deliveryPoint=1315 East West Hwy
  NC_GLOBAL#extent_description=North American-Pacific Region
  NC_GLOBAL#geospatial_lat_max=90
  NC_GLOBAL#geospatial_lat_min=-22
  NC_GLOBAL#geospatial_lon_max=350
  NC_GLOBAL#geospatial_lon_min=130
  NC_GLOBAL#institution=National Geodetic Survey, National Oceanic and Atmospheric Administration.
  NC_GLOBAL#interpolationCrsWkt=GEOGCRS["ITRF2020",
  DATUM["International Terrestrial Reference Frame 2020",
      ELLIPSOID["GRS 1980",6378137.0,298.2572221,LENGTHUNIT["metre",1]]],
  CS[ellipsoidal,2],
  AXIS["Geodetic latitude (Lat)",north],
  AXIS["Geodetic longitude (Lon)",east],
  ANGLEUNIT["degree",0.0174532925199433]]

  NC_GLOBAL#operationAccuracy=0.01
  NC_GLOBAL#parameters.0.parameterName=geoidHeight
  NC_GLOBAL#parameters.0.sourceCrsAxis=2
  NC_GLOBAL#parameters.0.unitName=metre
  NC_GLOBAL#parameters.0.unitSiRatio=1
  NC_GLOBAL#parameters.1.parameterName=geoidHeightUncertainty
  NC_GLOBAL#parameters.1.sourceCrsAxis=2
  NC_GLOBAL#parameters.1.uncertaintyMeasure=1SE
  NC_GLOBAL#parameters.1.unitName=metre
  NC_GLOBAL#parameters.1.unitSiRatio=1
  NC_GLOBAL#parameters.2.parameterName=geoidVelocity
  NC_GLOBAL#parameters.2.sourceCrsAxis=2
  NC_GLOBAL#parameters.2.unitName=m/yr
  NC_GLOBAL#parameters.2.unitSiRatio=3.16887651727315e-08
  NC_GLOBAL#parameters.3.parameterName=geoidVelocityUncertainty
  NC_GLOBAL#parameters.3.sourceCrsAxis=2
  NC_GLOBAL#parameters.3.uncertaintyMeasure=1SE
  NC_GLOBAL#parameters.3.unitName=m/yr
  NC_GLOBAL#parameters.3.unitSiRatio=3.16887651727315e-08
  NC_GLOBAL#parameters.count=4
  NC_GLOBAL#postalCode=20910
  NC_GLOBAL#producerDefinedContentCitation=
  NC_GLOBAL#product_version=v1.a
  NC_GLOBAL#publisher_url=
  NC_GLOBAL#sourceCrsWkt=GEOGCRS["ITRF2020",
  DATUM["International Terrestrial Reference Frame 2020",
      ELLIPSOID["GRS 1980",6378137.0,298.2572221,LENGTHUNIT["metre",1]]],
  CS[ellipsoidal,3],
  AXIS["Geodetic latitude (Lat)",north,
    ANGLEUNIT["degree",0.0174532925199433]],
  AXIS["Geodetic longitude (Lon)",east,
    ANGLEUNIT["degree",0.0174532925199433]],
  AXIS["Ellipsoidal height (h)",up,LENGTHUNIT["metre",1]]]

  NC_GLOBAL#source_file=GEOID2022.v1.a.ggxf
  NC_GLOBAL#summary=Gravimetric geoid model for the North American-Pacific region.

This file contains the GEOID2022 model which consists of 4 layers/parameters:
  
  geoidHeight: static geoid undulation (m)
  geoidHeightUncertainty: uncertainty in geoidHeight (m)
  geoidVelocity: change in geoid undulation surface in m/yr
  geoidVelocityUncertainty: uncertainty in the change in geoidVelocity (m/yr)

To get the geoid undulation value at a specific epoch, t, use the following equation:

  geoidHeight(t) = geoidHeight + geoidVelocity*(t-t0)

Where t0 is the reference epoch of the model, 2020.0.

Uncertainty at a specific epoch, t, can be calculated as follows:

  geoidHeightUncertainty(t) = sqrt(geoidHeightUncertainty**2 + (geoidVelocityUncertainty*(t-t0))**2)
  

  NC_GLOBAL#targetCrsWkt=VERTCRS["NAPGD2022",
  VDATUM["North American-Pacific Geopotential Datum of 2022"],
  CS[vertical,1],
  AXIS["Gravity-related height (H)",up],
  LENGTHUNIT["metre",1]]

  NC_GLOBAL#title=GEOID2022
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":"/SGEOID2022/North America/geoidHeight"
  SUBDATASET_1_DESC=[5401x10801] /SGEOID2022/North America/geoidHeight (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":"/SGEOID2022/Guam and Northern Mariana Islands/geoidHeight"
  SUBDATASET_2_DESC=[1021x2401] /SGEOID2022/Guam and Northern Mariana Islands/geoidHeight (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":"/SGEOID2022/American Samoa/geoidHeight"
  SUBDATASET_3_DESC=[1321x1261] /SGEOID2022/American Samoa/geoidHeight (32-bit floating-point)
  SUBDATASET_4_NAME=NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":"/SGEOID2022 Uncertainty/sgeoid2022-uncertainty/geoidHeightUncertainty"
  SUBDATASET_4_DESC=[2161x4321] /SGEOID2022 Uncertainty/sgeoid2022-uncertainty/geoidHeightUncertainty (32-bit floating-point)
  SUBDATASET_5_NAME=NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":/DGEOID2022/dgeoid2022/geoidVelocity
  SUBDATASET_5_DESC=[721x1441] /DGEOID2022/dgeoid2022/geoidVelocity (32-bit floating-point)
  SUBDATASET_6_NAME=NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":"/DGEOID2022 Uncertainty/dgeoid2022-uncertainty/geoidVelocityUncertainty"
  SUBDATASET_6_DESC=[721x1441] /DGEOID2022 Uncertainty/dgeoid2022-uncertainty/geoidVelocityUncertainty (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

$ docker run --rm -it -v $HOME:$HOME ghcr.io/osgeo/gdal gdalinfo NETCDF:"/home/even/Téléchargements/GEOID2022.v1.a.ggxf":"/SGEOID2022/North America/geoidHeight"
Driver: netCDF/Network Common Data Format
Files: /home/even/Téléchargements/GEOID2022.v1.a.ggxf
Size is 10801, 5401
Metadata:
  /DGEOID2022 Uncertainty/dgeoid2022-uncertainty/NC_GLOBAL#affineCoeffs={-90,0.25,0,0,0,0.25}
  /DGEOID2022 Uncertainty/dgeoid2022-uncertainty/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /DGEOID2022 Uncertainty/NC_GLOBAL#gridParameters=geoidVelocityUncertainty
  /DGEOID2022 Uncertainty/NC_GLOBAL#interpolationMethod=biquadratic
  /DGEOID2022/dgeoid2022/NC_GLOBAL#affineCoeffs={-90,0.25,0,0,0,0.25}
  /DGEOID2022/dgeoid2022/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /DGEOID2022/NC_GLOBAL#gridParameters=geoidVelocity
  /DGEOID2022/NC_GLOBAL#interpolationMethod=biquadratic
  /SGEOID2022 Uncertainty/NC_GLOBAL#gridParameters=geoidHeightUncertainty
  /SGEOID2022 Uncertainty/NC_GLOBAL#interpolationMethod=biquadratic
  /SGEOID2022 Uncertainty/sgeoid2022-uncertainty/NC_GLOBAL#affineCoeffs={-90,0.08333333333333,0,0,0,0.08333333333333}
  /SGEOID2022 Uncertainty/sgeoid2022-uncertainty/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/American Samoa/NC_GLOBAL#affineCoeffs={-22,0.016666666666667,0,172,0,0.016666666666667}
  /SGEOID2022/American Samoa/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/American Samoa/NC_GLOBAL#gridPriority=2
  /SGEOID2022/Guam and Northern Mariana Islands/NC_GLOBAL#affineCoeffs={5,0.016666666666667,0,130,0,0.016666666666667}
  /SGEOID2022/Guam and Northern Mariana Islands/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/Guam and Northern Mariana Islands/NC_GLOBAL#gridPriority=2
  /SGEOID2022/NC_GLOBAL#gridParameters=geoidHeight
  /SGEOID2022/NC_GLOBAL#interpolationMethod=biquadratic
  /SGEOID2022/North America/NC_GLOBAL#affineCoeffs={0,0.016666666666667,0,170,0,0.016666666666667}
  /SGEOID2022/North America/NC_GLOBAL#comment=grid starts in the bottom left (southwest) corner and works across (east) and up (north)

  /SGEOID2022/North America/NC_GLOBAL#gridPriority=1
  NC_GLOBAL#city=Silver Spring
  NC_GLOBAL#content=producerDefinedContent
  NC_GLOBAL#Conventions=GGXF-1.0, ACDD-1.3
  NC_GLOBAL#country=United States of America
  NC_GLOBAL#deliveryPoint=1315 East West Hwy
  NC_GLOBAL#extent_description=North American-Pacific Region
  NC_GLOBAL#geospatial_lat_max=90
  NC_GLOBAL#geospatial_lat_min=-22
  NC_GLOBAL#geospatial_lon_max=350
  NC_GLOBAL#geospatial_lon_min=130
  NC_GLOBAL#institution=National Geodetic Survey, National Oceanic and Atmospheric Administration.
  NC_GLOBAL#interpolationCrsWkt=GEOGCRS["ITRF2020",
  DATUM["International Terrestrial Reference Frame 2020",
      ELLIPSOID["GRS 1980",6378137.0,298.2572221,LENGTHUNIT["metre",1]]],
  CS[ellipsoidal,2],
  AXIS["Geodetic latitude (Lat)",north],
  AXIS["Geodetic longitude (Lon)",east],
  ANGLEUNIT["degree",0.0174532925199433]]

  NC_GLOBAL#operationAccuracy=0.01
  NC_GLOBAL#parameters.0.parameterName=geoidHeight
  NC_GLOBAL#parameters.0.sourceCrsAxis=2
  NC_GLOBAL#parameters.0.unitName=metre
  NC_GLOBAL#parameters.0.unitSiRatio=1
  NC_GLOBAL#parameters.1.parameterName=geoidHeightUncertainty
  NC_GLOBAL#parameters.1.sourceCrsAxis=2
  NC_GLOBAL#parameters.1.uncertaintyMeasure=1SE
  NC_GLOBAL#parameters.1.unitName=metre
  NC_GLOBAL#parameters.1.unitSiRatio=1
  NC_GLOBAL#parameters.2.parameterName=geoidVelocity
  NC_GLOBAL#parameters.2.sourceCrsAxis=2
  NC_GLOBAL#parameters.2.unitName=m/yr
  NC_GLOBAL#parameters.2.unitSiRatio=3.16887651727315e-08
  NC_GLOBAL#parameters.3.parameterName=geoidVelocityUncertainty
  NC_GLOBAL#parameters.3.sourceCrsAxis=2
  NC_GLOBAL#parameters.3.uncertaintyMeasure=1SE
  NC_GLOBAL#parameters.3.unitName=m/yr
  NC_GLOBAL#parameters.3.unitSiRatio=3.16887651727315e-08
  NC_GLOBAL#parameters.count=4
  NC_GLOBAL#postalCode=20910
  NC_GLOBAL#producerDefinedContentCitation=
  NC_GLOBAL#product_version=v1.a
  NC_GLOBAL#publisher_url=
  NC_GLOBAL#sourceCrsWkt=GEOGCRS["ITRF2020",
  DATUM["International Terrestrial Reference Frame 2020",
      ELLIPSOID["GRS 1980",6378137.0,298.2572221,LENGTHUNIT["metre",1]]],
  CS[ellipsoidal,3],
  AXIS["Geodetic latitude (Lat)",north,
    ANGLEUNIT["degree",0.0174532925199433]],
  AXIS["Geodetic longitude (Lon)",east,
    ANGLEUNIT["degree",0.0174532925199433]],
  AXIS["Ellipsoidal height (h)",up,LENGTHUNIT["metre",1]]]

  NC_GLOBAL#source_file=GEOID2022.v1.a.ggxf
  NC_GLOBAL#summary=Gravimetric geoid model for the North American-Pacific region.

This file contains the GEOID2022 model which consists of 4 layers/parameters:
  
  geoidHeight: static geoid undulation (m)
  geoidHeightUncertainty: uncertainty in geoidHeight (m)
  geoidVelocity: change in geoid undulation surface in m/yr
  geoidVelocityUncertainty: uncertainty in the change in geoidVelocity (m/yr)

To get the geoid undulation value at a specific epoch, t, use the following equation:

  geoidHeight(t) = geoidHeight + geoidVelocity*(t-t0)

Where t0 is the reference epoch of the model, 2020.0.

Uncertainty at a specific epoch, t, can be calculated as follows:

  geoidHeightUncertainty(t) = sqrt(geoidHeightUncertainty**2 + (geoidVelocityUncertainty*(t-t0))**2)
  

  NC_GLOBAL#targetCrsWkt=VERTCRS["NAPGD2022",
  VDATUM["North American-Pacific Geopotential Datum of 2022"],
  CS[vertical,1],
  AXIS["Gravity-related height (H)",up],
  LENGTHUNIT["metre",1]]

  NC_GLOBAL#title=GEOID2022
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 5401.0)
Upper Right (10801.0,    0.0)
Lower Right (10801.0, 5401.0)
Center      ( 5400.5, 2700.5)
Band 1 Block=10801x1 Type=Float32, ColorInterp=Undefined
  NoData Value=9.96921e+36
  Metadata:
    NETCDF_VARNAME=geoidHeight

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
feature request
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@rouault @busstoptaktik @jjimenezshaw @gmoney1729