0
# Geodesic Calculations
1

2
Great circle computations, distance calculations, and geodesic operations on ellipsoidal Earth models. Essential for navigation, surveying, and geographic analysis applications that require accurate distance and bearing calculations.
3

4
## Capabilities
5

6
### Geodesic Class
7

8
Core class for geodesic calculations on ellipsoidal Earth models using Proj geodesic functions.
9

10
```python { .api }
11
class Geodesic:
12
    """
13
    Define an ellipsoid on which to solve geodesic problems.
14
    
15
    Parameters:
16
    - radius: float, equatorial radius in metres (default: WGS84 semimajor axis 6378137.0)
17
    - flattening: float, flattening of ellipsoid (default: WGS84 flattening 1/298.257223563)
18
    """
19
    def __init__(self, radius=6378137.0, flattening=1/298.257223563): ...
20
    
21
    def direct(self, points, azimuths, distances):
22
        """
23
        Solve direct geodesic problem (distance and azimuth to endpoint).
24
        
25
        Parameters:
26
        - points: array_like, shape=(n or 1, 2), starting longitude-latitude points
27
        - azimuths: float or array_like, azimuth values in degrees
28
        - distances: float or array_like, distance values in metres
29
        
30
        Returns:
31
        numpy.ndarray, shape=(n, 3), longitudes, latitudes, and forward azimuths of endpoints
32
        """
33
    
34
    def inverse(self, points, endpoints):
35
        """
36
        Solve inverse geodesic problem (distance and azimuth between points).
37
        
38
        Parameters:
39
        - points: array_like, shape=(n or 1, 2), starting longitude-latitude points
40
        - endpoints: array_like, shape=(n or 1, 2), ending longitude-latitude points
41
        
42
        Returns:
43
        numpy.ndarray, shape=(n, 3), distances and forward azimuths of start/end points
44
        """
45
    
46
    def circle(self, lon, lat, radius, n_samples=180, endpoint=False):
47
        """
48
        Create geodesic circle of given radius at specified point.
49
        
50
        Parameters:
51
        - lon: float, longitude coordinate of center
52
        - lat: float, latitude coordinate of center  
53
        - radius: float, radius of circle in metres
54
        - n_samples: int, number of sample points on circle (default 180)
55
        - endpoint: bool, whether to repeat endpoint (default False)
56
        
57
        Returns:
58
        numpy.ndarray, shape=(n_samples, 2), longitude-latitude points on circle
59
        """
60
    
61
    def geometry_length(self, geometry):
62
        """
63
        Calculate physical length of Shapely geometry on ellipsoid.
64
        
65
        Parameters:
66
        - geometry: shapely.geometry.BaseGeometry, geometry in spherical coordinates
67
        
68
        Returns:
69
        float, length in metres
70
        """
71
```
72

73
## Usage Examples
74

75
### Basic Distance and Bearing Calculations
76

77
```python
78
from cartopy.geodesic import Geodesic
79
import numpy as np
80

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

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

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

95
print(f"Distance: {distance_m/1000:.0f} km")
96
print(f"Initial bearing: {initial_bearing:.1f}°")
97
print(f"Final bearing: {final_bearing:.1f}°")
98
```
99

100
### Direct Geodesic Problem
101

102
```python
103
from cartopy.geodesic import Geodesic
104

105
geod = Geodesic()
106

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

112
# Find endpoint
113
endpoint = geod.direct(start_point, bearing, distance)
114
end_lon = endpoint[0, 0]
115
end_lat = endpoint[0, 1]
116
final_bearing = endpoint[0, 2]
117

118
print(f"Endpoint: {end_lon:.4f}°, {end_lat:.4f}°")
119
print(f"Final bearing: {final_bearing:.1f}°")
120
```
121

122
### Creating Geodesic Circles
123

