Geodesic Calculations

class Geodesic:
    """
    Define an ellipsoid on which to solve geodesic problems.
    
    Parameters:
    - radius: float, equatorial radius in metres (default: WGS84 semimajor axis 6378137.0)
    - flattening: float, flattening of ellipsoid (default: WGS84 flattening 1/298.257223563)
    """
    def __init__(self, radius=6378137.0, flattening=1/298.257223563): ...
    
    def direct(self, points, azimuths, distances):
        """
        Solve direct geodesic problem (distance and azimuth to endpoint).
        
        Parameters:
        - points: array_like, shape=(n or 1, 2), starting longitude-latitude points
        - azimuths: float or array_like, azimuth values in degrees
        - distances: float or array_like, distance values in metres
        
        Returns:
        numpy.ndarray, shape=(n, 3), longitudes, latitudes, and forward azimuths of endpoints
        """
    
    def inverse(self, points, endpoints):
        """
        Solve inverse geodesic problem (distance and azimuth between points).
        
        Parameters:
        - points: array_like, shape=(n or 1, 2), starting longitude-latitude points
        - endpoints: array_like, shape=(n or 1, 2), ending longitude-latitude points
        
        Returns:
        numpy.ndarray, shape=(n, 3), distances and forward azimuths of start/end points
        """
    
    def circle(self, lon, lat, radius, n_samples=180, endpoint=False):
        """
        Create geodesic circle of given radius at specified point.
        
        Parameters:
        - lon: float, longitude coordinate of center
        - lat: float, latitude coordinate of center  
        - radius: float, radius of circle in metres
        - n_samples: int, number of sample points on circle (default 180)
        - endpoint: bool, whether to repeat endpoint (default False)
        
        Returns:
        numpy.ndarray, shape=(n_samples, 2), longitude-latitude points on circle
        """
    
    def geometry_length(self, geometry):
        """
        Calculate physical length of Shapely geometry on ellipsoid.
        
        Parameters:
        - geometry: shapely.geometry.BaseGeometry, geometry in spherical coordinates
        
        Returns:
        float, length in metres
        """

from cartopy.geodesic import Geodesic
import numpy as np

# Create geodesic calculator with WGS84 ellipsoid
geod = Geodesic()

# Calculate distance between two cities
# New York to London
ny_coords = [-74.0060, 40.7128]  # [lon, lat]
london_coords = [-0.1278, 51.5074]

# Inverse problem: find distance and bearing
result = geod.inverse(ny_coords, london_coords)
distance_m = result[0, 0]
initial_bearing = result[0, 1]
final_bearing = result[0, 2]

print(f"Distance: {distance_m/1000:.0f} km")
print(f"Initial bearing: {initial_bearing:.1f}°")
print(f"Final bearing: {final_bearing:.1f}°")

from cartopy.geodesic import Geodesic

geod = Geodesic()

# Start from a point and travel a specific distance and bearing
start_point = [0.0, 51.5]  # London [lon, lat]
bearing = 45.0  # Northeast
distance = 100000  # 100 km

# Find endpoint
endpoint = geod.direct(start_point, bearing, distance)
end_lon = endpoint[0, 0]
end_lat = endpoint[0, 1]
final_bearing = endpoint[0, 2]

print(f"Endpoint: {end_lon:.4f}°, {end_lat:.4f}°")
print(f"Final bearing: {final_bearing:.1f}°")

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.geodesic import Geodesic

# Create geodesic calculator
geod = Geodesic()

# Create circles around major cities
cities = {
    'London': [0.0, 51.5],
    'New York': [-74.0, 40.7],
    'Tokyo': [139.7, 35.7]
}

fig = plt.figure(figsize=(15, 10))
ax = plt.axes(projection=ccrs.PlateCarree())

for city_name, coords in cities.items():
    # Create 1000 km radius circle
    circle_points = geod.circle(coords[0], coords[1], 1000000, n_samples=100)
    
    # Plot circle
    ax.plot(circle_points[:, 0], circle_points[:, 1], 
            transform=ccrs.PlateCarree(), label=f'{city_name} (1000km)')
    
    # Mark city center
    ax.plot(coords[0], coords[1], 'ro', transform=ccrs.PlateCarree(), 
            markersize=8)

