0
# Geodetic Calculations
1
2
Geographic and coordinate system calculations for seismological applications including distance, azimuth, and coordinate transformations using spherical and ellipsoidal Earth models. These functions are essential for earthquake location, station spacing analysis, and seismological coordinate systems.
3
4
## Capabilities
5
6
### Distance and Azimuth Calculations
7
8
Precise geographic calculations using both spherical and ellipsoidal Earth models for seismological applications.
9
10
```python { .api }
11
# Import from obspy.geodetics
12
def gps2dist_azimuth(lat1: float, lon1: float, lat2: float, lon2: float,
13
a: float = 6378137.0, f: float = 1/298.257223563) -> tuple:
14
"""
15
Calculate distance and azimuth between two geographic points using ellipsoid.
16
17
Args:
18
lat1: Latitude of first point in degrees
19
lon1: Longitude of first point in degrees
20
lat2: Latitude of second point in degrees
21
lon2: Longitude of second point in degrees
22
a: Semi-major axis of ellipsoid in meters (WGS84 default)
23
f: Flattening of ellipsoid (WGS84 default)
24
25
Returns:
26
Tuple of (distance_in_meters, azimuth_deg, back_azimuth_deg)
27
28
Notes:
29
Uses Vincenty's formula for high precision on ellipsoid.
30
Azimuth is from point 1 to point 2 (0° = North, 90° = East).
31
Back-azimuth is from point 2 to point 1.
32
"""
33
34
def calc_vincenty_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
35
a: float = 6378137.0, f: float = 1/298.257223563,
36
max_iter: int = 200, tol: float = 1e-12) -> dict:
37
"""
38
Vincenty's inverse formula for ellipsoid distance calculation.
39
40
Args:
41
lat1: First point latitude in degrees
42
lon1: First point longitude in degrees
43
lat2: Second point latitude in degrees
44
lon2: Second point longitude in degrees
45
a: Semi-major axis in meters
46
f: Flattening parameter
47
max_iter: Maximum iterations for convergence
48
tol: Convergence tolerance
49
50
Returns:
51
Dictionary with keys:
52
- 'distance': Distance in meters
53
- 'azimuth': Forward azimuth in degrees
54
- 'reverse_azimuth': Reverse azimuth in degrees
55
- 'iterations': Number of iterations used
56
"""
57
58
def locations2degrees(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
59
"""
60
Calculate angular distance between two points on sphere.
61
62
Args:
63
lat1: First point latitude in degrees
64
lon1: First point longitude in degrees
65
lat2: Second point latitude in degrees
66
lon2: Second point longitude in degrees
67
68
Returns:
69
Angular distance in degrees
70
71
Notes:
72
Uses great circle distance on unit sphere.
73
Suitable for approximate seismological calculations.
74
"""
75
```
76
77
### Unit Conversions
78
79
Conversions between different distance units commonly used in seismology.
80
81
```python { .api }
82
def degrees2kilometers(degrees: float, radius: float = 6371.0) -> float:
83
"""
84
Convert degrees to kilometers using spherical Earth.
85
86
Args:
87
degrees: Angular distance in degrees
88
radius: Earth radius in kilometers (mean radius default)
89
90
Returns:
91
Distance in kilometers
92
"""
93
94
def kilometers2degrees(kilometers: float, radius: float = 6371.0) -> float:
95
"""
96
Convert kilometers to degrees using spherical Earth.
97
98
Args:
99
kilometers: Distance in kilometers
100
radius: Earth radius in kilometers
101
102
Returns:
103
Angular distance in degrees
104
"""
105
106
def kilometer2degrees(kilometers: float, radius: float = 6371.0) -> float:
107
"""
108
Convert kilometers to degrees (alias for backwards compatibility).
109
110
Args:
111
kilometers: Distance in kilometers
112
radius: Earth radius in kilometers
113
114
Returns:
115
Angular distance in degrees
116
"""
117
```
118
119
### Geographic Boundary Checking
120
121
Functions for determining if coordinates fall within specified geographic regions.
122
123
```python { .api }
124
def inside_geobounds(longitude: float, latitude: float,
125
min_longitude: float, max_longitude: float,
126
min_latitude: float, max_latitude: float) -> bool:
127
"""
128
Check if point lies within geographic boundaries.
129
130
Args:
131
longitude: Point longitude in degrees
132
latitude: Point latitude in degrees
133
min_longitude: Western boundary in degrees
134
max_longitude: Eastern boundary in degrees
135
min_latitude: Southern boundary in degrees
136
max_latitude: Northern boundary in degrees
137
138
Returns:
139
True if point is inside boundaries
140
141
Notes:
142
Handles longitude wrapping around ±180°.
143
Accounts for regions crossing the International Date Line.
144
"""
145
```
146
147
### Flinn-Engdahl Regionalization
148
149
Seismological regionalization system for identifying geographic regions by coordinates.
150
151
```python { .api }
152
class FlinnEngdahl:
153
"""
154
Flinn-Engdahl seismological regionalization.
155
156
Provides standardized geographic region identification used in
157
earthquake catalogs and seismological studies.
158
"""
159
160
def __init__(self):
161
"""Initialize Flinn-Engdahl regionalization system."""
162
163
def get_region(self, longitude: float, latitude: float) -> str:
164
"""
165
Get region name for geographic coordinates.
166
167
Args:
168
longitude: Longitude in degrees (-180 to 180)
169
latitude: Latitude in degrees (-90 to 90)
170
171
Returns:
172
Region name string (e.g., "Southern California", "Japan")
173
"""
174
175
def get_region_number(self, longitude: float, latitude: float) -> int:
176
"""
177
Get numeric region code for coordinates.
178
179
Args:
180
longitude: Longitude in degrees
181
latitude: Latitude in degrees
182
183
Returns:
184
Flinn-Engdahl region number (1-757)
185
"""
186
```
187
188
## Usage Examples
189
190
### Basic Distance and Azimuth Calculations
191
192
```python
193
from obspy.geodetics import gps2dist_azimuth, locations2degrees
194
195
# Calculate distance between two seismic stations
196
sta1_lat, sta1_lon = 34.0, -118.0 # Los Angeles
197
sta2_lat, sta2_lon = 37.8, -122.4 # San Francisco
198
199
# High-precision ellipsoid calculation
200
dist_m, az, baz = gps2dist_azimuth(sta1_lat, sta1_lon, sta2_lat, sta2_lon)
201
print(f"Distance: {dist_m/1000:.2f} km")
202
print(f"Azimuth LA→SF: {az:.1f}°")
203
print(f"Back-azimuth SF→LA: {baz:.1f}°")
204
205
# Quick spherical approximation
206
dist_deg = locations2degrees(sta1_lat, sta1_lon, sta2_lat, sta2_lon)
207
print(f"Angular distance: {dist_deg:.3f}°")
208
```
209
210
### Event-Station Distance Analysis
211
212
```python
213
from obspy.geodetics import gps2dist_azimuth, degrees2kilometers
214
from obspy import UTCDateTime
215
216
# Earthquake and station information
217
event_lat, event_lon = 35.0, 140.0 # Japan earthquake
218
event_time = UTCDateTime("2023-01-15T10:30:00")
219
220
stations = {
221
"Tokyo": (35.7, 139.7),
222
"Seoul": (37.6, 126.9),
223
"Beijing": (39.9, 116.4),
224
"Taipei": (25.0, 121.5),
225
"Manila": (14.6, 121.0)
226
}
227
228
print("Station distances and azimuths from earthquake:")
229
print("Station".ljust(10), "Distance (km)".ljust(12), "Azimuth (°)".ljust(12), "Category")
230
print("-" * 50)
231
232
for station, (sta_lat, sta_lon) in stations.items():
233
dist_m, az, baz = gps2dist_azimuth(event_lat, event_lon, sta_lat, sta_lon)
234
dist_km = dist_m / 1000.0
235
236
# Classify distance
237
if dist_km < 100:
238
category = "Local"
239
elif dist_km < 1000:
240
category = "Regional"
241
else:
242
category = "Teleseismic"
243
244
print(f"{station:<10} {dist_km:>8.1f} {az:>8.1f} {category}")
245
```
246
247
### Array Geometry Analysis
248
249
```python
250
from obspy.geodetics import gps2dist_azimuth
251
import numpy as np
252
253
# Seismic array station coordinates
254
array_center = (40.0, -105.0) # Colorado
255
stations = {
256
"A01": (40.00, -105.00), # Center station
257
"A02": (40.01, -105.00), # North
258
"A03": (40.00, -105.01), # West
259
"A04": (39.99, -105.00), # South
260
"A05": (40.00, -104.99), # East
261
}
262
263
print("Array geometry analysis:")
264
print("Station Distance (m) Azimuth (°)")
265
print("-" * 35)
266
267
distances = []
268
azimuths = []
269
270
for station, (lat, lon) in stations.items():
271
if station == "A01": # Skip center station
272
continue
273
274
dist_m, az, baz = gps2dist_azimuth(array_center[0], array_center[1], lat, lon)
275
distances.append(dist_m)
276
azimuths.append(az)
277
278
print(f"{station} {dist_m:>8.1f} {az:>7.1f}")
279
280
# Array statistics
281
print(f"\nArray statistics:")
282
print(f"Mean aperture: {np.mean(distances):.1f} m")
283
print(f"Aperture range: {np.min(distances):.1f} - {np.max(distances):.1f} m")
284
print(f"Azimuth coverage: {np.min(azimuths):.1f}° - {np.max(azimuths):.1f}°")
285
```
286
287
### Unit Conversion Utilities
288
289
```python
290
from obspy.geodetics import degrees2kilometers, kilometers2degrees
291
292
# Convert between angular and linear distance units
293
angular_distances = [1.0, 5.0, 10.0, 30.0, 90.0, 180.0] # degrees
294
295
print("Distance conversions:")
296
print("Degrees Kilometers")
297
print("-" * 22)
298
299
for deg in angular_distances:
300
km = degrees2kilometers(deg)
301
print(f"{deg:>7.1f} {km:>10.1f}")
302
303
print("\nReverse conversions:")
304
print("Kilometers Degrees")
305
print("-" * 22)
306
307
linear_distances = [100, 500, 1000, 5000, 10000] # km
308
for km in linear_distances:
309
deg = kilometers2degrees(km)
310
print(f"{km:>10.0f} {deg:>7.3f}")
311
312
# Earth radius variations
313
print("\nDistance conversions with different Earth radii:")
314
radii = {"Mean": 6371.0, "Equatorial": 6378.1, "Polar": 6356.8}
315
316
for name, radius in radii.items():
317
km = degrees2kilometers(10.0, radius=radius)
318
print(f"10° using {name} radius ({radius} km): {km:.1f} km")
319
```
320
321
### Geographic Region Classification
322
323
```python
324
from obspy.geodetics import FlinnEngdahl
325
326
# Initialize Flinn-Engdahl regionalization
327
fe = FlinnEngdahl()
328
329
# Classify earthquake locations
330
earthquakes = [
331
("M6.5 California", 34.0, -118.0),
332
("M7.2 Japan", 35.0, 140.0),
333
("M8.0 Chile", -30.0, -71.0),
334
("M6.8 Turkey", 39.0, 35.0),
335
("M7.5 Alaska", 61.0, -150.0),
336
("M6.9 Indonesia", -6.0, 106.0)
337
]
338
339
print("Earthquake regional classification:")
340
print("Event".ljust(20), "Coordinates".ljust(15), "Flinn-Engdahl Region")
341
print("-" * 70)
342
343
for event_name, lat, lon in earthquakes:
344
region = fe.get_region(lon, lat) # Note: longitude first
345
region_num = fe.get_region_number(lon, lat)
346
coords = f"({lat:.1f}, {lon:.1f})"
347
348
print(f"{event_name:<20} {coords:<15} {region} ({region_num})")
349
```
350
351
### Station Network Spacing Analysis
352
353
```python
354
from obspy.geodetics import gps2dist_azimuth
355
import numpy as np
356
357
# Global seismographic network stations (subset)
358
gsn_stations = {
359
"ANMO": (34.9459, -106.4572), # New Mexico, USA
360
"ANTO": (39.8682, 32.7934), # Turkey
361
"ASAR": (23.2729, 5.6309), # Algeria
362
"ASCN": (-7.9327, -14.3601), # Ascension Island
363
"BGCA": (52.3858, -113.3072), # Canada
364
"BOCO": (4.5869, -74.0432), # Colombia
365
"COCO": (-12.1901, 96.8349), # Cocos Islands
366
"CTAO": (-20.0882, 146.2545), # Australia
367
"FUNA": (-8.5259, 179.1966), # Tuvalu
368
"GUMO": (13.5893, 144.8684) # Guam
369
}
370
371
print("Global network station spacing analysis:")
372
print("-" * 50)
373
374
# Calculate minimum distances between stations
375
min_distances = {}
376
for sta1, (lat1, lon1) in gsn_stations.items():
377
min_dist = float('inf')
378
closest_station = ""
379
380
for sta2, (lat2, lon2) in gsn_stations.items():
381
if sta1 == sta2:
382
continue
383
384
dist_m, _, _ = gps2dist_azimuth(lat1, lon1, lat2, lon2)
385
if dist_m < min_dist:
386
min_dist = dist_m
387
closest_station = sta2
388
389
min_distances[sta1] = (min_dist/1000.0, closest_station)
390
391
# Sort by distance
392
sorted_stations = sorted(min_distances.items(), key=lambda x: x[1][0])
393
394
print("Station Closest Station Distance (km)")
395
print("-" * 45)
396
for station, (dist, closest) in sorted_stations:
397
print(f"{station:<10} {closest:<15} {dist:>8.1f}")
398
399
# Network statistics
400
distances = [dist for dist, _ in min_distances.values()]
401
print(f"\nNetwork spacing statistics:")
402
print(f"Mean minimum distance: {np.mean(distances):.1f} km")
403
print(f"Median minimum distance: {np.median(distances):.1f} km")
404
print(f"Distance range: {np.min(distances):.1f} - {np.max(distances):.1f} km")
405
```
406
407
### Boundary and Region Checking
408
409
```python
410
from obspy.geodetics import inside_geobounds
411
412
# Define study regions
413
regions = {
414
"California": {"min_lon": -125.0, "max_lon": -114.0,
415
"min_lat": 32.0, "max_lat": 42.0},
416
"Japan": {"min_lon": 129.0, "max_lon": 146.0,
417
"min_lat": 30.0, "max_lat": 46.0},
418
"Mediterranean": {"min_lon": -10.0, "max_lon": 45.0,
419
"min_lat": 30.0, "max_lat": 48.0}
420
}
421
422
# Test earthquake locations
423
test_events = [
424
("Northridge", 34.2, -118.5),
425
("Kobe", 34.6, 135.0),
426
("L'Aquila", 42.3, 13.4),
427
("Christchurch", -43.5, 172.6),
428
("Mexico City", 19.4, -99.1)
429
]
430
431
print("Event location classification:")
432
print("Event".ljust(15), "Coordinates".ljust(18), "Regions")
433
print("-" * 50)
434
435
for event_name, lat, lon in test_events:
436
coords = f"({lat:.1f}, {lon:.1f})"
437
regions_found = []
438
439
for region_name, bounds in regions.items():
440
if inside_geobounds(lon, lat,
441
bounds["min_lon"], bounds["max_lon"],
442
bounds["min_lat"], bounds["max_lat"]):
443
regions_found.append(region_name)
444
445
regions_str = ", ".join(regions_found) if regions_found else "None"
446
print(f"{event_name:<15} {coords:<18} {regions_str}")
447
```
448
449
## Types
450
451
```python { .api }
452
# Coordinate point structure
453
Point = {
454
'latitude': float, # Latitude in degrees
455
'longitude': float, # Longitude in degrees
456
'elevation': float # Elevation in meters (optional)
457
}
458
459
# Distance calculation result
460
DistanceResult = {
461
'distance_m': float, # Distance in meters
462
'distance_km': float, # Distance in kilometers
463
'distance_deg': float, # Angular distance in degrees
464
'azimuth': float, # Forward azimuth in degrees
465
'back_azimuth': float # Back azimuth in degrees
466
}
467
468
# Geographic boundary definition
469
GeoBounds = {
470
'min_latitude': float, # Southern boundary
471
'max_latitude': float, # Northern boundary
472
'min_longitude': float, # Western boundary
473
'max_longitude': float # Eastern boundary
474
}
475
476
# Flinn-Engdahl region information
477
RegionInfo = {
478
'number': int, # Region number (1-757)
479
'name': str, # Region name
480
'type': str # Region type (geographic/seismic)
481
}
482
```