00_ocentric_ographic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GDAL Map Projection Mars Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Requirements:\n",
    "> conda install -c conda-forge gdal=3 jupyter pyproj requests\n",
    "\n",
    "**Data for lesson**\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gathering information from GDAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GDAL's version is: 3.1.0\n"
     ]
    }
   ],
   "source": [
    "# First, from GDAL's Python API, let's check GDAL's version.\n",
    "# This will also double check our GDAL Python environment is working.\n",
    "\n",
    "# Import the \"gdal\" submodule from within the \"osgeo\" module\n",
    "from osgeo import gdal\n",
    "\n",
    "# We can check which version we're running by printing the \"__version__\" variable\n",
    "print(\"GDAL's version is: \" + gdal.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (2) Lets make a proj4 string for a spherical Mars (in degrees).\n",
    "\n",
    "Mars2000_Sphere = \"+proj=longlat +a=3396190 +b=3396190 +no_defs\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can find out what Lon/Lat maps into X/Y meters in the image's map projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs \n",
      "\n"
     ]
    }
   ],
   "source": [
    "import requests\n",
    "response = requests.get('https://spatialreference.org/ref/iau2000/49915/proj4/')\n",
    "IAU49915 = response.content.decode(\"utf-8\")\n",
    "print(IAU49915 + \"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 -2644524.7064723633 0.0\n"
     ]
    }
   ],
   "source": [
    "# Now in Python, you can also use TransformPoint\n",
    "from osgeo import osr\n",
    "\n",
    "s_SRS = osr.SpatialReference()  # makes an empty spatial ref object\n",
    "s_SRS.ImportFromProj4(Mars2000_Sphere)   # populates the spatial ref object with our WKT SRS\n",
    "\n",
    "t_SRS = osr.SpatialReference()  # makes an empty spatial ref object\n",
    "t_SRS.ImportFromProj4(IAU49915) # populates the spatial ref using Proj4 SRS\n",
    "\n",
    "# create the transform - this can be used repeatedly\n",
    "transform = osr.CoordinateTransformation(s_SRS, t_SRS)\n",
    "\n",
    "# transform the point. You can also create an ogr geometry and use the more generic `point.Transform()`\n",
    "                    # order:   Lon , Lat ,  Z\n",
    "X,Y,Z = transform.TransformPoint(0.0, -45.0, 0.0)\n",
    "\n",
    "print(X,Y,Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 -45.0 0.0\n"
     ]
    }
   ],
   "source": [
    "# Now go back to degrees\n",
    "\n",
    "# create the transform - this can be used repeatedly\n",
    "transform = osr.CoordinateTransformation(t_SRS, s_SRS)\n",
    "\n",
    "# transform the point. You can also create an ogr geometry and use the more generic `point.Transform()`\n",
    "                    # order:   Lon , Lat ,  Z\n",
    "X,Y,Z = transform.TransformPoint(X, Y, Z)\n",
    "\n",
    "print(X,Y,Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From degrees to geocent and then geocent_sphere\n",
      "2408546.880462214 0.0 2380276.874501003\n",
      "2401468.979197941 0.0 2401468.979197941\n",
      "\n",
      "From degrees to cart and then cart_sphere\n",
      "2408546.880462214 0.0 2380276.874501003\n",
      "2401468.979197941 0.0 2401468.979197941\n",
      "\n",
      " note geocent and cart ellipse/sphere are the same\n"
     ]
    }
   ],
   "source": [
    "from osgeo import osr\n",
    "\n",
    "s_SRS = osr.SpatialReference()  # makes an empty spatial ref object\n",
    "t_SRS = osr.SpatialReference()  # makes an empty spatial ref object\n",
    "\n",
    "Mars2000        = \"+proj=longlat +a=3396190 +b=3376200 +no_defs\"\n",
    "# Use of Mars2000_sphere doesn't impact anything\n",
    "Mars2000_sphere = \"+proj=longlat +a=3396190 +b=3396190 +no_defs\"\n",
    "\n",
    "geocent         = \"+proj=geocent +a=3396190 +b=3376200\"\n",
    "cart            = \"+proj=cart    +a=3396190 +b=3376200 +lon_0=0 +units=m +no_defs +wktext\"\n",
    "\n",
    "geocent_sphere  = \"+proj=geocent +a=3396190 +b=3396190\"\n",
    "cart_sphere     = \"+proj=cart    +a=3396190 +b=3396190 +lon_0=0 +units=m +no_defs +wktext\"\n",
    "\n",
    "s_SRS.ImportFromProj4(Mars2000)   # populates the spatial ref object with our WKT SRS\n",
    "t_SRS.ImportFromProj4(geocent)    # populates the spatial ref using Proj4 SRS\n",
    "# create the transform - this can be used repeatedly\n",
    "transform = osr.CoordinateTransformation(s_SRS, t_SRS)\n",
    "# transform the point. You can also create an ogr geometry and use the more generic `point.Transform()`\n",
    "                      # order:   Lon , Lat ,  Z\n",
    "X,Y,Z = transform.TransformPoint(0.0, 45.0, 0.0)\n",
    "print(\"From degrees to geocent and then geocent_sphere\")\n",
    "print(X,Y,Z)\n",
    "\n",
    "s_SRS.ImportFromProj4(Mars2000)       # populates the spatial ref object with our WKT SRS\n",
    "t_SRS.ImportFromProj4(geocent_sphere) # populates the spatial ref using Proj4 SRS\n",
    "transform = osr.CoordinateTransformation(s_SRS, t_SRS)\n",
    "X,Y,Z = transform.TransformPoint(0.0, 45.0, 0.0)\n",
    "print(X,Y,Z)\n",
    "print()\n",
    "\n",
    "s_SRS.ImportFromProj4(Mars2000)   # populates the spatial ref object with our WKT SRS\n",
    "t_SRS.ImportFromProj4(cart)       # populates the spatial ref using Proj4 SRS\n",
    "transform = osr.CoordinateTransformation(s_SRS, t_SRS)\n",
    "X,Y,Z = transform.TransformPoint(0.0, 45.0, 0.0)\n",
    "print(\"From degrees to cart and then cart_sphere\")\n",
    "print(X,Y,Z)\n",
    "\n",
    "s_SRS.ImportFromProj4(Mars2000)    # populates the spatial ref object with our WKT SRS\n",
    "t_SRS.ImportFromProj4(cart_sphere) # populates the spatial ref using Proj4 SRS\n",
    "transform = osr.CoordinateTransformation(s_SRS, t_SRS)\n",
    "X,Y,Z = transform.TransformPoint(0.0, 45.0, 0.0)\n",
    "print(X,Y,Z)\n",
    "print(\"\\n note geocent and cart ellipse/sphere are the same\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Rel. 7.0.0, March 1st, 2020\n",
      "<proj>: \n",
      "missing argument for -e\n",
      "program abnormally terminated\n"
     ]
    }
   ],
   "source": [
    "!proj -version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2.6.1.post1'"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pyproj\n",
    "pyproj.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 44.6617680466192\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Apps\\Anaconda3\\envs\\gdal3.1\\lib\\site-packages\\pyproj\\transformer.py:430: UserWarning: radian input with pipelines is not supported in pyproj 2. support for raidans will be added in pyproj 3.\n",
      "  self._transformer._transform(\n"
     ]
    }
   ],
   "source": [
    "# from Jay (and SnowMan2) on https://gis.stackexchange.com/questions/366676/proj-pyproj-convert-between-ocentric-and-ographic-latitudes\n",
    "import numpy as np\n",
    "\n",
    "trans = pyproj.transformer.Transformer.from_pipeline('+proj=pipeline +a=3396190 +b=3376200 +step +proj=geoc')\n",
    "lon, lat = trans.transform(0, np.radians(45.0), errcheck=True, radians=True)\n",
    "print(np.degrees(lon), np.degrees(lat))\n",
    "\n",
    "#from stack exchange\n",
    "#44.6617680466192"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the use of \"radians=True\" is very important here despite the warning.\n",
    "\n",
    "### Now for the inverse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 45.0\n"
     ]
    }
   ],
   "source": [
    "itrans = pyproj.transformer.Transformer.from_pipeline('+proj=pipeline +a=3396190 +b=3376200 +step +proj=geoc +inv')\n",
    "nlon, nlat = itrans.transform(lon, lat, errcheck=True, radians=True)\n",
    "print(np.degrees(nlon), np.degrees(nlat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To double check the answer, here is an alternative method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def oc2og(dlat, dMajorRadius = 3396190.0, dMinorRadius = 3376200.0):\n",
    "    try:\n",
    "        dlat = math.radians(dlat)\n",
    "        dlat = math.atan(((dMajorRadius / dMinorRadius)**2) * (math.tan(dlat)))\n",
    "        dlat = math.degrees(dlat)\n",
    "    except:\n",
    "        print (\"Error in oc2og conversion\")\n",
    "    return dlat\n",
    "\n",
    "def og2oc(dlat, dMajorRadius = 3396190.0, dMinorRadius = 3376200.0):\n",
    "    try:\n",
    "        dlat = math.radians(dlat)\n",
    "        dlat = math.atan((math.tan(dlat) / ((dMajorRadius / dMinorRadius)**2)))\n",
    "        dlat = math.degrees(dlat)\n",
    "    except:\n",
    "        print (\"Error in oc2og conversion\")\n",
    "    return dlat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44.6617680466192\n",
      "45.00000000000001\n"
     ]
    }
   ],
   "source": [
    "dlat = og2oc(45.0)\n",
    "print(dlat)\n",
    "\n",
    "nlat = oc2og(dlat)\n",
    "print(nlat)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}