3D Coordinate Reference System Transformation Resources

# NOTE: REQUIRES PROJ>=9.2 and GEOPANDAS>=0.13

#!pip install geopandas==0.13

import os
os.environ['USE_PYGEOS'] = '0'
os.environ['PROJ_DEBUG'] = '2' # For showing logs
#os.environ['PROJ_NETWORK'] = 'ON' # Ensure this is 'ON' to get shift grids over the internet
print(os.environ['PROJ_NETWORK']) 

import fiona
import geopandas as gpd
import shapely

ON

gpd.show_versions()


SYSTEM INFO
-----------
python     : 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:40:32) [GCC 12.3.0]
executable : /srv/conda/envs/notebook/bin/python
machine    : Linux-5.10.167-147.601.amzn2.x86_64-x86_64-with-glibc2.27

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : 3.11.2
GEOS lib   : None
GDAL       : 3.7.0
GDAL data dir: /srv/conda/envs/notebook/share/gdal
PROJ       : 9.2.1
PROJ data dir: /srv/conda/envs/notebook/share/proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.13.0
numpy      : 1.23.5
pandas     : 1.5.1
pyproj     : 3.6.0
shapely    : 2.0.1
fiona      : 1.9.4
geoalchemy2: None
geopy      : None
matplotlib : 3.6.2
mapclassify: 2.5.0
pygeos     : 0.14
pyogrio    : None
psycopg2   : 2.9.6 (dt dec pq3 ext lo64)
pyarrow    : 12.0.1
rtree      : 1.0.1

# from slidrule single granule, default parameters
#gf = icesat2.atl06p({}, resources=['ATL03_20181019224323_03250112_005_01.h5'])
gf = gpd.read_file('sliderule-output.gpkg',
                  rows=100)

gf.geometry.head(2)

0    POINT (-100.63246 -78.98922)
1    POINT (-100.63265 -78.98904)
Name: geometry, dtype: geometry

# Always check wehere you data is!
points = gf.reset_index()
points.loc[:, 'time'] = points.time.dt.strftime('%Y-%m-%d')
points.explore()

# Get bounding box of all of our points
w,s,e,n = gf.unary_union.bounds #W, S, E, N
print(w,s,e,n)

-100.94678306984892 -79.00104942913823 -100.63246339163007 -78.98240211911293

#!projinfo -s EPSG:7912 -t EPSG:9518 -o PROJ --hide-ballpark --spatial-test intersects
!projinfo -s EPSG:7912 -t EPSG:9518 -o PROJ --bbox {w},{s},{e},{n}  --hide-ballpark --spatial-test intersects

pj_open_lib(proj.ini): call fopen(/srv/conda/envs/notebook/share/proj/proj.ini) - succeeded
pj_open_lib(proj.db): call fopen(/srv/conda/envs/notebook/share/proj/proj.db) - succeeded
pj_open_lib(us_noaa_nadcon5_nad27_nad83_1986_conus.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_nadcon5_nad27_nad83_1986_conus.tif) - failed
pj_open_lib(proj.db): call fopen(/srv/conda/envs/notebook/share/proj/proj.db) - succeeded
pj_open_lib(us_noaa_nadcon5_nad83_1986_nad83_harn_conus.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_nadcon5_nad83_1986_nad83_harn_conus.tif) - failed
pj_open_lib(us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif) - failed
pj_open_lib(us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif) - failed
pj_open_lib(us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif) - failed
pj_open_lib(us_noaa_alhpgn.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_alhpgn.tif) - failed
pj_open_lib(alhpgn.gsb): call fopen(/srv/conda/envs/notebook/share/proj/alhpgn.gsb) - failed
pj_open_lib(us_noaa_conus.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_conus.tif) - failed
pj_open_lib(conus): call fopen(/srv/conda/envs/notebook/share/proj/conus) - failed
pj_open_lib(us_noaa_FL.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_FL.tif) - failed
pj_open_lib(FL): call fopen(/srv/conda/envs/notebook/share/proj/FL) - failed
pj_open_lib(us_noaa_lahpgn.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_lahpgn.tif) - failed
pj_open_lib(lahpgn.gsb): call fopen(/srv/conda/envs/notebook/share/proj/lahpgn.gsb) - failed
pj_open_lib(us_noaa_mshpgn.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_mshpgn.tif) - failed
pj_open_lib(mshpgn.gsb): call fopen(/srv/conda/envs/notebook/share/proj/mshpgn.gsb) - failed
pj_open_lib(us_noaa_ethpgn.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_noaa_ethpgn.tif) - failed
pj_open_lib(ethpgn.gsb): call fopen(/srv/conda/envs/notebook/share/proj/ethpgn.gsb) - failed
pj_open_lib(au_icsm_GDA94_GDA2020_conformal_and_distortion.tif): call fopen(/srv/conda/envs/notebook/share/proj/au_icsm_GDA94_GDA2020_conformal_and_distortion.tif) - failed
pj_open_lib(GDA94_GDA2020_conformal_and_distortion.gsb): call fopen(/srv/conda/envs/notebook/share/proj/GDA94_GDA2020_conformal_and_distortion.gsb) - failed
pj_open_lib(eur_nkg_nkgrf17vel.tif): call fopen(/srv/conda/envs/notebook/share/proj/eur_nkg_nkgrf17vel.tif) - failed
pj_open_lib(no_kv_NKGETRF14_EPSG7922_2000.tif): call fopen(/srv/conda/envs/notebook/share/proj/no_kv_NKGETRF14_EPSG7922_2000.tif) - failed
pj_open_lib(Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz): call fopen(/srv/conda/envs/notebook/share/proj/Und_min1x1_egm2008_isw=82_WGS84_TideFree.gz) - failed
pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed
pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed
pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed
pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed
Using https://cdn.proj.org/us_nga_egm08_25.tif
Candidate operations found: 1
-------------------------------------
Operation No. 1:

