plate-tectonics-geospatial

name	plate-tectonics-geospatial
description	Analyze plate tectonics data using GeoPandas, identify points within plates, and calculate distances to boundaries.

Plate Tectonics Geospatial Analysis

Overview

The PB2002 (Plate Boundaries 2002) dataset provides comprehensive data on tectonic plate boundaries and plate definitions. This skill covers loading, processing, and spatial analysis of plate data.

PB2002 Dataset Structure

Boundary Data (PB2002_boundaries.json)

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "LineString",
        "coordinates": [[lon1, lat1], [lon2, lat2], ...]
      },
      "properties": {
        "PLATE1": "Pacific",
        "PLATE2": "North American",
        "TYPE": "subduction zone",  // or "transform", "spreading", etc.
        "STEPOVER": ""
      }
    }
  ]
}

Plate Data (PB2002_plates.json)

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "MultiPolygon",  // or Polygon
        "coordinates": [...]
      },
      "properties": {
        "PlateName": "Pacific",
        "PlateID": 101
      }
    }
  ]
}

Loading and Processing

Load Plate Boundaries

import geopandas as gpd
import json

with open('/root/PB2002_boundaries.json', 'r') as f:
    boundaries_geojson = json.load(f)

boundaries_gdf = gpd.GeoDataFrame.from_features(
    boundaries_geojson['features'],
    crs='EPSG:4326'
)

# Filter for specific plate boundary
pacific_boundaries = boundaries_gdf[
    (boundaries_gdf['PLATE1'] == 'Pacific') |
    (boundaries_gdf['PLATE2'] == 'Pacific')
]

Load Plate Polygons

with open('/root/PB2002_plates.json', 'r') as f:
    plates_geojson = json.load(f)

plates_gdf = gpd.GeoDataFrame.from_features(
    plates_geojson['features'],
    crs='EPSG:4326'
)

# Get Pacific plate polygon
pacific_plate = plates_gdf[plates_gdf['PlateName'] == 'Pacific']

Spatial Operations

Identify Points Within Plate

# Use spatial join to find earthquakes within Pacific plate
earthquakes_in_pacific = gpd.sjoin(
    earthquakes_gdf,
    pacific_plate,
    how='inner',
    predicate='within'
)

Calculate Distance to Boundary

# For each earthquake, calculate minimum distance to boundary
def min_distance_to_boundary(point, boundary_lines_gdf):
    """Calculate minimum distance from point to any boundary line"""
    min_dist = float('inf')
    for idx, boundary in boundary_lines_gdf.iterrows():
        dist = point.distance(boundary['geometry'])
        if dist < min_dist:
            min_dist = dist
    return min_dist

# Apply to all earthquakes in the plate
# First, project to projected CRS for accurate distance in km
earthquakes_projected = earthquakes_in_pacific.to_crs('EPSG:3857')
boundaries_projected = pacific_boundaries.to_crs('EPSG:3857')

earthquakes_projected['distance_to_boundary'] = earthquakes_projected.geometry.apply(
    lambda point: min_distance_to_boundary(point, boundaries_projected)
)

# Convert from meters to kilometers
earthquakes_projected['distance_km'] = earthquakes_projected['distance_to_boundary'] / 1000

Find Maximum Distance

# Find earthquake furthest from boundary
furthest_idx = earthquakes_projected['distance_km'].idxmax()
furthest_earthquake = earthquakes_projected.loc[furthest_idx]

Important Considerations

Polygon Validation

# Ensure polygons are valid before spatial operations
if not plates_gdf.geometry.is_valid.all():
    plates_gdf['geometry'] = plates_gdf.geometry.buffer(0)

if not boundaries_gdf.geometry.is_valid.all():
    boundaries_gdf['geometry'] = boundaries_gdf.geometry.buffer(0)

Handling MultiPolygons

# If a plate is represented as MultiPolygon, create a union
pacific_plate_union = pacific_plate.unary_union

Projection for Accurate Distances

# Use Web Mercator (EPSG:3857) for global distances in meters
# Then convert to kilometers

# For regional analysis, consider UTM zones
# EPSG:32633 = UTM Zone 33N
# EPSG:32733 = UTM Zone 33S

Common Pitfalls

1. Not projecting before distance calculations - Distances in geographic CRS are meaningless
2. Invalid geometries - Use .buffer(0) to fix self-intersecting polygons
3. Antimeridian issues - Points near ±180° longitude may have wrapping issues
4. MultiPolygon confusion - Remember that a plate may consist of multiple disconnected polygons