Plate boundary data processing

Plate Boundary Data Processing

PB2002 Dataset Overview

The PB2002 (Plate Boundaries 2002) dataset contains:

• Boundaries file: LineStrings representing plate boundaries (ridges, trenches, transforms)
• Plates file: Polygon features representing individual plate boundaries and characteristics

Data Structure

PB2002_boundaries.json

• Contains LineString geometries for each plate boundary segment
• Properties typically include:
  • type: Boundary type (ridge, trench, transform, etc.)
  • plate1, plate2: Names of adjacent plates
  • age, length: Metadata about boundary

PB2002_plates.json

• Contains Polygon geometries for complete plates
• Properties typically include:
  • plate_name: Name of the plate (e.g., "Pacific")
  • area: Plate area
  • Other identifying characteristics

Loading and Exploring Data

Basic Loading

import geopandas as gpd
import json

# Load plate boundaries
boundaries_gdf = gpd.read_file('/root/PB2002_boundaries.json')
print(f"Loaded {len(boundaries_gdf)} boundary segments")
print(boundaries_gdf.head())
print(boundaries_gdf.columns)
print(boundaries_gdf.crs)

# Load plate definitions
plates_gdf = gpd.read_file('/root/PB2002_plates.json')
print(f"Loaded {len(plates_gdf)} plates")
print(plates_gdf['plate_name'].unique())

Understanding Data Types

# Check geometry types
print(boundaries_gdf.geometry.type.value_counts())
# Should mostly show LineString

print(plates_gdf.geometry.type.value_counts())
# Should show Polygon

Identifying the Pacific Plate

Method 1: Extract from Plates GeoDataFrame

# Find Pacific plate polygon
pacific_plate = plates_gdf[plates_gdf['plate_name'] == 'Pacific']

if len(pacific_plate) > 0:
    print(f"Found Pacific plate")
    print(f"Area: {pacific_plate.iloc[0].geometry.area}")
    pacific_boundary_polygon = pacific_plate.iloc[0].geometry
else:
    print("Pacific plate not found - check plate_name values")
    print(plates_gdf['plate_name'].unique())

Method 2: Extract Boundaries Belonging to Pacific

# If using plate1/plate2 columns, filter boundaries of the Pacific plate
pacific_boundaries = boundaries_gdf[
    (boundaries_gdf['plate1'] == 'Pacific') |
    (boundaries_gdf['plate2'] == 'Pacific')
]

print(f"Found {len(pacific_boundaries)} Pacific plate boundary segments")

Filtering Earthquakes Within a Plate

Using Spatial Join

# Earthquake GeoDataFrame already loaded in WGS84
earthquakes_in_pacific = gpd.sjoin(
    earthquakes_gdf,
    pacific_plate,
    how='inner',
    predicate='within'
)

print(f"Found {len(earthquakes_in_pacific)} earthquakes within Pacific plate")

Manual Point-in-Polygon

# Alternative if spatial join has issues
pacific_polygon = pacific_plate.iloc[0].geometry

earthquakes_in_pacific = earthquakes_gdf[
    earthquakes_gdf.geometry.within(pacific_polygon)
]

Extracting Boundaries for Distance Calculation

Get All Pacific Plate Boundary Lines

# Collect all boundary LineStrings for the Pacific plate
pacific_boundary_lines = boundaries_gdf[
    (boundaries_gdf['plate1'] == 'Pacific') |
    (boundaries_gdf['plate2'] == 'Pacific')
]['geometry'].tolist()

print(f"Pacific plate has {len(pacific_boundary_lines)} boundary segments")

Important Consideration: Boundary vs Plate Polygon

• Boundaries GeoDataFrame: Individual line segments that make up the plate edge
• Plates GeoDataFrame: Complete polygon representing the full plate area
• For distance calculation, use boundary lines (not the polygon boundary)
• The polygon exists but distance should be to actual plate boundary data

Handling Data Issues

Case Sensitivity

Plate names may vary in case:

# Make case-insensitive if needed
plates_gdf['plate_name_lower'] = plates_gdf['plate_name'].str.lower()
pacific = plates_gdf[plates_gdf['plate_name_lower'] == 'pacific']

Missing or Inconsistent Naming

# Explore actual values
print(plates_gdf.columns.tolist())
print(plates_gdf['plate_name'].value_counts())

# Check if 'Pacific' exists in plate1 or plate2 columns
if 'plate1' in boundaries_gdf.columns:
    print(boundaries_gdf['plate1'].unique())
    print(boundaries_gdf['plate2'].unique())

CRS Harmonization

# Ensure consistent CRS before spatial operations
if boundaries_gdf.crs != plates_gdf.crs:
    boundaries_gdf = boundaries_gdf.to_crs(plates_gdf.crs)

if earthquakes_gdf.crs != plates_gdf.crs:
    earthquakes_gdf = earthquakes_gdf.to_crs(plates_gdf.crs)

Common Workflow

# 1. Load all data
boundaries_gdf = gpd.read_file('/root/PB2002_boundaries.json')
plates_gdf = gpd.read_file('/root/PB2002_plates.json')
earthquakes_gdf = load_earthquake_data()  # Your function

# 2. Harmonize CRS
for gdf in [boundaries_gdf, plates_gdf, earthquakes_gdf]:
    if gdf.crs != 'EPSG:4326':
        gdf = gdf.to_crs('EPSG:4326')

# 3. Identify Pacific plate
pacific_plate = plates_gdf[plates_gdf['plate_name'] == 'Pacific']
pacific_boundaries = boundaries_gdf[
    (boundaries_gdf['plate1'] == 'Pacific') |
    (boundaries_gdf['plate2'] == 'Pacific')
]

# 4. Filter earthquakes
earthquakes_pacific = gpd.sjoin(
    earthquakes_gdf,
    pacific_plate,
    how='inner',
    predicate='within'
)

# 5. Calculate distances and find maximum
# (See geospatial-earthquake-analysis skill)

Verification Checklist

• ✓ Loaded boundaries and plates with same CRS
• ✓ Located Pacific plate in plates_gdf
• ✓ Identified Pacific boundaries in boundaries_gdf
• ✓ Found earthquakes within Pacific polygon
• ✓ Boundary lines are LineStrings, not just points
• ✓ Earthquake count reasonable (not 0 or unreasonably high)