Geodesic calculations

geodesics.py

# -*- coding: utf-8 -*-

# Copyright © 2022 by Gregory Werbin
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED “AS IS” AND ISC DISCLAIMS ALL WARRANTIES WITH REGARD
# TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
# FITNESS. IN NO EVENT SHALL ISC BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT,
# OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
# USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
# TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
# OF THIS SOFTWARE.

r"""Comparison of geodesic calculation and conversion tools.

Requirements:

  • CPython               : 3.8+
  • Proj                  : 9.1.0
  • PyProj                : 3.4.0
  • Geographiclib         : 1.49
  • Geographiclib (Python): 2.0
  • Nvector (Python)      : 0.7.7
"""

## Setup ################################################################# {{{

from __future__ import annotations
import logging
import shlex
import subprocess

import numpy as np
import nvector
import pyproj
from geographiclib.geodesic import Geodesic

logger = logging.getLogger(__name__)


def forward_geod(
    proj4_args: tuple[str, ...], lat1: float, lon1: float, az: float, dist: float,
) -> tuple[float, float]:
    cmd = ("geod", "-f", "%0.8f", "-F", "%0.8f", *proj4_args)
    input_data = f"{lat1:0.8f} {lon1:0.8f} {az:0.8f} {dist*1000:0.8f}"

    logger.info('%s <<< %s', shlex.join(cmd), shlex.quote(input_data))

    output = subprocess.check_output(cmd, input=input_data, text=True)
    return tuple(map(float, output.split(maxsplit=2)[:2]))


def forward_geodsolve(
    a: float, f: float, lat1: float, lon1: float, az: float, dist: float,
) -> tuple[float, float]:
    cmd = ("GeodSolve", "-e", format(a, "0.8f"), format(f, "0.8f"))
    input_data = f"{lat1:0.8f} {lon1:0.8f} {az:0.8f} {dist*1000:0.8f}"

    logger.info('%s <<< %s', shlex.join(cmd), shlex.quote(input_data))

    output = subprocess.check_output(cmd, input=input_data, text=True)
    return tuple(map(float, output.split(maxsplit=2)[:2]))


def forward_geographiclib(
    geodesic: Geodesic, lat1: float, lon1: float, az: float, dist: float,
) -> tuple[float, float]:
    out = geodesic.Direct(in_lat_deg, in_lon_deg, az, dist * 1000)
    return (out["lat2"], out["lon2"])


def forward_pyproj(
    geod: pyproj.Geod, lat1: float, lon1: float, az: float, dist: float,
) -> tuple[float, float]:
    out = geod.fwd(in_lat_deg, in_lon_deg, az, dist * 1000)
    return (out[0], out[1])


def latlon_to_nvect(lat, lon):
    return nvector.lat_lon2n_E(*map(np.radians, (lat, lon)))

def nvect_to_latlon(nvect):
    return tuple(map(np.degrees, nvector.n_E2lat_lon(nvect)))

def forward_nvect(
    a: float, f: float, lat: float, lon: float, az: float, dist: float
) -> tuple[float, float]:
    if f != 0.0:
        raise NotImplementedError("N-vector forward solution not yet implemented for ellipsoid.")

    nvect = latlon_to_nvect(lat, lon)
    az_rad = np.radians(az)
    dist_angle_rad = dist * 1000 / a

    out = nvector.n_EA_E_distance_and_azimuth2n_EB_E(nvect, dist_angle_rad, az_rad)
    return tuple(float(x[0]) for x in nvect_to_latlon(out))


def format_latlon(lat: float, lon: float, dec: int = 8) -> str:
    return f"{lat:0.{dec}f} °N, {lon:0.{dec}f} °E"

# }}}

## Settings ############################################################## {{{

# Copied from man page GeodSolve(1)
wgs84_equatorial_radius_m = 6378137
wgs84_flattening = 1 / 298.257223563