ax.coastlines()
ax.set_global()
ax.legend()
plt.title('Geodesic Circles Around Major Cities')
plt.show()

from cartopy.geodesic import Geodesic
import numpy as np

geod = Geodesic()

# Calculate distances from London to multiple cities
london = [0.0, 51.5]
cities = np.array([
    [-74.0, 40.7],    # New York  
    [139.7, 35.7],    # Tokyo
    [151.2, -33.9],   # Sydney
    [-43.2, -22.9]    # Rio de Janeiro
])

# Use broadcasting: single start point, multiple endpoints
results = geod.inverse(london, cities)
distances_km = results[:, 0] / 1000
bearings = results[:, 1]

city_names = ['New York', 'Tokyo', 'Sydney', 'Rio de Janeiro']
for i, city in enumerate(city_names):
    print(f"London to {city}: {distances_km[i]:.0f} km, bearing {bearings[i]:.1f}°")

from cartopy.geodesic import Geodesic
from shapely.geometry import LineString, Polygon
import numpy as np

geod = Geodesic()

# Calculate length of a route
route_coords = [
    [0.0, 51.5],      # London
    [2.3, 48.9],      # Paris
    [13.4, 52.5],     # Berlin
    [16.4, 48.2],     # Vienna
    [14.5, 46.0]      # Ljubljana
]

# Create LineString geometry
route = LineString(route_coords)

# Calculate total route length
route_length_km = geod.geometry_length(route) / 1000
print(f"Total route length: {route_length_km:.0f} km")

# Calculate polygon perimeter
# Create a triangle between three cities
triangle_coords = [
    [0.0, 51.5],      # London
    [2.3, 48.9],      # Paris  
    [4.9, 52.4],      # Amsterdam
    [0.0, 51.5]       # Back to London
]

triangle = Polygon(triangle_coords)
perimeter_km = geod.geometry_length(triangle) / 1000
print(f"Triangle perimeter: {perimeter_km:.0f} km")

from cartopy.geodesic import Geodesic

# Create geodesic calculator for different ellipsoids

# Sphere (flattening = 0)
sphere_geod = Geodesic(radius=6371000, flattening=0.0)

# GRS80 ellipsoid
grs80_geod = Geodesic(radius=6378137.0, flattening=1/298.257222101)

# Calculate same distance on different ellipsoids
start = [0.0, 0.0]
end = [1.0, 1.0]

wgs84_dist = Geodesic().inverse(start, end)[0, 0]
sphere_dist = sphere_geod.inverse(start, end)[0, 0]
grs80_dist = grs80_geod.inverse(start, end)[0, 0]

print(f"WGS84 distance: {wgs84_dist:.2f} m")
print(f"Sphere distance: {sphere_dist:.2f} m") 
print(f"GRS80 distance: {grs80_dist:.2f} m")

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.geodesic import Geodesic
import numpy as np

geod = Geodesic()

# Create great circle path between distant points
start_point = [-120.0, 35.0]  # California
end_point = [140.0, 35.0]     # Japan

# Calculate intermediate points along great circle
# Use multiple direct calculations
result = geod.inverse(start_point, end_point)
total_distance = result[0, 0]
initial_bearing = result[0, 1]

# Create points every 500 km along the great circle
distances = np.arange(0, total_distance, 500000)
path_points = geod.direct(start_point, initial_bearing, distances)

# Plot the great circle path
fig = plt.figure(figsize=(15, 8))
ax = plt.axes(projection=ccrs.PlateCarree())

# Plot path points
ax.plot(path_points[:, 0], path_points[:, 1], 'r-', 
        transform=ccrs.Geodetic(), linewidth=2, label='Great Circle')

# Mark start and end points
ax.plot(start_point[0], start_point[1], 'go', 
        transform=ccrs.PlateCarree(), markersize=10, label='Start')
ax.plot(end_point[0], end_point[1], 'ro', 
        transform=ccrs.PlateCarree(), markersize=10, label='End')

ax.coastlines()
ax.set_global()
ax.legend()
plt.title('Great Circle Path: California to Japan')
plt.show()

# Default WGS84 parameters
WGS84_RADIUS: float = 6378137.0  # metres
WGS84_FLATTENING: float = 1/298.257223563