tessl/pypi-obspy

ObsPy is a Python toolbox for seismology providing parsers for seismological data formats, clients for data centers, and signal processing routines for seismological time series analysis.

—

Pending

Overview

Eval results

Files

Geodetic Calculations

Name: tessl/pypi-obspy
Author: tessl

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.

Capabilities

Distance and Azimuth Calculations

Precise geographic calculations using both spherical and ellipsoidal Earth models for seismological applications.

# Import from obspy.geodetics
def gps2dist_azimuth(lat1: float, lon1: float, lat2: float, lon2: float,
                    a: float = 6378137.0, f: float = 1/298.257223563) -> tuple:
    """
    Calculate distance and azimuth between two geographic points using ellipsoid.
    
    Args:
        lat1: Latitude of first point in degrees
        lon1: Longitude of first point in degrees  
        lat2: Latitude of second point in degrees
        lon2: Longitude of second point in degrees
        a: Semi-major axis of ellipsoid in meters (WGS84 default)
        f: Flattening of ellipsoid (WGS84 default)
        
    Returns:
        Tuple of (distance_in_meters, azimuth_deg, back_azimuth_deg)
        
    Notes:
        Uses Vincenty's formula for high precision on ellipsoid.
        Azimuth is from point 1 to point 2 (0° = North, 90° = East).
        Back-azimuth is from point 2 to point 1.
    """

def calc_vincenty_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
                         a: float = 6378137.0, f: float = 1/298.257223563,
                         max_iter: int = 200, tol: float = 1e-12) -> dict:
    """
    Vincenty's inverse formula for ellipsoid distance calculation.
    
    Args:
        lat1: First point latitude in degrees
        lon1: First point longitude in degrees
        lat2: Second point latitude in degrees  
        lon2: Second point longitude in degrees
        a: Semi-major axis in meters
        f: Flattening parameter
        max_iter: Maximum iterations for convergence
        tol: Convergence tolerance
        
    Returns:
        Dictionary with keys:
        - 'distance': Distance in meters
        - 'azimuth': Forward azimuth in degrees
        - 'reverse_azimuth': Reverse azimuth in degrees  
        - 'iterations': Number of iterations used
    """

