07-Geodesy_Geographic_Geocentric_Coordinate_Transformation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "cb952e6e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-2381006.7075355262, 37.47944444076724, -3284342.9515491314, -3607.801921026781)\n"
     ]
    }
   ],
   "source": [
    "#Conversion from Geocentric(x,y,z)m to Geodetic(Grids)degrees\n",
    "#Clarke 1880 (6378249.145, 6356514.966)\n",
    "#first eccentricity (e2)\n",
    "#second eccentricity (e12)\n",
    "\n",
    "from math import sin, cos, atan, sqrt, radians, degrees\n",
    "clarke1880=(6378249.145, 6356514.966)\n",
    "\n",
    "def geocentric_geodetic (ell, x, y, z):\n",
    "    a, b= ell\n",
    "    e2 = (a**2 - b**2)/a**2\n",
    "    e12 = (a**2 - b**2)/b**2\n",
    "    p = sqrt (x**2 + y**2)\n",
    "    lat1 = degrees(atan(z*a/ p*b))\n",
    "    sin_lat = sin(lat1**3)\n",
    "    cos_lat = cos(lat1**3)\n",
    "    lon = degrees(atan(y/x))\n",
    "    lat = degrees(atan((z+e12*b**sin_lat))/ p-e2*a*cos_lat)\n",
    "    v = a / sqrt(1- e2 * sin(lat**2))\n",
    "    h= (z/sin(lat) - v * (1-e2))\n",
    "    h2 = ((p/cos(lat))-v)\n",
    "    return (lat, lon, h, h2 )\n",
    "x = 5058876.955\n",
    "y = 3878930.060\n",
    "z= 249870.507\n",
    "#fing h\n",
    "\n",
    "print (geocentric_geodetic(clarke1880, x,y,z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "7d51cf1f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4979478.88814634, 3818051.024692411, 246805.9040887472)\n"
     ]
    }
   ],
   "source": [
    "#conversion from geodetic to geocentric\n",
    "from math import sin ,  cos, atan, sqrt, radians, degrees\n",
    "wgs84 = (6278137, 298.257223563)\n",
    "\n",
    "def geodetic_geocentric (ell, latitude, longitude, height ):\n",
    "    a, rf = ell\n",
    "    f = 1/ rf\n",
    "    e2 = f * (1-f)\n",
    "    lat = radians(latitude)\n",
    "    lon = radians(longitude)\n",
    "    sin_lat = sin(lat**2)\n",
    "    h = height\n",
    "    v = a/ sqrt (1-e2 * sin_lat)\n",
    "    x = ((v+h)*cos(lat) * cos(lon))\n",
    "    y = ((v+h)*cos(lat) * sin(lon))\n",
    "    z= (v*(1-e2)+ h)*sin(lat)\n",
    "    \n",
    "    return (x,y,z)\n",
    "\n",
    "latitude = 2.26000000\n",
    "longitude = 37.47944444\n",
    "height = 1500\n",
    "\n",
    "print(geodetic_geocentric(wgs84, latitude, longitude, height))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "c93be492",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(89.99867895563354, 2.1309652352296395, 6392835.764589443)\n"
     ]
    }
   ],
   "source": [
    "#GEOCENTRIC TRANSLATION\n",
    "#Consider a North Sea Point...\n",
    "\n",
    "from math import sin, cos, atan, sqrt, radians, degrees\n",
    "#Conversion from wgs84 geodetic to geocentric\n",
    " \n",
    "wgs84 = (6378137, 6356752.3145)    \n",
    "def geodetic_geocentric (ell, latitude, longitude, height):\n",
    "    #conversion from wgs84 to geocentric\n",
    "    a,b = ell\n",
    "    e2 = (a**2 - b**2)/a**2\n",
    "    lat = radians(latitude)\n",
    "    lon = radians(longitude)\n",
    "    sin_lat = sin(lat**2)\n",
    "    h = height\n",
    "    v = a/ sqrt(1-e2 *sin_lat)\n",
    "    x = (v + h)* cos(lat) * cos(lon)\n",
    "    y = (v+h) * cos(lat) *sin(lon)\n",
    "    z= (v*(1-e2)+ h)*sin(lat)\n",
    "    \n",
    "    #conversion from wgs84 to ED50\n",
    "    \n",
    "    xt = x + 84.87\n",
    "    yt = y + 96.49\n",
    "    zt = z + 116.95\n",
    "    \n",
    "    #conversion of ED50 from geocentric to geodetic\n",
    "    a2 = 6378388\n",
    "    rf2 = 297.0\n",
    "    f2 = 1 / rf2\n",
    "    ed2 = f2*(2-f2)\n",
    "    b2 = sqrt (a2**2*(1-ed2))\n",
    "    ed12 = (a2**2-b2**2)/ b2**2\n",
    "    p2 = sqrt (xt**2 + yt**2)\n",
    "    lat11 = degrees(atan(zt*a2/p2*b2))\n",
    "    sin_lat11 = sin(lat**3)\n",
    "    cos_lat11 = cos(lat**3)\n",
    "    lon = degrees(atan (yt/xt))\n",
    "    lat = degrees(atan((zt*ed12*b2*sin_lat11)/(p2-ed2*a2*cos_lat11)))\n",
    "    \n",
    "    v2 = a2/ sqrt (1-ed2 * sin(lat**2))\n",
    "    h = h= (zt/sin(lat) - v2 * (1-ed2))\n",
    "    return (lat, lon, v2)\n",
    "latitude = 53.80939444\n",
    "longitude = 2.12955000\n",
    "height = 73\n",
    "print (geodetic_geocentric(wgs84, latitude, longitude, height ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7772059",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}