Developer Guide
This page serves to aid developers in correctly reading and implementing the component models of NAPGD2022.
Reading GGXF
NGS uses Geodetic data Grid eXchange Format (GGXF) as the file format for distribution of NAPGD2022 grids. GGXF uses NetCDF for its carrier encoding and thus may be read with appropriate NetCDF4 tools.
The following Python example illustrates how to explore and extract grid data from a GGXF file using the NetCDF4 library.
#Import Python NetCDF4
from netCDF4 import Dataset
import numpy as np
#Load the GGXF as a NetCDF object
ggxf_file = "path/to/ggxfs/GEOID2022.v1.ggxf"
rootgrp = Dataset(ggxf_file, mode="r")
#Subgroups by be explored by calling the group attribute, which gives a dicitonary of nested groups
'''
> rootgrp.groups
{'SGEOID2022': <class 'netCDF4._netCDF4.Group'>
group /SGEOID2022:
interpolationMethod: biquadratic
gridParameters: ['geoidHeight', 'geoidHeightUncertainty']
dimensions(sizes):
variables(dimensions):
groups: NA, GC, AS,
'DGEOID2022': <class 'netCDF4._netCDF4.Group'>
group /DGEOID2022:
interpolationMethod: biquadratic
gridParameters: ['geoidVelocity', 'geoidVelocityUncertainty']
dimensions(sizes):
variables(dimensions):
groups: GL}
'''
#The static geoid is under SGEOID2022 while the dynamic geoid is under DGEOID2022.
#DGEOID2022 is a single global grid, while SGEOID2022 is broken into three regions (NA, AS, and GC).
#In this example, we'll focus on the North America (NA) grid
region_path = 'SGEOID2022/NA"
#The affine coefficients and grid parameters may be extracted for a particular region
affineCoeffs = rootgrp[region_path].getncattr("affineCoeffs")
nlat = rootgrp[region_path ].dimensions['iNodeCount'].size
nlon = rootgrp[region_path ].dimensions['jNodeCount'].size
#The latitudes and longitudes of the grid corners may be constructed from the affine coefficients
lat = affineCoeffs[0] + np.arange(nlat)*affineCoeffs[1]
lon = affineCoeffs[3] + np.arange(nlon)*affineCoeffs[5]
#Finally, the grid data may be extracted using slicing
geoidHeight = rootgrp[region_path +"geoidHeight"][:]
#This returns a 2D NumPy array that may be plotted or interpolated with lat and lon as independent variables
#This extraction can be condensed into a single path
geoidHeight = rootgrp["SGEOID2022/NA/geoidHeight"][:]
A complete map of the NAPGD2022 attributes in the GGXF file structure is given below:
Time-dependent geoid undulation
GEOID2022.v1.ggxf
SGEOID2022
Unit: meter
Regions: NA, AS, GC
Resolution: 1’
Attributes: geoidHeight (geoid undulation at 2020.0)
geoidHeightUncertainty (1-sigma geoid uncertainty)
geoidQuasigeoidSeparation (geoid undulation minus quasigeoid)
DGEOID2022
Unit: meter per year
Regions: GL
Resolution: 15’
Attributes: geoidVelocity (geoid rate of change)
geoidVelocityUncertainty (1-sigma geoid rate uncertainty)
Time-dependent deflection of the vertical
DEFLEC2022.v1.ggxf
SDEFLEC2022
Unit: arcsecond
Regions: NA, AS, GC
Resolution: 1’
Attributes: deviationEast (prime vertical deflection, eta)
deviationNorth (meridional deflection, xi)
deviationEastUncertainty (1-sigma eta uncertainty)
deviationNorthUncertainty (1-sigma xi uncertainty)
DDEFLEC2022
Unit: arcsecond per year
Regions: GL
Resolution: 15’
Attributes: deviationEastVelocity (rate of change in eta)
deviationNorthVelocity (rate of change in xi)
Topographic surface orthometric height
DEM2022.v1.ggxf
DEM2022
Unit: meter
Regions: NA, AS, GC
Resolution: 1’
Attributes: elevation (orthometric height of topographic surface)
Refined Bouguer anomaly and terrain correction grids to predict gravity
GRAV2022.v1.ggxf
RBA
Unit: milligal
Regions: NA, AS, GC
Resolution: 1’
Attributes: gravityAnomaly (refined Bouguer anomaly)
gravityAnomalyUncertainty (1-sigma refined Bouguer anomaly uncertainty)
TC
Unit: milligal
Regions: NA, AS, GC
Resolution: 0.25’
Attributes: terrainCorrection (terrain correction)
For convenience, here is a complete list of paths to variables in the GGXF files
SGEOID2022/NA/geoidHeight
SGEOID2022/NA/geoidHeightUncertainty
SGEOID2022/NA/geoidQuasigeoidSeparation
SGEOID2022/AS/geoidHeight
SGEOID2022/AS/geoidHeightUncertainty
SGEOID2022/AS/geoidQuasigeoidSeparation
SGEOID2022/GC/geoidHeight
SGEOID2022/GC/geoidHeightUncertainty
SGEOID2022/GC/geoidQuasigeoidSeparation
DGEOID2022/GL/geoidVelocity
DGEOID2022/GL/geoidVelocityUncertainty
SDEFLEC2022/NA/deviationEast
SDEFLEC2022/NA/deviationNorth
SDEFLEC2022/NA/deviationEastUncertainty
SDEFLEC2022/NA/deviationNorthUncertainty
SDEFLEC2022/AS/deviationEast
SDEFLEC2022/AS/deviationNorth
SDEFLEC2022/AS/deviationEastUncertainty
SDEFLEC2022/AS/deviationNorthUncertainty
SDEFLEC2022/GC/deviationEast
SDEFLEC2022/GC/deviationNorth
SDEFLEC2022/GC/deviationEastUncertainty
SDEFLEC2022/GC/deviationNorthUncertainty
DDEFLEC2022/GL/deviationEastVelocity
DDEFLEC2022/GL/deviationNorthVelocity
DEM2022/NA/elevation
DEM2022/AS/elevation
DEM2022/GC/elevation
RBA/NA/gravityAnomaly
RBA/NA/gravityAnomalyUncertainty
TC/NA/terrainCorrection
RBA/AS/gravityAnomaly
RBA/AS/gravityAnomalyUncertainty
TC/AS/terrainCorrection
RBA/GC/gravityAnomaly
RBA/GC/gravityAnomalyUncertainty
TC/GC/terrainCorrection
Reading .b files
NGS has a long history of using a binary gridded data format with ‘.b’ extension. More information about this format for non-Fortran users can be found here: NOAA Technical Report NOS NGS 63. Fortran code to read these grids is as follows:
real*8 xlatmin,xlonmin,dphi,dlam
real*4 undu(5401,10801)
read (*) xlatmin,xlonmin,dphi,dlam,ilat,ilon
Do I=1,ilat
xlat=xlatmin+(i-1)*dphi
read(*) (undu(i,j),j=1,ilon)
do j=1,ilon
xlon=xlonmin+(j-1)*dlam
ENDDO
Interpolation
The prescribed method of interpolation for all the NAPGD2022 products is to use biquadratic interpolation. More information can be found here: NOAA Technical Report NOS NGS 84.
Time-dependency - Combining the static and dynamic fields
The static grids in NAPGD2022 (e.g., SGEOID2022) represent the geopotential field at the 2020.0 (January 1st, 2020) epoch. The dynamic components of GEOID2022 and DEFLEC2022 were computed using rates derived from grids timestamped in Julian years, which are exactly 365.25 days. Thus, evaluating NAPGD2022 at a particular epoch requires the user to compute the exact number of Julian years elapsed since the 2020.0 epoch. This may be computed using modified Julian dates to obtain the exact number of days since the epoch and dividing the result by 365.25 days per year.
To compute the time-dependent geoid undulation, both the SGEOID and DGEOID grids must be interpolated to a user-specified latitude and longitude and then combined according to the integer modified Julian day (mjdn) of the computation time:
GEOID = SGEOID + DGEOID*(mjdn-58849)/365.25
The modified Julian day is determined from the year, month, and day of the date of interest according to the following algorithm:
a = (14 - month)//12
y = year + 4800 - a
m = month + 12*a - 3
jdn = day + (153*m+2)//5 + 365*y + y//4 - y//100 + y//400 - 32045
mjdn = jdn - 2400001
The // symbol denotes integer division, e.g. 4//3 = 1 and 2//3 = 0. MJD 58849 corresponds to midnight, January 1st, 2020 or the 2020.0 epoch of the static geoid
Uncertainty propagation with time-dependent components of NAPGD2022 is performed by adding the static and time-dependent uncertainty contributions in quadrature.
delta_t = (mjdn-58849)/365.25
sigmaN = sqrt(geoidHeightUncertainty**2 +
(delta_t*geoidVelocityUncertainty)**2)
Test cases
A spreadsheet of test evaluations of the component models at selected points has been included so users can verify their correct implementation of these models or answer common questions about the datum. These test points include grid boundaries, major cities, and points of geographic and scientific interest.