def locations2degrees(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    """  
    Calculate angular distance between two points on sphere.
    
    Args:
        lat1: First point latitude in degrees
        lon1: First point longitude in degrees
        lat2: Second point latitude in degrees
        lon2: Second point longitude in degrees
        
    Returns:
        Angular distance in degrees
        
    Notes:
        Uses great circle distance on unit sphere.
        Suitable for approximate seismological calculations.
    """

Unit Conversions

Conversions between different distance units commonly used in seismology.

def degrees2kilometers(degrees: float, radius: float = 6371.0) -> float:
    """
    Convert degrees to kilometers using spherical Earth.
    
    Args:
        degrees: Angular distance in degrees
        radius: Earth radius in kilometers (mean radius default)
        
    Returns:
        Distance in kilometers
    """

def kilometers2degrees(kilometers: float, radius: float = 6371.0) -> float:
    """
    Convert kilometers to degrees using spherical Earth.
    
    Args:
        kilometers: Distance in kilometers
        radius: Earth radius in kilometers
        
    Returns:
        Angular distance in degrees
    """

def kilometer2degrees(kilometers: float, radius: float = 6371.0) -> float:
    """
    Convert kilometers to degrees (alias for backwards compatibility).
    
    Args:
        kilometers: Distance in kilometers
        radius: Earth radius in kilometers
        
    Returns:
        Angular distance in degrees
    """

Geographic Boundary Checking

Functions for determining if coordinates fall within specified geographic regions.

def inside_geobounds(longitude: float, latitude: float, 
                    min_longitude: float, max_longitude: float,
                    min_latitude: float, max_latitude: float) -> bool:
    """
    Check if point lies within geographic boundaries.
    
    Args:
        longitude: Point longitude in degrees
        latitude: Point latitude in degrees
        min_longitude: Western boundary in degrees
        max_longitude: Eastern boundary in degrees  
        min_latitude: Southern boundary in degrees
        max_latitude: Northern boundary in degrees
        
    Returns:
        True if point is inside boundaries
        
    Notes:
        Handles longitude wrapping around ±180°.
        Accounts for regions crossing the International Date Line.
    """

Flinn-Engdahl Regionalization

Seismological regionalization system for identifying geographic regions by coordinates.

class FlinnEngdahl:
    """
    Flinn-Engdahl seismological regionalization.
    
    Provides standardized geographic region identification used in
    earthquake catalogs and seismological studies.
    """
    
    def __init__(self):
        """Initialize Flinn-Engdahl regionalization system."""
    
    def get_region(self, longitude: float, latitude: float) -> str:
        """
        Get region name for geographic coordinates.
        
        Args:
            longitude: Longitude in degrees (-180 to 180)
            latitude: Latitude in degrees (-90 to 90)
            
        Returns:
            Region name string (e.g., "Southern California", "Japan")
        """
    
    def get_region_number(self, longitude: float, latitude: float) -> int:
        """
        Get numeric region code for coordinates.
        
        Args:
            longitude: Longitude in degrees
            latitude: Latitude in degrees
            
        Returns:
            Flinn-Engdahl region number (1-757)
        """

Usage Examples

Basic Distance and Azimuth Calculations

from obspy.geodetics import gps2dist_azimuth, locations2degrees

# Calculate distance between two seismic stations
sta1_lat, sta1_lon = 34.0, -118.0  # Los Angeles
sta2_lat, sta2_lon = 37.8, -122.4   # San Francisco

# High-precision ellipsoid calculation  
dist_m, az, baz = gps2dist_azimuth(sta1_lat, sta1_lon, sta2_lat, sta2_lon)
print(f"Distance: {dist_m/1000:.2f} km")
print(f"Azimuth LA→SF: {az:.1f}°")
print(f"Back-azimuth SF→LA: {baz:.1f}°")

# Quick spherical approximation
dist_deg = locations2degrees(sta1_lat, sta1_lon, sta2_lat, sta2_lon)  
print(f"Angular distance: {dist_deg:.3f}°")

Event-Station Distance Analysis

from obspy.geodetics import gps2dist_azimuth, degrees2kilometers
from obspy import UTCDateTime

# Earthquake and station information
event_lat, event_lon = 35.0, 140.0  # Japan earthquake
event_time = UTCDateTime("2023-01-15T10:30:00")

stations = {
    "Tokyo": (35.7, 139.7),
    "Seoul": (37.6, 126.9), 
    "Beijing": (39.9, 116.4),
    "Taipei": (25.0, 121.5),
    "Manila": (14.6, 121.0)
}

print("Station distances and azimuths from earthquake:")
print("Station".ljust(10), "Distance (km)".ljust(12), "Azimuth (°)".ljust(12), "Category")
print("-" * 50)

for station, (sta_lat, sta_lon) in stations.items():
    dist_m, az, baz = gps2dist_azimuth(event_lat, event_lon, sta_lat, sta_lon)
    dist_km = dist_m / 1000.0
    
    # Classify distance
    if dist_km < 100:
        category = "Local"
    elif dist_km < 1000:
        category = "Regional"
    else:
        category = "Teleseismic"
    
    print(f"{station:<10} {dist_km:>8.1f}    {az:>8.1f}      {category}")

Array Geometry Analysis

from obspy.geodetics import gps2dist_azimuth
import numpy as np

# Seismic array station coordinates
array_center = (40.0, -105.0)  # Colorado
stations = {
    "A01": (40.00, -105.00),  # Center station
    "A02": (40.01, -105.00),  # North
    "A03": (40.00, -105.01),  # West  
    "A04": (39.99, -105.00),  # South
    "A05": (40.00, -104.99),  # East
}

print("Array geometry analysis:")
print("Station  Distance (m)  Azimuth (°)")
print("-" * 35)

distances = []
azimuths = []

for station, (lat, lon) in stations.items():
    if station == "A01":  # Skip center station
        continue
        
    dist_m, az, baz = gps2dist_azimuth(array_center[0], array_center[1], lat, lon)
    distances.append(dist_m)
    azimuths.append(az)
    
    print(f"{station}     {dist_m:>8.1f}     {az:>7.1f}")

# Array statistics
print(f"\nArray statistics:")
print(f"Mean aperture: {np.mean(distances):.1f} m")
print(f"Aperture range: {np.min(distances):.1f} - {np.max(distances):.1f} m") 
print(f"Azimuth coverage: {np.min(azimuths):.1f}° - {np.max(azimuths):.1f}°")

Unit Conversion Utilities

from obspy.geodetics import degrees2kilometers, kilometers2degrees

# Convert between angular and linear distance units
angular_distances = [1.0, 5.0, 10.0, 30.0, 90.0, 180.0]  # degrees

print("Distance conversions:")
print("Degrees    Kilometers")
print("-" * 22)

for deg in angular_distances:
    km = degrees2kilometers(deg)
    print(f"{deg:>7.1f}    {km:>10.1f}")

print("\nReverse conversions:")
print("Kilometers    Degrees")
print("-" * 22)

linear_distances = [100, 500, 1000, 5000, 10000]  # km
for km in linear_distances:
    deg = kilometers2degrees(km)
    print(f"{km:>10.0f}    {deg:>7.3f}")

# Earth radius variations
print("\nDistance conversions with different Earth radii:")
radii = {"Mean": 6371.0, "Equatorial": 6378.1, "Polar": 6356.8}

for name, radius in radii.items():
    km = degrees2kilometers(10.0, radius=radius)
    print(f"10° using {name} radius ({radius} km): {km:.1f} km")

Geographic Region Classification

from obspy.geodetics import FlinnEngdahl

# Initialize Flinn-Engdahl regionalization
fe = FlinnEngdahl()

# Classify earthquake locations
earthquakes = [
    ("M6.5 California", 34.0, -118.0),
    ("M7.2 Japan", 35.0, 140.0),
    ("M8.0 Chile", -30.0, -71.0),
    ("M6.8 Turkey", 39.0, 35.0),
    ("M7.5 Alaska", 61.0, -150.0),
    ("M6.9 Indonesia", -6.0, 106.0)
]

print("Earthquake regional classification:")
print("Event".ljust(20), "Coordinates".ljust(15), "Flinn-Engdahl Region")
print("-" * 70)

for event_name, lat, lon in earthquakes:
    region = fe.get_region(lon, lat)  # Note: longitude first
    region_num = fe.get_region_number(lon, lat)
    coords = f"({lat:.1f}, {lon:.1f})"
    
    print(f"{event_name:<20} {coords:<15} {region} ({region_num})")

Station Network Spacing Analysis

from obspy.geodetics import gps2dist_azimuth
import numpy as np

# Global seismographic network stations (subset)
gsn_stations = {
    "ANMO": (34.9459, -106.4572),  # New Mexico, USA
    "ANTO": (39.8682, 32.7934),    # Turkey
    "ASAR": (23.2729, 5.6309),     # Algeria  
    "ASCN": (-7.9327, -14.3601),   # Ascension Island
    "BGCA": (52.3858, -113.3072),  # Canada
    "BOCO": (4.5869, -74.0432),    # Colombia
    "COCO": (-12.1901, 96.8349),   # Cocos Islands
    "CTAO": (-20.0882, 146.2545),  # Australia
    "FUNA": (-8.5259, 179.1966),   # Tuvalu
    "GUMO": (13.5893, 144.8684)    # Guam
}

print("Global network station spacing analysis:")
print("-" * 50)

# Calculate minimum distances between stations
min_distances = {}
for sta1, (lat1, lon1) in gsn_stations.items():
    min_dist = float('inf')
    closest_station = ""
    
    for sta2, (lat2, lon2) in gsn_stations.items():
        if sta1 == sta2:
            continue
            
        dist_m, _, _ = gps2dist_azimuth(lat1, lon1, lat2, lon2)
        if dist_m < min_dist:
            min_dist = dist_m
            closest_station = sta2
    
    min_distances[sta1] = (min_dist/1000.0, closest_station)

# Sort by distance
sorted_stations = sorted(min_distances.items(), key=lambda x: x[1][0])

print("Station    Closest Station    Distance (km)")
print("-" * 45)
for station, (dist, closest) in sorted_stations:
    print(f"{station:<10} {closest:<15}    {dist:>8.1f}")

# Network statistics
distances = [dist for dist, _ in min_distances.values()]
print(f"\nNetwork spacing statistics:")
print(f"Mean minimum distance: {np.mean(distances):.1f} km")
print(f"Median minimum distance: {np.median(distances):.1f} km")
print(f"Distance range: {np.min(distances):.1f} - {np.max(distances):.1f} km")

Boundary and Region Checking

from obspy.geodetics import inside_geobounds

# Define study regions
regions = {
    "California": {"min_lon": -125.0, "max_lon": -114.0, 
                  "min_lat": 32.0, "max_lat": 42.0},
    "Japan": {"min_lon": 129.0, "max_lon": 146.0,
             "min_lat": 30.0, "max_lat": 46.0}, 
    "Mediterranean": {"min_lon": -10.0, "max_lon": 45.0,
                     "min_lat": 30.0, "max_lat": 48.0}
}

# Test earthquake locations
test_events = [
    ("Northridge", 34.2, -118.5),
    ("Kobe", 34.6, 135.0),
    ("L'Aquila", 42.3, 13.4),
    ("Christchurch", -43.5, 172.6),
    ("Mexico City", 19.4, -99.1)
]

print("Event location classification:")
print("Event".ljust(15), "Coordinates".ljust(18), "Regions")
print("-" * 50)

for event_name, lat, lon in test_events:
    coords = f"({lat:.1f}, {lon:.1f})"
    regions_found = []
    
    for region_name, bounds in regions.items():
        if inside_geobounds(lon, lat, 
                           bounds["min_lon"], bounds["max_lon"],
                           bounds["min_lat"], bounds["max_lat"]):
            regions_found.append(region_name)
    
    regions_str = ", ".join(regions_found) if regions_found else "None"
    print(f"{event_name:<15} {coords:<18} {regions_str}")

Types

# Coordinate point structure
Point = {
    'latitude': float,        # Latitude in degrees
    'longitude': float,       # Longitude in degrees
    'elevation': float        # Elevation in meters (optional)
}

# Distance calculation result
DistanceResult = {
    'distance_m': float,      # Distance in meters
    'distance_km': float,     # Distance in kilometers  
    'distance_deg': float,    # Angular distance in degrees
    'azimuth': float,         # Forward azimuth in degrees
    'back_azimuth': float     # Back azimuth in degrees
}

# Geographic boundary definition
GeoBounds = {
    'min_latitude': float,    # Southern boundary
    'max_latitude': float,    # Northern boundary
    'min_longitude': float,   # Western boundary
    'max_longitude': float    # Eastern boundary
}

# Flinn-Engdahl region information
RegionInfo = {
    'number': int,            # Region number (1-757)
    'name': str,              # Region name
    'type': str               # Region type (geographic/seismic)
}

Install with Tessl CLI