124
```python
125
import matplotlib.pyplot as plt
126
import cartopy.crs as ccrs
127
from cartopy.geodesic import Geodesic
128

129
# Create geodesic calculator
130
geod = Geodesic()
131

132
# Create circles around major cities
133
cities = {
134
    'London': [0.0, 51.5],
135
    'New York': [-74.0, 40.7],
136
    'Tokyo': [139.7, 35.7]
137
}
138

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

142
for city_name, coords in cities.items():
143
    # Create 1000 km radius circle
144
    circle_points = geod.circle(coords[0], coords[1], 1000000, n_samples=100)
145
    
146
    # Plot circle
147
    ax.plot(circle_points[:, 0], circle_points[:, 1], 
148
            transform=ccrs.PlateCarree(), label=f'{city_name} (1000km)')
149
    
150
    # Mark city center
151
    ax.plot(coords[0], coords[1], 'ro', transform=ccrs.PlateCarree(), 
152
            markersize=8)
153

154
ax.coastlines()
155
ax.set_global()
156
ax.legend()
157
plt.title('Geodesic Circles Around Major Cities')
158
plt.show()
159
```
160

161
### Multiple Point Calculations with Broadcasting
162

163
```python
164
from cartopy.geodesic import Geodesic
165
import numpy as np
166

167
geod = Geodesic()
168

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

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

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

188
### Geometry Length Calculations
189

190
```python
191
from cartopy.geodesic import Geodesic
192
from shapely.geometry import LineString, Polygon
193
import numpy as np
194

195
geod = Geodesic()
196

197
# Calculate length of a route
198
route_coords = [
199
    [0.0, 51.5],      # London
200
    [2.3, 48.9],      # Paris
201
    [13.4, 52.5],     # Berlin
202
    [16.4, 48.2],     # Vienna
203
    [14.5, 46.0]      # Ljubljana
204
]
205

206
# Create LineString geometry
207
route = LineString(route_coords)
208

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

213
# Calculate polygon perimeter
214
# Create a triangle between three cities
215
triangle_coords = [
216
    [0.0, 51.5],      # London
217
    [2.3, 48.9],      # Paris  
218
    [4.9, 52.4],      # Amsterdam
219
    [0.0, 51.5]       # Back to London
220
]
221

222
triangle = Polygon(triangle_coords)
223
perimeter_km = geod.geometry_length(triangle) / 1000
224
print(f"Triangle perimeter: {perimeter_km:.0f} km")
225
```
226

227
### Custom Ellipsoid Calculations
228

229
```python
230
from cartopy.geodesic import Geodesic
231

232
# Create geodesic calculator for different ellipsoids
233

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

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

240
# Calculate same distance on different ellipsoids
241
start = [0.0, 0.0]
242
end = [1.0, 1.0]
243

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

248
print(f"WGS84 distance: {wgs84_dist:.2f} m")
249
print(f"Sphere distance: {sphere_dist:.2f} m") 
250
print(f"GRS80 distance: {grs80_dist:.2f} m")
251
```
252

253
### Great Circle Paths for Plotting
254

255
```python
256
import matplotlib.pyplot as plt
257
import cartopy.crs as ccrs
258
from cartopy.geodesic import Geodesic
259
import numpy as np
260

261
geod = Geodesic()
262

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

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

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

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

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

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

291
ax.coastlines()
292
ax.set_global()
293
ax.legend()
294
plt.title('Great Circle Path: California to Japan')
295
plt.show()
296
```
297

298
## Mathematical Background
299

300
The geodesic calculations in Cartopy use the algorithms developed by Charles Karney, which provide accurate solutions to geodesic problems on ellipsoids. These methods are implemented via the Proj library's geodesic functions.
301

302
### Key Concepts:
303

304
- **Direct Problem**: Given a starting point, initial bearing, and distance, find the endpoint
305
- **Inverse Problem**: Given two points, find the distance and bearings between them
306
- **Geodesic**: The shortest path between two points on an ellipsoid (generalization of great circles)
307
- **Azimuth**: The angle measured clockwise from north to the geodesic direction
308

309
### Accuracy:
310
The algorithms provide accuracy better than 15 nanometers for distances up to 20,000 km on the WGS84 ellipsoid.
311

312
## Constants
313

314
```python { .api }
315
# Default WGS84 parameters
316
WGS84_RADIUS: float = 6378137.0  # metres
317
WGS84_FLATTENING: float = 1/298.257223563
318
```