This notebook illustrates reprojecting (Lon,Lat,Elevation) points in one 3D CRS to another.
Import libraries and configure logging¶
# This environment variable must be set before importing geopandas for logging
import os
os.environ['PROJ_DEBUG'] = '2'
# Ensure this is 'ON' to get shift grids over the internet
print(os.environ['PROJ_NETWORK'])
import loggingON
import fiona
import geopandas as gpd# It's good to keep track of versions of geospatial libraries and dependencies
gpd.show_versions()Output
SYSTEM INFO
-----------
python : 3.12.11 | packaged by conda-forge | (main, Jun 4 2025, 14:45:31) [GCC 13.3.0]
executable : /home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/bin/python
machine : Linux-6.11.0-1015-azure-x86_64-with-glibc2.39
GEOS, GDAL, PROJ INFO
---------------------
GEOS : 3.13.1
GEOS lib : None
GDAL : 3.10.3
GDAL data dir: /home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/gdal/
PROJ : 9.6.2
PROJ data dir: /home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj
PYTHON DEPENDENCIES
-------------------
geopandas : 1.1.0
numpy : 2.3.0
pandas : 2.3.0
pyproj : 3.7.1
shapely : 2.1.1
pyogrio : 0.11.0
geoalchemy2: None
geopy : None
matplotlib : 3.10.3
mapclassify: 2.9.0
fiona : 1.10.1
psycopg : None
psycopg2 : None
pyarrow : None
Load 3D points¶
# ICESat-2 data saved from sliderule:
#gf = icesat2.atl06p({}, resources=['ATL03_20181019224323_03250112_005_01.h5'])
#gf[:100].to_file('ATL03_20181019224323_03250112_005_01.geojson', driver='GeoJSON')
gf = gpd.read_file('ATL03_20181019224323_03250112_005_01.geojson')gf.head(2)Loading...
gf.crs<Geographic 3D CRS: EPSG:7912>
Name: ITRF2014
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: International Terrestrial Reference Frame 2014
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwichpoints = gf.reset_index()
points.loc[:, 'time'] = points.time.dt.strftime('%Y-%m-%d')
points.explore(zoom_start=2, column='h_mean')Loading...
# Get bounding box of all of our points
w,s,e,n = gf.union_all().bounds #W, S, E, N
print(w,s,e,n)-100.94678307364592 -79.00104942722542 -100.63246345881362 -78.98240211911416
Check reprojection options¶
The projinfo command is very helpful to see which algorithms (or ‘pipelines’) could be used to go from one CRS to another. Below we see Candidate operations found: 75 indicating there are a total of 75 options, which are ordered top to bottom in decreasing preference!
If logging is enabled (PROJ_DEBUG=2) you will see many lines like pj_open_lib(us_noaa_FL.tif) which correspond to PROJ checking for availability of shift grids required for any of the possible transforms. These may either be files in local directories, or retrieved over the network from https://
!PROJ_DEBUG=0 projinfo -s EPSG:7912 -t EPSG:9518 -o PROJ --hide-ballpark --spatial-test intersects | grep CandidateCandidate operations found: 75
!projinfo -s EPSG:7912 -t EPSG:9518 -o PROJ --grid-check none --bbox {w},{s},{e},{n} --hide-ballpark --spatial-test intersects | head -n 20pj_open_lib(proj.ini): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/proj.ini) - succeeded
pj_open_lib(proj.db): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/proj.db) - succeeded
pj_open_lib(us_nga_egm08_25.tif): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/us_nga_egm08_25.tif) - failed
pj_open_lib(egm08_25.gtx): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/egm08_25.gtx) - failed
pj_open_lib(us_nga_egm08_25.tif): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/us_nga_egm08_25.tif) - failed
pj_open_lib(egm08_25.gtx): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/egm08_25.gtx) - failed
Using https://cdn.proj.org/us_nga_egm08_25.tif
pj_open_lib(Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz) - failed
pj_open_lib(Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz) - failed
Cannot open https://cdn.proj.org/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz: HTTP error 404: <?xml version="1.0" encoding="UTF-8"?>
<Error><Code>NoSuchKey</Code><Message>The specified key does not exist.</Message><Key>Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz</Key><RequestId>A6STBKVAPVDGG05J</RequestId><HostId>eOdGT7FBgkcdfHfSGC+zGYPMWrR2Ncoav1/loPEp6MRLV1uF7wMt/ZyvGWRVUlT4Rlkt8SdMvUo=</HostId></Error>
Candidate operations found: 13
-------------------------------------
Operation No. 1:
unknown id, Inverse of Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D) + Inverse of WGS 84 (G2139) to ITRF2014 (1) + Inverse of Conversion from WGS 84 (G2139) (geog3D) to WGS 84 (G2139) (geocentric) + WGS 84 (G2139) to EGM2008 height (from WGS 84 to EGM2008 height (1)) + Inverse of WGS 84 to WGS 84 (G2139), 3.01 m, World
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +inv +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
-------------------------------------
Operation No. 2:
unknown id, Inverse of Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D) + ITRF2014 to ITRF2020 (1) + Inverse of WGS 84 (G2296) to ITRF2020 (1) + Inverse of Conversion from WGS 84 (G2296) (geog3D) to WGS 84 (G2296) (geocentric) + WGS 84 (G2296) to EGM2008 height (from WGS 84 to EGM2008 height (1)) + Inverse of WGS 84 to WGS 84 (G2296), 3.011 m, World
PROJ string:
Reproject data¶
By default geopandas will use the first operation reported by projinfo. In this case:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +inv +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1☝️ +proj=vgridshift +grids=us_nga_egm08_25.tif will apply interpolated vertical offsets corresponding to the repojection of ellipsoid height to geoid height. For the above transform we do not have horizontal position changes, only vertical.
points3D = gpd.points_from_xy(gf.geometry.x, gf.geometry.y, gf.h_mean)
gf3D = gpd.GeoDataFrame(geometry=points3D, crs='EPSG:7912')gf3D[:1].get_coordinates(include_z=True)Loading...
# Reprojection happens here
gfGeoid = gf3D.to_crs(epsg=9518)
gfGeoid[:1].get_coordinates(include_z=True)Loading...
Validate results¶
Once you’re satisfied that the results reflect the expected magnitude of difference, it can be a good idea to add tests to your code to ensure there are no issues in the future. A common gotcha is if you do not have a vertical shift grid locally and there are network connectivity issues the reproject may not actually result in reprojected values!
# X&Y coordinates should be the same, but Z should be different
gpd.pd.testing.assert_frame_equal(gfGeoid.get_coordinates(), gf3D.get_coordinates())
max_dz = (gfGeoid.get_coordinates(include_z=True).z - gf3D.get_coordinates(include_z=True).z).max().astype('int16')
assert max_dz == 29Avoiding bogus reprojection¶
# Here we'll use a target CRS valid only for the United States, but our data is in Antarctica!
!projinfo EPSG:2927+5703 -o WKT2:2019 --single-linepj_open_lib(proj.ini): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/proj.ini) - succeeded
pj_open_lib(proj.db): call fopen(/home/runner/work/3D_CRS_Transformation_Resources/3D_CRS_Transformation_Resources/.pixi/envs/default/share/proj/proj.db) - succeeded
WKT2:2019 string:
COMPOUNDCRS["NAD83(HARN) / Washington South (ftUS) + NAVD88 height",PROJCRS["NAD83(HARN) / Washington South (ftUS)",BASEGEOGCRS["NAD83(HARN)",DATUM["NAD83 (High Accuracy Reference Network)",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4152]],CONVERSION["SPCS83 Washington South zone (US survey foot)",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",45.3333333333333,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",-120.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",47.3333333333333,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",45.8333333333333,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",1640416.667,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["US survey foot",0.304800609601219],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["US survey foot",0.304800609601219]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["US survey foot",0.304800609601219]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["United States (USA) - Washington - counties of Adams; Asotin; Benton; Clark; Columbia; Cowlitz; Franklin; Garfield; Grant south of approximately 47°30'N; Grays Harbor; Kittitas; Klickitat; Lewis; Mason; Pacific; Pierce; Skamania; Thurston; Wahkiakum; Walla Walla; Whitman; Yakima."],BBOX[45.54,-124.4,47.61,-116.91]],ID["EPSG",2927]],VERTCRS["NAVD88 height",VDATUM["North American Vertical Datum 1988"],CS[vertical,1],AXIS["gravity-related height (H)",up,LENGTHUNIT["metre",1]],USAGE[SCOPE["Geodesy, engineering survey, topographic mapping."],AREA["Mexico - onshore. United States (USA) - CONUS and Alaska - onshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming."],BBOX[14.51,172.42,71.4,-66.91]],ID["EPSG",5703]]]
# That is a complicated transform!
!PROJ_DEBUG=0 projinfo -s EPSG:7912 -t EPSG:2927+5703 -q -o PROJ+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=cart +ellps=GRS80
+step +proj=helmert +x=1.0053 +y=-1.90921 +z=-0.54157 +rx=0.02678138
+ry=-0.00042027 +rz=0.01093206 +s=0.00036891 +dx=0.00079 +dy=-0.0006
+dz=-0.00144 +drx=6.667e-05 +dry=-0.00075744 +drz=-5.133e-05
+ds=-7.201e-05 +t_epoch=2010 +convention=coordinate_frame
+step +inv +proj=cart +ellps=GRS80
+step +inv +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1
+step +inv +proj=gridshift +no_z_transform
+grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif
+step +inv +proj=gridshift +no_z_transform
+grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif
+step +inv +proj=gridshift +no_z_transform
+grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif
+step +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333
+lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80
+step +proj=unitconvert +xy_in=m +xy_out=us-ft
logging.basicConfig(level=logging.DEBUG)
with fiona.Env(CPL_DEBUG=True):
bogus = gf3D.to_crs(epsg="2927+5703")DEBUG:fiona._env:GDAL_DATA found in environment.
DEBUG:fiona._env:PROJ_DATA found in environment.
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g2018u0.tif): call fopen(/home/runner/.local/share/proj/us_noaa_g2018u0.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g2018u0.gtx): call fopen(/home/runner/.local/share/proj/g2018u0.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid03_conus.tif): call fopen(/home/runner/.local/share/proj/us_noaa_geoid03_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid03_conus.gtx): call fopen(/home/runner/.local/share/proj/geoid03_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid09_conus.tif): call fopen(/home/runner/.local/share/proj/us_noaa_geoid09_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid09_conus.gtx): call fopen(/home/runner/.local/share/proj/geoid09_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_vertconw.tif): call fopen(/home/runner/.local/share/proj/us_noaa_vertconw.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(vertconw.gtx): call fopen(/home/runner/.local/share/proj/vertconw.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g1999u01.tif): call fopen(/home/runner/.local/share/proj/us_noaa_g1999u01.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g1999u01.gtx): call fopen(/home/runner/.local/share/proj/g1999u01.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g2018u0.tif): call fopen(/home/runner/.local/share/proj/us_noaa_g2018u0.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g2018u0.gtx): call fopen(/home/runner/.local/share/proj/g2018u0.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g1999u01.tif): call fopen(/home/runner/.local/share/proj/us_noaa_g1999u01.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g1999u01.gtx): call fopen(/home/runner/.local/share/proj/g1999u01.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_1993_nad83_1997_prvi.tif): call fopen(/home/runner/.local/share/proj/us_noaa_nadcon5_nad83_1993_nad83_1997_prvi.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_1997_nad83_2002_prvi.tif): call fopen(/home/runner/.local/share/proj/us_noaa_nadcon5_nad83_1997_nad83_2002_prvi.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid03_conus.tif): call fopen(/home/runner/.local/share/proj/us_noaa_geoid03_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid03_conus.gtx): call fopen(/home/runner/.local/share/proj/geoid03_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid09_conus.tif): call fopen(/home/runner/.local/share/proj/us_noaa_geoid09_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid09_conus.gtx): call fopen(/home/runner/.local/share/proj/geoid09_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: Using coordinate operation Inverse of Transformation from NAVD88 height to ITRF2014 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of Ballpark geographic offset from NAD83(HARN) to ITRF2014 + SPCS83 Washington South zone (US survey foot)
bogus.head()Loading...