Calculate Geoidal Undulation

This (notebook) computes and plots geoid undulations from a gravity model

Python Dependencies

• numpy: Scientific Computing Tools For Python

• matplotlib: Python 2D plotting library

• cartopy: Python package designed for geospatial data processing

• lxml: processing XML and HTML in Python

Program Dependencies

• utilities.py: download and management utilities for syncing files

• geoid_undulation.py: geoidal undulation at a given latitude and longitude

• read_ICGEM_harmonics.py: reads the coefficients for a given gravity model file

• calculate_tidal_offset.py: calculates the C20 offset for a tidal system

• real_potential.py: real potential at a latitude and height for gravity model

• norm_potential.py: normal potential of an ellipsoid at a latitude and height

• norm_gravity.py: normal gravity of an ellipsoid at a latitude and height

• ref_ellipsoid.py: Computes parameters for a reference ellipsoid

• gauss_weights.py: Computes Gaussian weights as a function of degree

import numpy as np
import matplotlib

matplotlib.rcParams['axes.linewidth'] = 2.0
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = ['Helvetica']
matplotlib.rcParams['mathtext.default'] = 'regular'
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
import matplotlib.ticker as ticker
import cartopy.crs as ccrs
import ipywidgets as widgets
import geoid_toolkit.utilities
from geoid_toolkit.read_ICGEM_harmonics import read_ICGEM_harmonics
from geoid_toolkit.geoid_undulation import geoid_undulation

Choose gravity model

# list gfc models from GFZ ICGEM
MODELS = geoid_toolkit.utilities.icgem_list()
modelDropdown = widgets.Dropdown(
    options=sorted(MODELS.keys()),
    value='GGM05C',
    description='Model:',
    disabled=False,
)
display(modelDropdown)

Read gravity model coefficients

# gfc file with spherical harmonic coefficients
MODEL = MODELS[modelDropdown.value]
GRAVITY = geoid_toolkit.utilities.get_data_path(['data', MODEL[-1]])
MD5 = geoid_toolkit.utilities.get_hash(GRAVITY)
# Download coefficients from GFZ ICGEM server
geoid_toolkit.utilities.from_http(
    ['http://icgem.gfz-potsdam.de', *MODEL],
    local=GRAVITY,
    hash=MD5,
    verbose=True,
)
# use maximum degree and order of model
LMAX = None
# use original tide system
TIDE = None
# read gravity model Ylms and change tide if specified
Ylms = read_ICGEM_harmonics(GRAVITY, LMAX=LMAX, TIDE=TIDE)
# extract parameters
R = np.float64(Ylms['radius'])
GM = np.float64(Ylms['earth_gravity_constant'])
LMAX = np.int64(Ylms['max_degree']) if not LMAX else LMAX

Calculate map of geoid height

# PURPOSE: calculate geoid heights at a set of latitudes and longitudes
dlon, dlat = (1.0, 1.0)
lon = np.arange(-180 + dlon / 2.0, 180 + dlon / 2.0, dlon)
lat = np.arange(-90 + dlat / 2.0, 90 + dlat / 2.0, dlat)
nlon = len(lon)
nlat = len(lat)
# reference to WGS84 ellipsoid
REFERENCE = 'WGS84'
# Gaussian Smoothing Radius in km (default is no filtering)
# no gaussian smoothing
GAUSS = 0
# calculate geoid at coordinates
N = np.zeros((nlat, nlon))
for i in range(nlat):
    N[i, :] = geoid_undulation(
        np.ones((nlon)) * lat[i],
        lon,
        REFERENCE,
        Ylms['clm'],
        Ylms['slm'],
        LMAX,
        R,
        GM,
        GAUSS=GAUSS,
    )

Create output plot

# setup Plate Carree projection
fig, ax1 = plt.subplots(
    num=1,
    nrows=1,
    ncols=1,
    figsize=(10.375, 6.625),
    subplot_kw=dict(projection=ccrs.PlateCarree()),
)

# contours
PRANGE = (-80, 80, 20)
levels = np.arange(PRANGE[0], PRANGE[1] + PRANGE[2], PRANGE[2])
norm = colors.Normalize(vmin=PRANGE[0], vmax=PRANGE[1])

# plot image with transparency using normalization
im = ax1.imshow(
    N,
    interpolation='nearest',
    cmap=cm.viridis_r,
    extent=(lon.min(), lon.max(), lat.min(), lat.max()),
    norm=norm,
    alpha=1.0,
    transform=ccrs.PlateCarree(),
    origin='lower',
)

# add generic coastlines
ax1.coastlines()

# draw lat/lon grid lines
GRID = [15, 15]
grid_meridians = np.arange(0, 360 + GRID[0], GRID[0])
grid_parallels = np.arange(-90, 90 + GRID[1], GRID[1])
gl = ax1.gridlines(
    crs=ccrs.PlateCarree(),
    draw_labels=False,
    linewidth=0.1,
    color='0.25',
    linestyle='-',
)
gl.xlocator = ticker.FixedLocator(grid_meridians)
gl.ylocator = ticker.FixedLocator(grid_parallels)

# Add horizontal colorbar and adjust size
# extend = add extension triangles to upper and lower bounds
# options: neither, both, min, max
# pad = distance from main plot axis
# shrink = percent size of colorbar
# aspect = lengthXwidth aspect of colorbar
cbar = plt.colorbar(
    im,
    ax=ax1,
    extend='both',
    extendfrac=0.0375,
    orientation='horizontal',
    pad=0.025,
    shrink=0.90,
    aspect=22,
    drawedges=False,
)
# rasterized colorbar to remove lines
cbar.solids.set_rasterized(True)
# Add label to the colorbar
cbar.ax.set_xlabel('Geoidal Undulation', labelpad=10, fontsize=20)
cbar.ax.set_ylabel('m', fontsize=20, rotation=0)
cbar.ax.yaxis.set_label_coords(1.04, 0.15)
# Set the tick levels for the colorbar
cbar.set_ticks(levels)
cbar.set_ticklabels([f'{ct:d}' for ct in levels])
# ticks lines all the way across
cbar.ax.tick_params(
    which='both', width=1, length=27, labelsize=20, direction='in'
)

# axis = equal
ax1.set_aspect('equal', adjustable='box')
# no ticks on the x and y axes
ax1.get_xaxis().set_ticks([])
ax1.get_yaxis().set_ticks([])
# add main title
ax1.set_title(modelDropdown.value, fontsize=24)
ax1.title.set_y(1.01)

# stronger linewidth on frame
ax1.outline_patch.set_linewidth(2.0)
ax1.outline_patch.set_capstyle('projecting')
# output to file
fig.subplots_adjust(left=0.04, right=0.96, bottom=0.05, top=0.96)
plt.show()