# Copied from H3 source code: https://github.com/uber/h3/blob/v3.7.2/src/h3lib/include/constants.h#L58-L59
wgs84_authalic_radius_m = 6371.007180918475 * 1000

a_f_ellpse = (wgs84_equatorial_radius_m , wgs84_flattening)
a_f_sphere = (wgs84_authalic_radius_m   ,              0.0)

geodesic_ellpse = Geodesic(*a_f_ellpse)
geodesic_sphere = Geodesic(*a_f_sphere)

# WTF: proj(1) doesn't work when I use +datum: "Ellipse setup failure"
proj4_ellpse = ("+proj=latlon", "+ellps=WGS84")
proj4_sphere = ("+proj=latlon", "+ellps=WGS84", "+R_A")
#proj4_datume = ("+proj=latlon", "+datum=WGS84")
#proj4_sphere = ("+proj=latlon", "+datum=WGS84", "+R_A")

geod_ellpse = pyproj.Geod(a=wgs84_equatorial_radius_m, f=wgs84_flattening)
geod_sphere = pyproj.Geod(a=wgs84_authalic_radius_m, f=0)

crs_ellpse = pyproj.CRS(proj="latlon", datum="WGS84")
crs_sphere = pyproj.CRS(proj="latlon", datum="WGS84", R_A=True)
#crs_ellpse = pyproj.CRS(proj="latlon", a=wgs84_equatorial_radius_m, f=wgs84_flattening)
#crs_sphere = pyproj.CRS(proj="latlon", a=wgs84_authalic_radius_m, f=0)

geod_ellpse_derived = crs_ellpse.get_geod()
geod_sphere_derived = crs_sphere.get_geod()

# WTF: pyproj.Transformer().from_crs has no effect when I use +ellps instead of +datum: outputs are identical to inputs!
trans_ellps_sphere = pyproj.Transformer.from_proj(
    "+proj=latlon +datum=WGS84",
    "+proj=latlon +datum=WGS84 +R_A",
)
# trans_ellps_sphere = pyproj.Transformer.from_proj(
#     "+proj=latlon +datum=WGS84",
#     "+proj=latlon +datum=WGS84 +R_A",
# )
# trans_ellps_sphere = pyproj.Transformer.from_crs(crs_ellps, crs_sphere)
# trans_ellps_sphere = pyproj.Transformer.from_crs(
#     pyproj.CRS.from_proj4("+proj=latlon +datum=WGS84"),
#     pyproj.CRS.from_proj4("+proj=latlon +datum=WGS84 +R_A"),
# )
# trans_ellps_sphere = pyproj.Transformer.from_crs(
#     pyproj.CRS.from_proj4("+proj=latlon +ellps=WGS84"),
#     pyproj.CRS.from_proj4("+proj=latlon +ellps=WGS84 +R_A"),
# )

# }}}

## Input data ############################################################ {{{

# Boston City Hall Plaza
in_lat_deg  = 42.36035
in_lon_deg = -75.05918

bearing_deg = 210
distance_km = 500

# }}}

## Run ################################################################### {{{

logging.basicConfig(
    level=logging.INFO,
    format="░▒▓ [%(levelname)s | %(filename)s : %(funcName)15s]\t%(message)s"
)

