I use pyproj to calculate the geoid undulation at some coordinates where I measured some data. Today I noticed that the undulation is always calculated to be zero when using version 3.7.2 while it looks as expected for version 3.71. varying between -40 and 60 meters over Europe.
Attached you'll find a minimal working example and pictures.
Please excuse If I provide too much or too little information. This is my first bug report ever...
[geoid_undulation_bug_report_MWE.py](https://github.com/user-attachments/files/27515221/geoid_undulation_bug_report_MWE.py)
import numpy as np
import matplotlib.pyplot as plt
from pyproj import CRS, Transformer
from mpl_toolkits.basemap import Basemap
def calcOrthometricHeights(lat, lon, ellipsoidal_height):
'''
Calculates the orthometric height from a set of GPS coordinates made of lat, lon and ellipsoidal WGS84 height
Returns the orthometric height (height above mean sea level) using EGM2008
'''
wgs84 = CRS("EPSG:4979")
egm08 = CRS("EPSG:3855")
# create a transformer from WGS84 to EGM08
transformer = Transformer.from_crs(wgs84, egm08, always_xy=True)
# Compute orthometric heights
orthometric_heights = []
for lon, lat, h_ell in zip(lon, lat, ellipsoidal_height):
# Transform (lon, lat, ellipsoidal_height) → (lon, lat, orthometric_height)
_, _, h_ortho = transformer.transform(lon, lat, h_ell)
orthometric_heights.append(h_ortho)
orthometric_height = np.array(orthometric_heights)
geoid_undulation = ellipsoidal_height - orthometric_height
return orthometric_height, geoid_undulation
# Calculate the geoid undulation
max_lat = 65
min_lat = 20
max_lon = 60
min_lon = -30
lower_left_corner = (min_lat, min_lon) # lat, lon
upper_right_corner = (max_lat, max_lon)
grid_lat = np.arange(lower_left_corner[0], upper_right_corner[0], 0.1)
grid_lon = np.arange(lower_left_corner[1], upper_right_corner[1], 0.1)
lats, lons = np.meshgrid(grid_lat, grid_lon)
# Calculate the orthometric height at ellipsoidal height 0
ellipsodial_heigth = np.zeros((lats.shape[0], lons.shape[1]))
orthometric_height, undulation = calcOrthometricHeights(lats, lons, ellipsoidal_height=ellipsodial_heigth)
print(undulation)
# Plot the results
lat_min, lat_max = lower_left_corner[0], upper_right_corner[0]
lon_min, lon_max = lower_left_corner[1], upper_right_corner[1]
fig, axs = plt.subplots(1, 1, figsize=(15, 10), constrained_layout=True)
ax = axs
m = Basemap(projection='cyl',
llcrnrlat=lat_min, urcrnrlat=lat_max,
llcrnrlon=lon_min, urcrnrlon=lon_max,
resolution='i', ax=ax)
m.drawcoastlines()
m.drawcountries()
m.drawparallels(np.arange(lat_min, lat_max+1, 5), labels=[1,0,0,0], fontsize=8)
m.drawmeridians(np.arange(lon_min, lon_max+1, 10), labels=[0,0,0,1], fontsize=8)
cs = m.pcolormesh(lons, lats, undulation, shading='auto', cmap='turbo', latlon=True)
ax.set_title('Geoid undulation (EGM2008)', fontsize=12)
cbar = fig.colorbar(cs, ax=ax, orientation='vertical', pad=0.02, aspect=30, shrink=0.8)
cbar.set_label('Geoid undulation (m)', fontsize=10)
cbar.ax.tick_params(labelsize=8)
plt.show()
For version 3.7.1 the undulation looks as expected, varying between -40 and 60 meters in Europe:
For version 3.7.2 the undulation is always 0, which is wrong.
$ conda list pyproj
# Name Version Build Channel
pyproj 3.7.2 py313hbf73894_3 conda-forge
Dear developers,
I use pyproj to calculate the geoid undulation at some coordinates where I measured some data. Today I noticed that the undulation is always calculated to be zero when using version 3.7.2 while it looks as expected for version 3.71. varying between -40 and 60 meters over Europe.
Attached you'll find a minimal working example and pictures.
Please excuse If I provide too much or too little information. This is my first bug report ever...
Thank you!
Code Example
Problem description
For version 3.7.1 the undulation looks as expected, varying between -40 and 60 meters in Europe:
For version 3.7.2 the undulation is always 0, which is wrong.
Expected Output
The undulation should look like for version 3.7.1
Installation method
Miniforge
Conda environment information (if you installed with conda):