pyproj and geopandas

geopandas-3Dtransforms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8036a086-fa08-42a2-8fd1-af9a5951e18c",
   "metadata": {},
   "source": [
    "# Proj pipelines with geopandas\n",
    "\n",
    "Goal: 3D geodetic transforms with geopandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "57449c56-cf84-459f-8294-7e406a9f8823",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "SYSTEM INFO\n",
      "-----------\n",
      "python     : 3.11.3 | packaged by conda-forge | (main, Apr  6 2023, 08:57:19) [GCC 11.3.0]\n",
      "executable : /home/jovyan/.local/envs/pandas2/bin/python\n",
      "machine    : Linux-5.4.238-148.346.amzn2.x86_64-x86_64-with-glibc2.35\n",
      "\n",
      "GEOS, GDAL, PROJ INFO\n",
      "---------------------\n",
      "GEOS       : 3.11.2\n",
      "GEOS lib   : None\n",
      "GDAL       : 3.6.4\n",
      "GDAL data dir: /home/jovyan/.local/envs/pandas2/share/gdal\n",
      "PROJ       : 9.2.0\n",
      "PROJ data dir: /home/jovyan/.local/envs/pandas2/share/proj\n",
      "\n",
      "PYTHON DEPENDENCIES\n",
      "-------------------\n",
      "geopandas  : 0.13.0\n",
      "numpy      : 1.24.3\n",
      "pandas     : 2.0.1\n",
      "pyproj     : 3.5.0\n",
      "shapely    : 2.0.1\n",
      "fiona      : 1.9.3\n",
      "geoalchemy2: None\n",
      "geopy      : None\n",
      "matplotlib : 3.7.1\n",
      "mapclassify: 2.5.0\n",
      "pygeos     : None\n",
      "pyogrio    : 0.6.0\n",
      "psycopg2   : None\n",
      "pyarrow    : 12.0.0\n",
      "rtree      : 1.0.1\n"
     ]
    }
   ],
   "source": [
    "from pyproj.transformer import TransformerGroup, Transformer\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "gpd.show_versions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d08902a5-655b-4d5b-9678-2b4e67b9da1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<TransformerGroup: best_available=True>\n",
       "- transformers: 5\n",
       "- unavailable_operations: 0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get pipeline you're interested in\n",
    "\n",
    "tg = TransformerGroup(\"EPSG:7912\", \"EPSG:2927+5703\", \n",
    "                              always_xy= True, \n",
    "                              allow_ballpark= False\n",
    "                             )\n",
    "tg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "470bf307-e45a-46bc-93c3-3ad7759889c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.105"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Select a specific available transform\n",
    "transformer = tg.transformers[0] # 1st = best accuracy\n",
    "transformer.accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "300baaab-c19e-4771-a4c3-2efe4d921186",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=pipeline \n",
      "+step +proj=unitconvert +xy_in=deg +xy_out=rad \n",
      "+step +proj=cart +ellps=GRS80 \n",
      "+step +proj=helmert +x=1.0053 +y=-1.9092 +z=-0.5416 +rx=0.0267814 +ry=-0.0004203 +rz=0.0109321 +s=0.00037 +dx=0.0008 +dy=-0.0006 +dz=-0.0014 +drx=6.67e-05 +dry=-0.0007574 +drz=-5.13e-05 +ds=-7e-05 +t_epoch=2010 +convention=coordinate_frame \n",
      "+step +inv +proj=cart +ellps=GRS80 \n",
      "+step +inv +proj=vgridshift +grids=us_noaa_g2018u0.tif +multiplier=1 \n",
      "+step +inv +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif \n",
      "+step +inv +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif \n",
      "+step +inv +proj=gridshift +no_z_transform +grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif \n",
      "+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 \n",
      "+step +proj=unitconvert +xy_in=m +xy_out=us-ft\n"
     ]
    }
   ],
   "source": [
    "pipeline = transformer.to_proj4()\n",
    "print(pipeline.replace('+step', '\\n+step'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a5d7e6b3-cb63-438d-9092-554274e6a561",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.       ,  0.       ,  0.0392002],\n",
       "       [ 0.       ,  0.       , -0.0672624]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example of forcing specific shift grid \n",
    "\n",
    "pipeline = '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=1.0053 +y=-1.9092 +z=-0.5416 +rx=0.0267814 +ry=-0.0004203 +rz=0.0109321 +s=0.00037 +dx=0.0008 +dy=-0.0006 +dz=-0.0014 +drx=6.67e-05 +dry=-0.0007574 +drz=-5.13e-05 +ds=-7e-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'\n",
    "custom_pipeline = pipeline.replace('us_noaa_g2018u0.tif', 'us_noaa_g2012bu0.tif')\n",
    "\n",
    "transformer = Transformer.from_pipeline(pipeline)\n",
    "custom_transformer = Transformer.from_pipeline(custom_pipeline)\n",
    "\n",
    "gf = gpd.GeoDataFrame(geometry=gpd.points_from_xy([-120.4, -120.52502], [48.6, 48.48523], [1400, 2234.139058]),\n",
    "                  crs='EPSG:7912'\n",
    "                 )\n",
    "\n",
    "\n",
    "points2012 = custom_transformer.transform(gf.geometry.x.values, gf.geometry.y.values, gf.geometry.z.values)\n",
    "points2018 = transformer.transform(gf.geometry.x.values, gf.geometry.y.values, gf.geometry.z.values)\n",
    "\n",
    "np.stack(points2012).T - np.stack(points2018).T  #x,y,z points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "667b452e-3d10-46b4-8500-c9d4cc32d260",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:.local-pandas2]",
   "language": "python",
   "name": "conda-env-.local-pandas2-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}