TP03 - Conversão de Coordenadas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PTR 5003 - Fundamentos de Informações Espaciais\n",
    "\n",
    "_Nome:_ *Fernando Gomes*\n",
    "\n",
    "_NUSP:_ *11863953*\n",
    "\n",
    "> Trabalho Prático - TP 03\n",
    "## Conversão de coordenadas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "São dadas as coordenadas cartesianas de um ponto P1 pertencente a uma rede GPS, no\n",
    "sistema WGS-84:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1_x = 4010243.76\n",
    "p1_y = -4261161.918\n",
    "p1_z = -2532046.653"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bem como as coordenadas geodésicas de um vértice localizado na cidade de Franca –\n",
    "SP, no sistema SAD-69:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "fi = \"-20 34 53.2594\"\n",
    "la = \"-47 22 49.9192\"\n",
    "h_franca = 1013.8947"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## São fornecidos:\n",
    "> Os parâmetros do elipsóide WGS-84:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "wgs84_a = 6378137\n",
    "wgs84_e_2 = 0.00669437999"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Os parâmetros do elipsóide de Referência 1967 (SAD-69):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "sad69_a = 6378160\n",
    "sad69_e_2 = 0.00669454185"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> *Os parâmetros de translação do sistema WGS-84 para o SAD-69:*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "tx = 66.35\n",
    "ty = -03.88\n",
    "tz = 38.22"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pede-se:\n",
    "\n",
    "### as coordenadas cartesianas de P1 referidas ao centro do elipsóide de referência 1967 (SAD-69)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Primeiramente temos que obter os ângulos $\\lambda$ e $\\varphi$  então calcular `X`, `Y` e `Z` referidos ao centro da elipsoide de referência de 1967 (SAD-69)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lambda de P1 é -46.737569700392285\n",
      "Phi de P1 com h = 0 é -23.539722690745275\n",
      "Na primeira aproximação temos h = 1050.4077828126028\n",
      "Phi de P1 com h = 1050.4077828126028 é -23.53969942223693\n",
      "Na segunda aproximação temos h = 1049.2849502796307\n",
      "Phi de P1 com h = 1049.2849502796307 é -23.53969944710565\n",
      "Na terceira aproximação temos h = 1049.2861503306776\n"
     ]
    }
   ],
   "source": [
    "from math import degrees, atan2, sqrt, atan, sin, cos\n",
    "\n",
    "p1_lambda = atan2(p1_y, p1_x)\n",
    "print(f'Lambda de P1 é {degrees(p1_lambda)}')\n",
    "\n",
    "# obtendo p1_phi para h = 0\n",
    "tan_p1_phi = (p1_z / sqrt(p1_x**2 + p1_y**2))*((1 - wgs84_e_2)**-1)\n",
    "p1_phi = atan(tan_p1_phi)\n",
    "print(f'Phi de P1 com h = 0 é {degrees(p1_phi)}')\n",
    "\n",
    "# com o p1_phi obtido calcula-se o N e depois o h\n",
    "N = wgs84_a / sqrt(1 - wgs84_e_2 * sin(p1_phi) * sin(p1_phi))\n",
    "h = sqrt(p1_x**2 + p1_y**2) / cos(p1_phi) - N\n",
    "print(f'Na primeira aproximação temos h = {h}')\n",
    "\n",
    "# com o novo valor de h vamos calcular novamente o p1_phi\n",
    "tan_p1_phi = (p1_z / sqrt(p1_x**2 + p1_y**2))*((1 - wgs84_e_2*(N/(N+h)))**-1)\n",
    "p1_phi = atan(tan_p1_phi)\n",
    "print(f'Phi de P1 com h = {h} é {degrees(p1_phi)}')\n",
    "\n",
    "# com o p1_phi obtido calcula-se o N e depois o h\n",
    "N = wgs84_a / sqrt(1 - wgs84_e_2 * sin(p1_phi) * sin(p1_phi))\n",
    "h = sqrt(p1_x**2 + p1_y**2) / cos(p1_phi) - N\n",
    "print(f'Na segunda aproximação temos h = {h}')\n",
    "\n",
    "# com o novo valor de h vamos calcular novamente o p1_phi\n",
    "tan_p1_phi = (p1_z / sqrt(p1_x**2 + p1_y**2))*((1 - wgs84_e_2*(N/(N+h)))**-1)\n",
    "p1_phi = atan(tan_p1_phi)\n",
    "print(f'Phi de P1 com h = {h} é {degrees(p1_phi)}')\n",
    "\n",
    "# com o p1_phi obtido calcula-se o N e depois o h\n",
    "N = wgs84_a / sqrt(1 - wgs84_e_2 * sin(p1_phi) * sin(p1_phi))\n",
    "h = sqrt(p1_x**2 + p1_y**2) / cos(p1_phi) - N\n",
    "print(f'Na terceira aproximação temos h = {h}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos chegar ao mesmo resultado usando a biblioteca `pyproj`vamos ver como proceder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-46.737569700392285, -23.539699447082114, 1049.286149124615)\n"
     ]
    }
   ],
   "source": [
    "import pyproj\n",
    "from math import sqrt\n",
    "\n",
    "# Determinando o raio menor\n",
    "wgs84_b = wgs84_a*sqrt(1 - wgs84_e_2)\n",
    "\n",
    "ecef = pyproj.Proj(f\"+proj=geocent +a={wgs84_a} +b={wgs84_b} +units=m\")\n",
    "lla = pyproj.Proj(f'+proj=latlong +a={wgs84_a} +b={wgs84_b}')\n",
    "xyz_wgs84_pyproj = pyproj.transform(ecef, lla, p1_x, p1_y, p1_z)\n",
    "print(xyz_wgs84_pyproj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agora que temos as coordenadas geodésicas em latitude e longitude podemos obter as coordenadas cartesianas em referência a elipsoide de referência de 1967 (SAD69)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temos phi = -0.4108452602818824, lambda = -0.8157244756521851, N = 6381545.03746993 e h = 1049.2861503306776\n"
     ]
    }
   ],
   "source": [
    "print(f'Temos phi = {p1_phi}, lambda = {p1_lambda}, N = {N} e h = {h}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RESPOSTA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As cordenadas cartesianas de P1 referidas ao centro da elipsoide de 1967 (SAD69) são 4010258.2706521163, -4261177.3365735095 e -2532055.4024095153\n"
     ]
    }
   ],
   "source": [
    "from math import cos, sin, radians, degrees\n",
    "\n",
    "p1_phi = degrees(p1_phi)\n",
    "p1_lambda = degrees(p1_lambda)\n",
    "\n",
    "N = sad69_a / sqrt(1 - sad69_e_2 * sin(radians(p1_phi)) * sin(radians(p1_phi)))\n",
    "\n",
    "X = (N + h) * cos(radians(p1_phi)) * cos(radians(p1_lambda))\n",
    "Y = (N + h) * cos(radians(p1_phi)) * sin(radians(p1_lambda))\n",
    "Z = ((1 - sad69_e_2) * N + h) * sin(radians(p1_phi))\n",
    "\n",
    "print(f\"As cordenadas cartesianas de P1 referidas ao centro da elipsoide de 1967 (SAD69) são {X}, {Y} e {Z}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos chegar ao mesmo resultado usando a biblioteca `pyproj`vamos ver como proceder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4010258.2706520706, -4261177.336573461, -2532055.402406647)"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pyproj\n",
    "from math import sqrt\n",
    "\n",
    "# Determinando o raio menor\n",
    "wgs84_b = wgs84_a*sqrt(1 - wgs84_e_2)\n",
    "sad69_b = sad69_a*sqrt(1 - sad69_e_2)\n",
    "\n",
    "ecef_wgs84 = pyproj.Proj(f\"+proj=geocent +a={wgs84_a} +b={wgs84_b} +units=m\")\n",
    "lla_wgs84 = pyproj.Proj(f'+proj=latlong +a={wgs84_a} +b={wgs84_b}')\n",
    "ecef_sad69 = pyproj.Proj(f\"+proj=geocent +a={sad69_a} +b={sad69_b} +units=m\")\n",
    "lla_sad69 = pyproj.Proj(f'+proj=latlong +a={sad69_a} +b={sad69_b}')\n",
    "\n",
    "p1_lat, p1_long, p1_altura = pyproj.transform(ecef_wgs84, lla_wgs84, p1_x, p1_y, p1_z)\n",
    "\n",
    "pyproj.transform(lla_sad69, ecef_sad69, p1_lat, p1_long, p1_altura)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### as coordenadas geodésicas de P1 no sistema SAD-69."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como chagamos ao valores de X, Y e Z referidos ao SAD69 podemos simplesmente obter os valores de $\\varphi$, $\\lambda$ e `h` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Temos os valores de X, Y e Z sendo respectivamente: 4010258.2706521163, -4261177.3365735095 e -2532055.4024095153 no sistema cartesiado referenciado a elipsoide SAD69\n"
     ]
    }
   ],
   "source": [
    "print(f'Temos os valores de X, Y e Z sendo respectivamente: {X}, {Y} e {Z} no sistema cartesiado referenciado a elipsoide SAD69')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RESPOSTA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lambda de P1 em SAD69 é -46.737569700392285\n",
      "Phi de P1 em SAD69 com h = 0 é -23.539722691253537\n",
      "Na primeira aproximação temos h = 1050.4078112421557 para SAD69\n",
      "Phi de P1  em SAD69 com h = 1050.4078112421557 é -23.53969942226237\n",
      "Na segunda aproximação temos h = 1049.2849515024573 para SAD69\n",
      "Phi de P1 em SAD69 com h = 1049.2849515024573 é -23.539699447132207\n",
      "Na terceira aproximação temos h = 1049.2861516121775 para SAD69\n"
     ]
    }
   ],
   "source": [
    "from math import degrees, atan2, sqrt, atan, sin, cos\n",
    "\n",
    "p1_lambda = atan2(Y, X)\n",
    "print(f'Lambda de P1 em SAD69 é {degrees(p1_lambda)}')\n",
    "\n",
    "# obtendo p1_phi para h = 0\n",
    "tan_p1_phi = (Z/ sqrt(X**2 + Y**2))*((1 - sad69_e_2)**-1)\n",
    "p1_phi = atan(tan_p1_phi)\n",
    "print(f'Phi de P1 em SAD69 com h = 0 é {degrees(p1_phi)}')\n",
    "\n",
    "# com o p1_phi obtido calcula-se o N e depois o h\n",
    "N = sad69_a / sqrt(1 - sad69_e_2 * sin(p1_phi) * sin(p1_phi))\n",
    "h = sqrt(X**2 + Y**2) / cos(p1_phi) - N\n",
    "print(f'Na primeira aproximação temos h = {h} para SAD69')\n",
    "\n",
    "# com o novo valor de h vamos calcular novamente o p1_phi\n",
    "tan_p1_phi = (Z / sqrt(X**2 + Y**2))*((1 - sad69_e_2*(N/(N+h)))**-1)\n",
    "p1_phi = atan(tan_p1_phi)\n",
    "print(f'Phi de P1  em SAD69 com h = {h} é {degrees(p1_phi)}')\n",
    "\n",
    "# com o p1_phi obtido calcula-se o N e depois o h\n",
    "N = sad69_a / sqrt(1 - sad69_e_2 * sin(p1_phi) * sin(p1_phi))\n",
    "h = sqrt(X**2 + Y**2) / cos(p1_phi) - N\n",
    "print(f'Na segunda aproximação temos h = {h} para SAD69')\n",
    "\n",
    "# com o novo valor de h vamos calcular novamente o p1_phi\n",
    "tan_p1_phi = (Z / sqrt(X**2 + Y**2))*((1 - sad69_e_2*(N/(N+h)))**-1)\n",
    "p1_phi = atan(tan_p1_phi)\n",
    "print(f'Phi de P1 em SAD69 com h = {h} é {degrees(p1_phi)}')\n",
    "\n",
    "# com o p1_phi obtido calcula-se o N e depois o h\n",
    "N = sad69_a / sqrt(1 - sad69_e_2 * sin(p1_phi) * sin(p1_phi))\n",
    "h = sqrt(X**2 + Y**2) / cos(p1_phi) - N\n",
    "print(f'Na terceira aproximação temos h = {h} para SAD69')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-46.737569700392285, -23.539699447105658, 1049.2861503325403)"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pyproj\n",
    "import math\n",
    "\n",
    "# Determinando o raio menor\n",
    "sad69_b = sad69_a*math.sqrt(1 - sad69_e_2)\n",
    "\n",
    "ecef = pyproj.Proj(f\"+proj=geocent +a={sad69_a} +b={sad69_b} +units=m\")\n",
    "lla = pyproj.Proj(f'+proj=latlong +a={sad69_a} +b={sad69_b}')\n",
    "pyproj.transform(ecef, lla, X, Y, Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### as coordenadas cartesianas de Franca referidas ao WGS-84."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nesse caso primeiro temos que converter os ângulos que estão em GMS (Graus, minutos e segundos) para decimal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Para o ângulo GMS -20 34 53.2594 agora temos o ângulo decimal -20.581460944444444\n",
      "Para o ângulo GMS -47 22 49.9192 agora temos o ângulo decimal -47.38053311111111\n"
     ]
    }
   ],
   "source": [
    "import numpy as np \n",
    "\n",
    "# função de converção de GMS (Grau, Minuto, Segundo) para ângulo decimal\n",
    "def gms_para_decimal(angulo):\n",
    "    divisor = np.asarray((1, 60, 3600))\n",
    "    gms = np.fromstring(angulo, dtype=float, sep=' ')\n",
    "    operador = gms[0] / abs(gms[0])\n",
    "    return np.sum(np.abs(gms / divisor)) * operador\n",
    "\n",
    "fi_dec = gms_para_decimal(fi)\n",
    "la_dec = gms_para_decimal(la)\n",
    "\n",
    "print(f\"Para o ângulo GMS {fi} agora temos o ângulo decimal {fi_dec}\")\n",
    "print(f\"Para o ângulo GMS {la} agora temos o ângulo decimal {la_dec}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dessa forma agora podemos utilizar as seguintes fórmulas para calcular `X`, `Y` e `Z`\n",
    "$$X = (N + h)\\cos\\varphi\\cos\\lambda$$\n",
    "$$Y = (N + h)\\cos\\varphi\\sin\\lambda$$\n",
    "$$Z = [(1 - e²)N + h]\\sin\\varphi$$\n",
    "Onde para finalidade desse código $\\varphi=$ `fi_dec` e $\\lambda=$ `la_dec`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RESPOSTA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As cordenadas cartesianas de P1 referidas ao WGS84 são 4045478.59534743, -4396425.999580176 e -2228438.368709564\n"
     ]
    }
   ],
   "source": [
    "from math import cos, sin, radians\n",
    "\n",
    "N = sad69_a / math.sqrt(1 - sad69_e_2 * sin(radians(fi_dec)) * sin(radians(fi_dec)))\n",
    "\n",
    "X = (N + h_franca) * cos(radians(fi_dec)) * cos(radians(la_dec))\n",
    "Y = (N + h_franca) * cos(radians(fi_dec)) * sin(radians(la_dec))\n",
    "Z = ((1 - sad69_e_2) * N + h_franca) * sin(radians(fi_dec))\n",
    "\n",
    "print(f\"As cordenadas cartesianas de P1 referidas ao WGS84 são {X}, {Y} e {Z}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos chegar ao mesmo resultado usando a biblioteca `pyproj`vamos ver como proceder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyproj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4045478.59534743, -4396425.999580176, -2228438.368709564\n"
     ]
    }
   ],
   "source": [
    "print(f'{X}, {Y}, {Z}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As coordenadas cartesianas do ponto em Franca referida a SAD69 são 4045478.5953474296, -4396425.999580176 e -2228438.368709567\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "# Determinando o raio menor\n",
    "sad69_b = sad69_a*math.sqrt(1 - sad69_e_2)\n",
    "wgs84_b = wgs84_a*math.sqrt(1 - wgs84_e_2)\n",
    "\n",
    "ecef_wgs84 = pyproj.Proj(f\"+proj=geocent +a={wgs84_a} +b={wgs84_b} +units=m\")\n",
    "lla_wgs84 = pyproj.Proj(f'+proj=latlong +a={wgs84_a} +b={wgs84_b}')\n",
    "ecef_sad69 = pyproj.Proj(f\"+proj=geocent +a={sad69_a} +b={sad69_b} +units=m\")\n",
    "lla_sad69 = pyproj.Proj(f'+proj=latlong +a={sad69_a} +b={sad69_b}')\n",
    "\n",
    "f1_x, f1_y, f1_z = pyproj.transform(lla_sad69, ecef_wgs84, la_dec, fi_dec, h_franca)\n",
    " \n",
    "print(f'As coordenadas cartesianas do ponto em Franca referida a SAD69 são {f1_x}, {f1_y} e {f1_z}')"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}