0
# Geodesic Operations
1
2
Geodesic operations compute the shortest path (great circle) between points on a sphere or ellipsoid. PyProj provides comprehensive geodesic calculations including distance, azimuth, area, and intermediate point computations using the GeographicLib algorithms for high accuracy.
3
4
## Capabilities
5
6
### Geod Creation
7
8
Create Geod objects for geodesic computations on specific ellipsoids with various initialization options.
9
10
```python { .api }
11
class Geod:
12
def __init__(
13
self,
14
initstring: str | None = None,
15
a: float | None = None,
16
b: float | None = None,
17
f: float | None = None,
18
es: float | None = None,
19
rf: float | None = None,
20
ellps: str | None = None,
21
sphere: bool = False,
22
**kwargs
23
) -> None:
24
"""
25
Create a Geod object for geodesic computations.
26
27
Args:
28
initstring: Initialization string (e.g., 'ellps=WGS84')
29
a: Semi-major axis of ellipsoid in meters
30
b: Semi-minor axis of ellipsoid in meters
31
f: Flattening of ellipsoid
32
es: Eccentricity squared of ellipsoid
33
rf: Reciprocal flattening of ellipsoid
34
ellps: Ellipsoid name (e.g., 'WGS84', 'GRS80')
35
sphere: Use spherical approximation
36
**kwargs: Additional ellipsoid parameters
37
38
Raises:
39
GeodError: If ellipsoid parameters are invalid or inconsistent
40
"""
41
42
@property
43
def a(self) -> float:
44
"""Semi-major axis of ellipsoid in meters."""
45
46
@property
47
def b(self) -> float:
48
"""Semi-minor axis of ellipsoid in meters."""
49
50
@property
51
def f(self) -> float:
52
"""Flattening of ellipsoid."""
53
54
@property
55
def es(self) -> float:
56
"""Eccentricity squared of ellipsoid."""
57
58
@property
59
def initstring(self) -> str:
60
"""Initialization string used to create Geod object."""
61
```
62
63
### Forward Geodesic Calculations
64
65
Compute positions and azimuths from starting points, initial bearings, and distances.
66
67
```python { .api }
68
class Geod:
69
def fwd(
70
self,
71
lons,
72
lats,
73
az,
74
dist,
75
radians: bool = False,
76
return_back_azimuth: bool = False
77
) -> tuple:
78
"""
79
Forward geodesic computation from points, azimuths, and distances.
80
81
Args:
82
lons: Longitude(s) of starting point(s) in degrees
83
lats: Latitude(s) of starting point(s) in degrees
84
az: Azimuth(s) in degrees from north (0-360)
85
dist: Distance(s) in meters
86
radians: Input/output coordinates in radians
87
return_back_azimuth: Return back azimuth at destination
88
89
Returns:
90
If return_back_azimuth is False: (end_lons, end_lats)
91
If return_back_azimuth is True: (end_lons, end_lats, back_azimuths)
92
93
Raises:
94
GeodError: If computation fails or inputs are invalid
95
"""
96
97
def fwd_intermediate(
98
self,
99
lon1: float,
100
lat1: float,
101
azi1: float,
102
dist: float,
103
npts: int,
104
del_s: float | None = None,
105
initial_idx: int = 0,
106
terminus_idx: int = 0,
107
flags: GeodIntermediateFlag = GeodIntermediateFlag.DEFAULT,
108
radians: bool = False,
109
return_back_azimuth: bool = False
110
) -> GeodIntermediateReturn:
111
"""
112
Forward geodesic with intermediate points.
113
114
Args:
115
lon1: Starting longitude in degrees
116
lat1: Starting latitude in degrees
117
azi1: Starting azimuth in degrees
118
dist: Total distance in meters
119
npts: Number of intermediate points
120
del_s: Distance between points (alternative to npts)
121
initial_idx: Starting index for point generation
122
terminus_idx: Ending index for point generation
123
flags: Flags for intermediate calculation behavior
124
radians: Input/output in radians
125
return_back_azimuth: Return back azimuths
126
127
Returns:
128
GeodIntermediateReturn object with points and metadata
129
130
Raises:
131
GeodError: If parameters are inconsistent or calculation fails
132
"""
133
```
134
135
### Inverse Geodesic Calculations
136
137
Compute distances, azimuths, and paths between pairs of points.
138
139
```python { .api }
140
class Geod:
141
def inv(
142
self,
143
lons1,
144
lats1,
145
lons2,
146
lats2,
147
radians: bool = False,
148
return_back_azimuth: bool = False
149
) -> tuple:
150
"""
151
Inverse geodesic computation between point pairs.
152
153
Args:
154
lons1: Longitude(s) of first point(s) in degrees
155
lats1: Latitude(s) of first point(s) in degrees
156
lons2: Longitude(s) of second point(s) in degrees
157
lats2: Latitude(s) of second point(s) in degrees
158
radians: Input coordinates in radians
159
return_back_azimuth: Return back azimuth from point 2 to point 1
160
161
Returns:
162
If return_back_azimuth is False: (forward_azimuths, distances)
163
If return_back_azimuth is True: (forward_azimuths, back_azimuths, distances)
164
165
Raises:
166
GeodError: If computation fails or coordinates are invalid
167
"""
168
169
def inv_intermediate(
170
self,
171
lon1: float,
172
lat1: float,
173
lon2: float,
174
lat2: float,
175
npts: int,
176
del_s: float | None = None,
177
initial_idx: int = 0,
178
terminus_idx: int = 0,
179
flags: GeodIntermediateFlag = GeodIntermediateFlag.DEFAULT,
180
radians: bool = False,
181
return_back_azimuth: bool = False
182
) -> GeodIntermediateReturn:
183
"""
184
Inverse geodesic with intermediate points.
185
186
Args:
187
lon1: First point longitude in degrees
188
lat1: First point latitude in degrees
189
lon2: Second point longitude in degrees
190
lat2: Second point latitude in degrees
191
npts: Number of intermediate points
192
del_s: Distance between points (alternative to npts)
193
initial_idx: Starting index for point generation
194
terminus_idx: Ending index for point generation
195
flags: Flags for intermediate calculation behavior
196
radians: Input/output in radians
197
return_back_azimuth: Return back azimuths
198
199
Returns:
200
GeodIntermediateReturn object with points and metadata
201
202
Raises:
203
GeodError: If points are identical or calculation fails
204
"""
205
206
def npts(
207
self,
208
lon1: float,
209
lat1: float,
210
lon2: float,
211
lat2: float,
212
npts: int,
213
radians: bool = False
214
) -> list[tuple[float, float]]:
215
"""
216
Compute intermediate points along geodesic.
217
218
Args:
219
lon1: First point longitude in degrees
220
lat1: First point latitude in degrees
221
lon2: Second point longitude in degrees
222
lat2: Second point latitude in degrees
223
npts: Number of intermediate points (excluding endpoints)
224
radians: Input/output coordinates in radians
225
226
Returns:
227
List of (longitude, latitude) tuples for intermediate points
228
229
Raises:
230
GeodError: If points are identical or npts is invalid
231
"""
232
```
233
234
### Line and Path Calculations
235
236
Compute lengths, distances, and cumulative measurements along paths defined by coordinate sequences.
237
238
```python { .api }
239
class Geod:
240
def line_length(
241
self,
242
lons: Any,
243
lats: Any,
244
radians: bool = False
245
) -> float:
246
"""
247
Calculate total length of line defined by coordinate sequence.
248
249
Args:
250
lons: Array-like of longitude coordinates
251
lats: Array-like of latitude coordinates
252
radians: Coordinates are in radians
253
254
Returns:
255
Total line length in meters
256
257
Raises:
258
GeodError: If coordinate arrays have different lengths or < 2 points
259
"""
260
261
def line_lengths(
262
self,
263
lons: Any,
264
lats: Any,
265
radians: bool = False
266
) -> Any:
267
"""
268
Calculate segment lengths along line defined by coordinate sequence.
269
270
Args:
271
lons: Array-like of longitude coordinates
272
lats: Array-like of latitude coordinates
273
radians: Coordinates are in radians
274
275
Returns:
276
Array of segment lengths in meters (length = n_points - 1)
277
278
Raises:
279
GeodError: If coordinate arrays have different lengths or < 2 points
280
"""
281
```
282
283
### Polygon Area and Perimeter
284
285
Compute area and perimeter of polygons defined by coordinate sequences using accurate geodesic methods.
286
287
```python { .api }
288
class Geod:
289
def polygon_area_perimeter(
290
self,
291
lons,
292
lats,
293
radians: bool = False
294
) -> tuple[float, float]:
295
"""
296
Calculate area and perimeter of polygon.
297
298
Args:
299
lons: Array-like of polygon vertex longitudes
300
lats: Array-like of polygon vertex latitudes
301
radians: Coordinates are in radians
302
303
Returns:
304
Tuple of (area_square_meters, perimeter_meters)
305
306
Note:
307
Polygon is automatically closed if first != last point.
308
Positive area indicates counter-clockwise vertex order.
309
Negative area indicates clockwise vertex order.
310
311
Raises:
312
GeodError: If coordinate arrays have different lengths or < 3 points
313
"""
314
```
315
316
### Geometry Operations
317
318
Work with geometry objects (requires Shapely) for convenient geodesic calculations on geometric features.
319
320
```python { .api }
321
class Geod:
322
def geometry_length(
323
self,
324
geometry,
325
radians: bool = False
326
) -> float:
327
"""
328
Calculate length of geometry object.
329
330
Args:
331
geometry: Shapely geometry object (LineString, MultiLineString, etc.)
332
radians: Coordinates are in radians
333
334
Returns:
335
Total length in meters
336
337
Note:
338
Requires Shapely to be installed.
339
Supports LineString, MultiLineString, LinearRing.
340
341
Raises:
342
ImportError: If Shapely is not installed
343
GeodError: If geometry type is not supported
344
"""
345
346
def geometry_area_perimeter(
347
self,
348
geometry,
349
radians: bool = False
350
) -> tuple[float, float]:
351
"""
352
Calculate area and perimeter of geometry object.
353
354
Args:
355
geometry: Shapely geometry object (Polygon, MultiPolygon, etc.)
356
radians: Coordinates are in radians
357
358
Returns:
359
Tuple of (area_square_meters, perimeter_meters)
360
361
Note:
362
Requires Shapely to be installed.
363
Supports Polygon, MultiPolygon.
364
365
Raises:
366
ImportError: If Shapely is not installed
367
GeodError: If geometry type is not supported
368
"""
369
```
370
371
## Usage Examples
372
373
### Basic Geodesic Calculations
374
375
```python
376
from pyproj import Geod
377
378
# Create Geod object for WGS84 ellipsoid
379
geod = Geod(ellps='WGS84')
380
381
# Forward calculation: from point + bearing + distance -> destination
382
lon1, lat1 = -74.0, 40.7 # New York
383
azimuth = 45.0 # Northeast
384
distance = 100000 # 100 km
385
386
lon2, lat2 = geod.fwd(lon1, lat1, azimuth, distance)
387
print(f"Destination: {lon2:.6f}, {lat2:.6f}")
388
389
# Inverse calculation: between two points -> azimuth + distance
390
lon2, lat2 = -73.0, 41.0 # Another point
391
azimuth, back_azimuth, distance = geod.inv(
392
lon1, lat1, lon2, lat2,
393
return_back_azimuth=True
394
)
395
print(f"Distance: {distance:.2f} m, Azimuth: {azimuth:.2f}°")
396
```
397
398
### Path Analysis
399
400
```python
401
from pyproj import Geod
402
import numpy as np
403
404
geod = Geod(ellps='WGS84')
405
406
# Define a path (e.g., flight route)
407
lons = np.array([-74.0, -70.0, -65.0, -60.0]) # New York to Atlantic
408
lats = np.array([40.7, 42.0, 43.0, 44.0])
409
410
# Calculate total path length
411
total_length = geod.line_length(lons, lats)
412
print(f"Total path length: {total_length/1000:.2f} km")
413
414
# Calculate segment lengths
415
segments = geod.line_lengths(lons, lats)
416
print(f"Segment lengths: {segments/1000}") # Convert to km
417
418
# Get intermediate points along path
419
intermediate = geod.npts(lons[0], lats[0], lons[-1], lats[-1], npts=10)
420
print(f"Intermediate points: {len(intermediate)}")
421
```
422
423
### Polygon Calculations
424
425
```python
426
from pyproj import Geod
427
import numpy as np
428
429
geod = Geod(ellps='WGS84')
430
431
# Define a polygon (e.g., state boundary)
432
# Triangle around New York area
433
lons = np.array([-74.0, -73.0, -75.0, -74.0]) # Close polygon
434
lats = np.array([40.7, 41.5, 41.0, 40.7])
435
436
# Calculate area and perimeter
437
area, perimeter = geod.polygon_area_perimeter(lons, lats)
438
439
print(f"Area: {abs(area)/1e6:.2f} km²") # Convert to km²
440
print(f"Perimeter: {perimeter/1000:.2f} km") # Convert to km
441
442
# Check if polygon is clockwise or counter-clockwise
443
if area > 0:
444
print("Counter-clockwise orientation")
445
else:
446
print("Clockwise orientation")
447
```
448
449
### Working with Different Ellipsoids
450
451
```python
452
from pyproj import Geod
453
454
# Different ellipsoids for regional accuracy
455
wgs84_geod = Geod(ellps='WGS84')
456
grs80_geod = Geod(ellps='GRS80')
457
clarke_geod = Geod(ellps='clrk66') # Clarke 1866
458
459
# Custom ellipsoid parameters
460
custom_geod = Geod(a=6378137.0, f=1/298.257223563)
461
462
# Compare distance calculations
463
lon1, lat1 = -100.0, 45.0
464
lon2, lat2 = -95.0, 50.0
465
466
distances = []
467
for name, geod in [('WGS84', wgs84_geod), ('GRS80', grs80_geod), ('Clarke', clarke_geod)]:
468
_, _, dist = geod.inv(lon1, lat1, lon2, lat2)
469
distances.append((name, dist))
470
print(f"{name}: {dist:.3f} m")
471
```
472
473
### Intermediate Points and Great Circle Routes
474
475
```python
476
from pyproj import Geod, GeodIntermediateFlag
477
478
geod = Geod(ellps='WGS84')
479
480
# Great circle route from London to Tokyo
481
lon1, lat1 = 0.0, 51.5 # London
482
lon2, lat2 = 139.7, 35.7 # Tokyo
483
484
# Generate route with specific number of points
485
route_points = geod.npts(lon1, lat1, lon2, lat2, npts=50)
486
487
# Get detailed intermediate information
488
intermediate = geod.inv_intermediate(
489
lon1, lat1, lon2, lat2,
490
npts=20,
491
flags=GeodIntermediateFlag.AZIS_KEEP, # Keep azimuth info
492
return_back_azimuth=True
493
)
494
495
print(f"Total points in route: {len(intermediate.lons)}")
496
print(f"Distance intervals: {intermediate.del_s:.2f} m")
497
498
# Forward intermediate calculation
499
forward_route = geod.fwd_intermediate(
500
lon1, lat1,
501
45.0, # Initial azimuth (northeast)
502
1000000, # 1000 km
503
npts=10
504
)
505
```
506
507
### Working with Shapely Geometries
508
509
```python
510
from pyproj import Geod
511
from shapely.geometry import LineString, Polygon
512
513
geod = Geod(ellps='WGS84')
514
515
# LineString length calculation
516
line = LineString([(-74.0, 40.7), (-73.0, 41.0), (-72.0, 41.5)])
517
length = geod.geometry_length(line)
518
print(f"Line length: {length/1000:.2f} km")
519
520
# Polygon area calculation
521
coords = [(-74.0, 40.5), (-73.0, 40.5), (-73.0, 41.0), (-74.0, 41.0), (-74.0, 40.5)]
522
polygon = Polygon(coords)
523
area, perimeter = geod.geometry_area_perimeter(polygon)
524
print(f"Polygon area: {area/1e6:.2f} km²")
525
print(f"Polygon perimeter: {perimeter/1000:.2f} km")
526
```
527
528
## Types
529
530
```python { .api }
531
# Geodesic intermediate calculation results
532
class GeodIntermediateReturn:
533
"""Return object for intermediate geodesic calculations."""
534
npts: int # Number of points generated
535
del_s: float # Distance between points in meters
536
lons: list[float] # Longitude coordinates
537
lats: list[float] # Latitude coordinates
538
azis: list[float] | None # Azimuth values (if requested)
539
540
# Flags for intermediate calculations
541
class GeodIntermediateFlag(Enum):
542
"""Flags controlling intermediate point calculation behavior."""
543
DEFAULT = 0
544
NPTS_ROUND = 1 # Round number of points to nearest integer
545
NPTS_CEIL = 2 # Round number of points up
546
NPTS_TRUNC = 3 # Truncate number of points
547
DEL_S_RECALC = 4 # Recalculate distance between points
548
DEL_S_NO_RECALC = 5 # Don't recalculate distance
549
AZIS_DISCARD = 6 # Don't store azimuth information
550
AZIS_KEEP = 7 # Store azimuth information
551
552
# Global geodesic constants
553
geodesic_version_str: str # Version of GeographicLib
554
pj_ellps: dict[str, dict] # Dictionary of ellipsoid parameters
555
```