proto_dist.ipynb

proto_dist.ipynb

{
  "cells": [
    {
      "metadata": {
        "trusted": false
      },
      "id": "9022b250",
      "cell_type": "code",
      "source": "import numpy as np\n\n\nlat_arr = np.array([52.11, -52.31, 24, -24])\nlon_arr = np.array([312.44, 132.54, 313, 300])",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "id": "a2a30b27",
      "cell_type": "code",
      "source": "def current_great_circle_distance(lat_arr, lon_arr):\n    from pyproj import Geod\n    dist = np.ma.zeros(lon_arr.size, dtype=np.float64)\n    g = Geod(ellps='WGS84')\n    _, _, dist[1:] = g.inv(lon_arr[:-1], lat_arr[:-1], lon_arr[1:], lat_arr[1:])\n    return dist",
      "execution_count": 2,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "id": "6bb0dd02",
      "cell_type": "code",
      "source": "def new_great_circle_distance(lat_arr, lon_arr):\n    from geographiclib.geodesic import Geodesic\n    dist = np.ma.zeros(lon_arr.size, dtype=np.float64)\n    positions = list(zip(lat_arr, lon_arr))\n    positions_pairs = list(zip(positions[:-1], positions[1:]))\n    dist[1:] = [Geodesic.WGS84.Inverse(*pair[0], *pair[1])[\"s12\"] for pair in positions_pairs]\n    return dist",
      "execution_count": 3,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "id": "6ce588fe",
      "cell_type": "code",
      "source": "current = current_great_circle_distance(lat_arr, lon_arr)\nnew = new_great_circle_distance(lat_arr, lon_arr)\n\nnp.testing.assert_array_almost_equal(new, current, decimal=9)  # They are not equal to 10 decimal places.",
      "execution_count": 4,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "id": "a14886e9",
      "cell_type": "code",
      "source": "new",
      "execution_count": 5,
      "outputs": [
        {
          "data": {
            "text/plain": "masked_array(data=[       0.        , 19981187.20119626,\n                   16861142.39264485,  5493010.77922261],\n             mask=False,\n       fill_value=1e+20)"
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {
        "trusted": false
      },
      "id": "bfa8048d",
      "cell_type": "code",
      "source": "current",
      "execution_count": 6,
      "outputs": [
        {
          "data": {
            "text/plain": "masked_array(data=[       0.        , 19981187.20119626,\n                   16861142.39264485,  5493010.77922261],\n             mask=False,\n       fill_value=1e+20)"
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3 (ipykernel)",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.10.5",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "proto_dist.ipynb",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}