Plate Tectonics Geospatial

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)

jsonCopy code
{
  "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)

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

Loading and Processing

Load Plate Boundaries

pythonCopy code
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

pythonCopy code
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

pythonCopy code
# 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

pythonCopy code
# 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

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

Important Considerations

Polygon Validation

pythonCopy code
# 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

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

Projection for Accurate Distances

pythonCopy code
# 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