geoutils.py

import numpy as np
from dataclasses import dataclass

from earth_models import EarthModel, WGS84, GRS80


@dataclass
class GeoPoint(Datum):
    lon: LatLonType = field(default_factory=LatLonType)
    lat: LatLonType = field(default_factory=LatLonType)
    height: LatLonType = field(default_factory=LatLonType)


def degrees2radians(x):
    return x * np.pi / 180.0


def radians2degrees(x):
    return x * 180.0 / np.pi


@dataclass
class XYZ:
    x: int
    y: int
    z: int


def geo2xyz(latitude, longitude, altitude, geo_system: str) -> XYZ:
    """
    Convert geodetic coordinate into cartesian XYZ coordinate with specified geodetic system.

    # param latitude  The latitude of a given pixel (in degree).
    # param longitude The longitude of the given pixel (in degree).
    # param altitude  The altitude of the given pixel (in m)
    # param xyz       The xyz coordinates of the given pixel.
    # param geo_system The geodetic system.
    :return:
    """

    if geo_system == EarthModel.WGS84.name:
        a = WGS84.a
        e2 = WGS84.e2

    elif geo_system == EarthModel.GRS80.name:
        a = GRS80.a
        e2 = GRS80.e2

    else:
        raise ValueError(
            "Invalid geodetic system. Expected: {}. Got: {}".format(
                ", ".join(gs.value for gs in EarthModel), geo_system
            )
        )

    lat = degrees2radians(latitude)
    lon = degrees2radians(longitude)

    # x, y, z in meters
    sin_lat = np.sin(lat)
    cos_lat = np.cos(lat)
    N = a / np.sqrt(1 - e2 * sin_lat * sin_lat)
    x = (N + altitude) * cos_lat * np.cos(lon)
    y = (N + altitude) * cos_lat * np.sin(lon)
    z = ((1 - e2) * N + altitude) * sin_lat

    return XYZ(x, y, z)


def xyz_to_geo_wgs84(xyz: XYZ, geo_point: GeoPoint):
    """
    Convert cartesian XYZ coordinate into geodetic coordinate with specified geodetic system.

    :param xyz    The xyz coordinate of the given pixel.
    :param geo_point: The geodetic coordinate of the given pixel.
    :return:
    """

    x = xyz.x
    y = xyz.y
    z = xyz.z
    s = np.sqrt(x * x + y * y)
    theta = np.atan(z * WGS84.a / (s * WGS84.b))

    geo_point.lon = radians2degrees(np.atan(y / x))

    if geo_point.lon < 0.0 and y >= 0.0:
        geo_point.lon += 180.0

    elif geo_point.lon > 0.0 and y < 0.0:
        geo_point.lon -= 180.0

    geo_point.lat = np.atan(
        (z + WGS84.ep2 * WGS84.b * np.pow(np.sin(theta), 3))
        / radians2degrees(s - WGS84.e2 * WGS84.a * np.pow(np.cos(theta), 3))
    )


def polar2cartesian(latitude, longitude, radius) -> XYZ:
    lat_rad = degrees2radians(latitude)
    lon_rad = degrees2radians(longitude)
    sin_lat = np.sin(lat_rad)
    cos_lat = np.cos(lat_rad)
    return XYZ(
        x=radius * cos_lat * np.cos(lon_rad),
        y=radius * cos_lat * np.sin(lon_rad),
        z=radius * sin_lat,
    )