Geospatial earthquake analysis

Geospatial Earthquake and Plate Boundary Analysis

Overview

This skill provides a complete workflow for analyzing earthquake locations relative to tectonic plate boundaries using GeoPandas. It handles data loading, coordinate projection management, and spatial distance calculations.

Key Concepts

Coordinate Projections

• Earthquake data typically comes in WGS84 (EPSG:4326) - latitude/longitude
• Distance calculations require a projected coordinate system (meters/kilometers)
• Use appropriate azimuthal equidistant or equal-area projections for accurate distance measurements
• Project data before calculating distances; transform results back if needed

Plate Boundary Data Structure

• Boundaries are typically represented as LineStrings (boundary lines)
• Plates are regions bounded by these lines
• A point is "within a plate" if it's inside the polygon formed by the plate's boundaries
• "Distance to boundary" = minimum distance from a point to any boundary line of that plate

Workflow

1. Loading Earthquake Data

import geopandas as gpd
import json
from shapely.geometry import Point

# Load earthquake JSON
with open('/root/earthquakes_2024.json', 'r') as f:
    earthquakes = json.load(f)

# Convert to GeoDataFrame
features = []
for eq in earthquakes['features']:
    props = eq['properties']
    geom = Point(eq['geometry']['coordinates'])
    features.append({'geometry': geom, 'id': props.get('id'), 'magnitude': props.get('mag'), ...})

earthquakes_gdf = gpd.GeoDataFrame(features, crs='EPSG:4326')

2. Loading Plate Boundary Data

# Load boundaries as LineStrings
boundaries_gdf = gpd.read_file('/root/PB2002_boundaries.json')  # Has CRS info
boundaries_gdf = boundaries_gdf.to_crs('EPSG:4326')  # Ensure consistent CRS

# Load plate definitions
plates_gdf = gpd.read_file('/root/PB2002_plates.json')
plates_gdf = plates_gdf.to_crs('EPSG:4326')

3. Identifying Earthquakes Within a Specific Plate

# Method: Use spatial join to find earthquakes within plate polygons
earthquakes_in_plate = gpd.sjoin(
    earthquakes_gdf,
    plates_gdf[plates_gdf['plate_name'] == 'Pacific'],
    how='inner',
    predicate='within'
)

4. Calculating Distance to Plate Boundaries

For each earthquake within the Pacific plate, calculate its minimum distance to the Pacific plate boundary:

# Project to a suitable projection for distance (meters)
# Azimuthal equidistant projection centered near the Pacific
pacific_proj = 'EPSG:3832'  # Or use a custom projection

earthquakes_projected = earthquakes_gdf.to_crs(pacific_proj)
boundaries_projected = boundaries_gdf.to_crs(pacific_proj)

# Calculate distance from each earthquake point to the nearest boundary
def distance_to_boundary(point, boundary_lines):
    min_dist = float('inf')
    for line in boundary_lines.geometry:
        dist = point.distance(line)
        min_dist = min(min_dist, dist)
    return min_dist / 1000  # Convert to km

distances = []
for idx, eq in earthquakes_in_plate.iterrows():
    dist = distance_to_boundary(eq.geometry, boundaries_projected)
    distances.append(dist)

earthquakes_in_plate['distance_to_boundary_km'] = distances

5. Finding the Furthest Earthquake

furthest = earthquakes_in_plate.loc[earthquakes_in_plate['distance_to_boundary_km'].idxmax()]
print(f"Furthest from boundary: {furthest['place']}, {furthest['distance_to_boundary_km']:.2f} km")

Important Notes

• CRS Consistency: Always ensure all layers use the same CRS before spatial operations
• Projection Selection: Use projections appropriate for your region (azimuthal equidistant minimizes distortion)
• Boundary Definition: The Pacific plate boundary includes all associated trenches and ridges
• Distance Metric: Minimum distance accounts for complex geometry; straight-line distance is not appropriate
• Performance: For large datasets, use spatial indexing and consider chunking operations

Output Format

Return earthquake details with distance information:

result = {
    'id': earthquake_id,
    'place': earthquake_location,
    'time': earthquake_time_iso8601,
    'magnitude': earthquake_mag,
    'latitude': eq_lat,
    'longitude': eq_lon,
    'distance_km': distance_rounded_2dp
}