print("Input:", in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(flush=True)

## Ellipsoidal {{{

print("Ellipsoidal:")

out_lat_deg, out_lon_deg = forward_geodsolve(*a_f_ellpse, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  GeodSolve(1):  {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_geographiclib(geodesic_ellpse, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  Geographiclib: {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_geod(proj4_ellpse, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  geod(1):       {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_pyproj(geod_ellpse, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  Pyproj:        {format_latlon(out_lat_deg, out_lon_deg)}")

print(flush=True)

## }}}

## Spherical {{{

print("Spherical:")

out_lat_deg, out_lon_deg = forward_geodsolve(*a_f_sphere, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  GeodSolve(1):  {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_geographiclib(geodesic_sphere, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  Geographiclib: {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_geod(proj4_sphere, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  geod(1):       {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_pyproj(geod_sphere, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  Pyproj:        {format_latlon(out_lat_deg, out_lon_deg)}")

out_lat_deg, out_lon_deg = forward_nvect(*a_f_sphere, in_lat_deg, in_lon_deg, bearing_deg, distance_km)
print(f"  Nvector:       {format_latlon(out_lat_deg, out_lon_deg)}")

print(flush=True)

## }}}

## Transformation: ellipsoidal → spherical {{{

print("ellipsoidal → spherical:")
print("  pyproj.Transformer:")
out_lat_deg, out_lon_deg = trans_ellps_sphere.transform(in_lat_deg, in_lon_deg)
out_lat_deg = round(out_lat_deg, 5)
out_lon_deg = round(out_lon_deg, 5)
print(f"    (input)  {format_latlon(in_lat_deg, in_lon_deg)}")
print(f"    (output) {format_latlon(out_lat_deg, out_lon_deg)}")

print("  cs2cs(1):")
cmd = (
    "cs2cs", "-f", "%0.7f",
    "+proj=latlon", "+datum=WGS84", "+to",
    "+proj=latlon", "+datum=WGS84", "+R_A",
)
input_data = format_latlon(in_lat_deg, in_lon_deg)
logger.info('%s <<< %s', shlex.join(cmd), shlex.quote(input_data))
out = subprocess.check_output(cmd, input=input_data, text=True)
out_lat_deg, out_lon_deg = (round(float(val), 7) for val in out.split(maxsplit=2)[:2])
print(f"    (input)  {format_latlon(in_lat_deg, in_lon_deg)}")
print(f"    (output) {format_latlon(out_lat_deg, out_lon_deg)}")

## }}}

## }}}

# vim: set shiftwidth=4 expandtab foldmethod=marker fileencoding=utf8:

output.txt

Input: 42.36035 -75.05918 210 500

Ellipsoidal:
░▒▓ [INFO | geodesics.py : forward_geodsolve]   GeodSolve -e 6378137.00000000 0.00335281 <<< '42.36035000 -75.05918000 210.00000000 500000.00000000'
  GeodSolve(1):  38.42388682 °N, -77.92032251 °E
  Geographiclib: 38.42388681 °N, -77.92032251 °E
░▒▓ [INFO | geodesics.py :    forward_geod]     geod -f %0.8f -F %0.8f +proj=latlon +ellps=WGS84 <<< '42.36035000 -75.05918000 210.00000000 500000.00000000'
  geod(1):       38.42388681 °N, -77.92032251 °E
  Pyproj:        30.85176677 °N, -78.71870807 °E

Spherical:
░▒▓ [INFO | geodesics.py : forward_geodsolve]   GeodSolve -e 6371007.18091847 0.00000000 <<< '42.36035000 -75.05918000 210.00000000 500000.00000000'
  GeodSolve(1):  38.42920316 °N, -77.92744560 °E
  Geographiclib: 38.42920316 °N, -77.92744560 °E
░▒▓ [INFO | geodesics.py :    forward_geod]     geod -f %0.8f -F %0.8f +proj=latlon +ellps=WGS84 +R_A <<< '42.36035000 -75.05918000 210.00000000 500000.00000000'
  geod(1):       38.42920316 °N, -77.92744560 °E
  Pyproj:        30.78742318 °N, -78.73196191 °E
  Nvector:       38.42920316 °N, -77.92744560 °E

ellipsoidal → spherical:
  pyproj.Transformer:
    (input)  42.36035000 °N, -75.05918000 °E
    (output) 42.36035000 °N, -74.96303000 °E
  cs2cs(1):
░▒▓ [INFO | geodesics.py :        <module>]     cs2cs -f %0.7f +proj=latlon +datum=WGS84 +to +proj=latlon +datum=WGS84 +R_A <<< '42.36035000 °N, -75.05918000 °E'
    (input)  42.36035000 °N, -75.05918000 °E
    (output) 42.36035000 °N, 0.00000000 °E