pj_open_lib(us_nga_egm08_25.tif): call fopen(/srv/conda/envs/notebook/share/proj/us_nga_egm08_25.tif) - failed
pj_open_lib(egm08_25.gtx): call fopen(/srv/conda/envs/notebook/share/proj/egm08_25.gtx) - failed
Using https://cdn.proj.org/us_nga_egm08_25.tif
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) (geog2D) to WGS 84 (G2139) (geocentric) + Inverse of WGS 84 to WGS 84 (G2139) + Inverse of Null geographic offset from WGS 84 (geog3D) to WGS 84 (geog2D) + WGS 84 to EGM2008 height (1), 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

# NOTE: only vshift so only expect vertical coordinate change
points3D = gpd.points_from_xy(gf.geometry.x, gf.geometry.y, gf.h_mean)
gf3D = gpd.GeoDataFrame(geometry=points3D, crs='EPSG:7912')
gf3D.head(2)

gf3D.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: Greenwich

shapely.get_coordinates(gf3D.iloc[0], include_z=True)

array([[-100.63246339, -78.98921541, 1875.21727596]])

gfGeoid = gf3D.to_crs(epsg=9518)
shapely.get_coordinates(gfGeoid.iloc[0], include_z=True)

array([[-100.63246339, -78.98921541, 1904.67380708]])

#gfGeoid.head()

# SHow debug log

# CAREFUL, geopandas will happily do bogus transforms for you

import logging
logging.basicConfig(level=logging.DEBUG)

with fiona.Env(CPL_DEBUG=True):
    test = 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_nadcon5_nad83_harn_nad83_fbn_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g2018u0.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_g2018u0.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g2018u0.gtx): call fopen(/home/jovyan/.local/share/proj/g2018u0.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid03_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_geoid03_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid03_conus.gtx): call fopen(/home/jovyan/.local/share/proj/geoid03_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid09_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_geoid09_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid09_conus.gtx): call fopen(/home/jovyan/.local/share/proj/geoid09_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_vertconw.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_vertconw.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(vertconw.gtx): call fopen(/home/jovyan/.local/share/proj/vertconw.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g1999u01.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_g1999u01.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g1999u01.gtx): call fopen(/home/jovyan/.local/share/proj/g1999u01.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g2018u0.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_g2018u0.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g2018u0.gtx): call fopen(/home/jovyan/.local/share/proj/g2018u0.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_g2018u0.tif
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_1993_nad83_1997_prvi.tif): call fopen(/home/jovyan/.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/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_1997_nad83_2002_prvi.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid03_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_geoid03_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid03_conus.gtx): call fopen(/home/jovyan/.local/share/proj/geoid03_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_geoid03_conus.tif
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_geoid09_conus.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_geoid09_conus.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(geoid09_conus.gtx): call fopen(/home/jovyan/.local/share/proj/geoid09_conus.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_geoid09_conus.tif
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(us_noaa_g1999u01.tif): call fopen(/home/jovyan/.local/share/proj/us_noaa_g1999u01.tif) - failed
DEBUG:pyproj:PROJ_DEBUG: pj_open_lib(g1999u01.gtx): call fopen(/home/jovyan/.local/share/proj/g1999u01.gtx) - failed
DEBUG:pyproj:PROJ_DEBUG: Using https://cdn.proj.org/us_noaa_g1999u01.tif
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 feet)

test.head()

Examples

3D transformation of rasters using GDAL and proj libraries

Global DEMs

Brief overview of CRSs used for global datasets