Exploring interconversion between pyproj and rasterio CRS objects

pyproj_rasterio_crs_interconversion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploring CRS interconversion between `rasterio.crs.CRS` and `pyproj.crs.CRS`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rasterio Version: 1.0.23\n",
      "PyProj Version: 2.5.0\n"
     ]
    }
   ],
   "source": [
    "import rasterio\n",
    "import pyproj\n",
    "print(f'Rasterio Version: {rasterio.__version__}')\n",
    "print(f'PyProj Version: {pyproj.__version__}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epsg:4326\n"
     ]
    }
   ],
   "source": [
    "pp_crs = pyproj.CRS.from_epsg(4326)\n",
    "print(pp_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EPSG:4326\n"
     ]
    }
   ],
   "source": [
    "rio_crs = rasterio.crs.CRS.from_epsg(4326)\n",
    "print(rio_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_crs == rio_crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's great that these can be directly compared to one another. The problems arise when we try to interconvert between the classes, which is going to be necessary for standardized CRS management in `solaris`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "int() argument must be a string, a bytes-like object or a number, not 'CRS'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-da50520c04b5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mp_to_r_crs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcrs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCRS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_epsg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpp_crs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/miniconda3/envs/solaris/lib/python3.6/site-packages/rasterio/crs.py\u001b[0m in \u001b[0;36mfrom_epsg\u001b[0;34m(cls, code)\u001b[0m\n\u001b[1;32m    302\u001b[0m         \"\"\"\n\u001b[1;32m    303\u001b[0m         \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 304\u001b[0;31m         \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_crs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_CRS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_epsg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    305\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    306\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mrasterio/_crs.pyx\u001b[0m in \u001b[0;36mrasterio._crs._CRS.from_epsg\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: int() argument must be a string, a bytes-like object or a number, not 'CRS'"
     ]
    }
   ],
   "source": [
    "p_to_r_crs = rasterio.crs.CRS.from_epsg(pp_crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That doesn't work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'CRS' object is not iterable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-63da7816a873>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mp_to_r_crs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcrs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCRS\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpp_crs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/miniconda3/envs/solaris/lib/python3.6/site-packages/rasterio/crs.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, initialdata, **kwargs)\u001b[0m\n\u001b[1;32m     65\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0minitialdata\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     66\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 67\u001b[0;31m             \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minitialdata\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     68\u001b[0m             \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     69\u001b[0m             \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mall_proj_keys\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: 'CRS' object is not iterable"
     ]
    }
   ],
   "source": [
    "p_to_r_crs = rasterio.crs.CRS(pp_crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Neither does that. As far as I can tell, the best way to interconvert is to use WKT strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS[\"ellipsoidal\",2],AXIS[\"geodetic latitude (Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433]],AXIS[\"geodetic longitude (Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]],USAGE[SCOPE[\"unknown\"],AREA[\"World\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]\n"
     ]
    }
   ],
   "source": [
    "p_to_r_crs_wkt = rasterio.crs.CRS.from_wkt(pp_crs.to_wkt())\n",
    "print(p_to_r_crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reads a lot differently from the earlier `rio_crs` object. Will it be identified as the same?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rio_crs == p_to_r_crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nope. It's also worth noting that `rasterio` can't identify the EPSG code of this object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rio_crs EPSG code: 4326\n",
      "p_to_r_crs EPSG code: None\n"
     ]
    }
   ],
   "source": [
    "print(f'rio_crs EPSG code: {rio_crs.to_epsg()}')\n",
    "print(f'p_to_r_crs EPSG code: {p_to_r_crs.to_epsg()}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EPSG:4326\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_to_r_epsg = rasterio.crs.CRS.from_epsg(pp_crs.to_epsg())\n",
    "print(p_to_r_epsg)\n",
    "p_to_r_epsg == rio_crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That works, but it means we _must_ use CRSs/SRSs that have EPSG codes.\n",
    "\n",
    "What if we do the conversion the other way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]\n"
     ]
    }
   ],
   "source": [
    "r_to_p_crs = pyproj.CRS.from_wkt(rio_crs.to_wkt())\n",
    "print(r_to_p_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r_to_p_crs == pp_crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This seems to work better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/envs/solaris/lib/python3.6/site-packages/pyproj/crs/crs.py:55: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n",
      "  return _prepare_from_string(\" \".join(pjargs))\n"
     ]
    }
   ],
   "source": [
    "crs_dict = {'init': 'epsg:4326'}\n",
    "pp_crs_d = pyproj.CRS(crs_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs_str = 'epsg:4326'\n",
    "pp_crs_epsgstr = pyproj.CRS(crs_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_crs_code = pyproj.CRS(4326)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_crs_d == pp_crs_epsgstr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is annoying - things that use the (admittedly deprecated) '+init' format aren't identified as identical to the preferred syntax for the same object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Geographic 2D CRS: +init=epsg:4326 +type=crs>\n",
       "Name: WGS 84\n",
       "Axis Info [ellipsoidal]:\n",
       "- lon[east]: Longitude (degree)\n",
       "- lat[north]: Latitude (degree)\n",
       "Area of Use:\n",
       "- name: World\n",
       "- bounds: (-180.0, -90.0, 180.0, 90.0)\n",
       "Datum: World Geodetic System 1984\n",
       "- Ellipsoid: WGS 84\n",
       "- Prime Meridian: Greenwich"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_crs_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Geographic 2D CRS: EPSG:4326>\n",
       "Name: WGS 84\n",
       "Axis Info [ellipsoidal]:\n",
       "- Lat[north]: Geodetic latitude (degree)\n",
       "- Lon[east]: Geodetic longitude (degree)\n",
       "Area of Use:\n",
       "- name: World\n",
       "- bounds: (-180.0, -90.0, 180.0, 90.0)\n",
       "Datum: World Geodetic System 1984\n",
       "- Ellipsoid: WGS 84\n",
       "- Prime Meridian: Greenwich"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_crs_epsgstr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the best part as I see it is that you can pass just about anything into the `__init__` method of a pyproj.CRS object and get the same thing out (with the obvious exception shown above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_crs_str = pyproj.CRS('epsg:4326')\n",
    "pp_crs_rio = pyproj.CRS(rio_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp_crs == pp_crs_str and pp_crs == pp_crs_rio and pp_crs == pp_crs_epsgstr and pp_crs == pp_crs_code"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "solaris",
   "language": "python",
   "name": "solaris"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}