/dem-pipeline - DEM Lineament Analysis Skill

Purpose

Run Digital Elevation Model (DEM) analysis pipeline to identify unmapped fault structures through edge detection, drainage pattern analysis, and topographic lineament extraction. Validated in Italy where DEM revealed 25+ km unmapped fault.

Usage

/dem-pipeline <region> [--dem-path PATH] [--output-dir DIR] [--overlay GEOJSON]

Examples:

/dem-pipeline italy --dem-path dem_tiles/tinitaly_basura.tif
/dem-pipeline san-diego --overlay data/fault_databases/scec_cfm_sandiego.geojson
/dem-pipeline turkey --dem-path dem_tiles/sofular_srtm.tif

Prerequisites

Required Data

Python Dependencies

numpy, scipy, matplotlib, rasterio, geopandas,
scikit-image (edge detection), pysheds (drainage)

Pipeline Steps

Step 1: Load and Preprocess DEM

# Load DEM
dem = rasterio.open(dem_path)
elevation = dem.read(1)

# Resample if needed (target: 10-30m resolution)
if resolution > 30:
    elevation = resample_dem(elevation, target_res=30)

# Fill voids/nodata
elevation = fill_voids(elevation)

Step 2: Generate Hillshade

# Calculate hillshade for visualization
hillshade = calculate_hillshade(
    elevation,
    azimuth=315,  # NW illumination
    altitude=45
)

# Save as GeoTIFF
save_geotiff(hillshade, 'hillshade.tif')

Step 3: Edge Detection

Apply multiple edge detection algorithms:

# Sobel edge detection (gradient-based)
sobel_edges = sobel(elevation)

# Canny edge detection (multi-scale)
canny_edges = canny(elevation, sigma=2)

# Save results
save_geotiff(sobel_edges, 'edge_sobel.tif')
save_geotiff(canny_edges, 'edge_canny.tif')

Look for:

• Linear features at consistent azimuths
• Edges that align with orphan earthquake clusters
• Topographic breaks/scarps

Step 4: Drainage Network Analysis

# Calculate flow direction
flow_dir = calculate_flow_direction(elevation)

# Calculate flow accumulation
flow_acc = calculate_flow_accumulation(flow_dir)

# Extract drainage network (threshold: 1000 cells)
drainage = extract_drainage(flow_acc, threshold=1000)

# Save results
save_geotiff(drainage, 'drainage.tif')

Look for:

• Drainage anomalies (right-angle bends, beheaded streams)
• Linear valley alignments
• Offset drainage patterns (strike-slip indicators)

Step 5: Lineament Extraction

# Automated lineament extraction
lineaments = extract_lineaments(
    sobel_edges,
    min_length=5000,  # 5 km minimum
    max_gap=500       # 500m max gap
)

# Calculate azimuths
for lineament in lineaments:
    lineament.azimuth = calculate_azimuth(lineament)
    lineament.length = calculate_length(lineament)

# Group by azimuth
azimuth_groups = group_by_azimuth(lineaments, bin_size=15)

Step 6: Overlay Analysis

# Load fault database
faults = gpd.read_file(fault_geojson)

# Load microseismicity
earthquakes = pd.read_csv(eq_csv)
orphans = earthquakes[earthquakes['is_orphan'] == True]

# Create overlay figure
fig, ax = plt.subplots(figsize=(12, 10))
ax.imshow(hillshade, cmap='gray')
faults.plot(ax=ax, color='red', linewidth=2, label='Mapped faults')
ax.scatter(orphans.lon, orphans.lat, c='yellow', s=10, label='Orphan EQs')
lineaments.plot(ax=ax, color='cyan', linewidth=1, label='DEM lineaments')

plt.savefig('overlay_analysis.png', dpi=300)

Output Products

1. GeoTIFF Rasters

hillshade.tif

edge_sobel.tif

edge_canny.tif

drainage.tif

slope.tif

aspect.tif

2. Publication Figures (300 DPI PNG)

fig_hillshade.png

fig_edge_detection.png

fig_drainage.png

fig_overlay.png

fig_azimuth_rose.png

3. Analysis Report

## DEM Lineament Analysis: [Region]

**DEM source**: [TINITALY/SRTM/etc.]
**Resolution**: [X] m
**Analysis date**: [date]

---

### Lineament Statistics

| Azimuth Bin | Count | Total Length | Avg Length |
|-------------|-------|--------------|------------|
| N-S (0-15°) | [N] | [X] km | [X] km |
| NNE (15-30°) | [N] | [X] km | [X] km |
| ... | ... | ... | ... |

### Dominant Lineament Sets

1. **[Azimuth]° trend**: [N] lineaments, [X] km total
   - Correlation with: [orphan cluster / mapped fault / etc.]

2. **[Azimuth]° trend**: ...

---

### Comparison with Mapped Faults

| Mapped Fault | Strike | DEM Lineament Match? |
|--------------|--------|---------------------|
| [fault name] | [X]° | [YES/NO] |

### Unmapped Structures Identified

| Structure | Azimuth | Length | Evidence |
|-----------|---------|--------|----------|
| NNW lineament | 340° | 25+ km | Edge + drainage + 85 orphan EQs |

---

### Interpretation

[2-3 paragraphs on geological interpretation]

---

### Files Created

| File | Path |
|------|------|
| Hillshade | `dem_tiles/[region]/hillshade.tif` |
| Edge detection | `dem_tiles/[region]/edge_sobel.tif` |
| Overlay figure | `dem_tiles/[region]/fig_overlay.png` |

Example: Italy Bàsura Results

DEM Lineament Analysis: Ligurian Alps (Italy)

DEM: TINITALY 10m
Area: 50 km radius around Bàsura Cave

Key Finding: NNW-SSE Lineament Set

- Azimuth: 335-345° (NNW)
- Length: 25+ km
- Evidence:
  - Strong edge detection signal
  - Anomalous drainage (right-angle bends)
  - 85 orphan earthquakes aligned
  - NOT in DISS, ITHACA, or EFSM20

Interpretation:
Previously unmapped fault structure, primary candidate
for 1285/1394 dark earthquakes. Strike consistent with
regional stress field. Recommend InSAR/GPS validation.

Integration with Other Skills

/orphan-analysis

/verify-dark

/regional-doc

DEM_LINEAMENT_FINDINGS.md

Scripts Reference

Existing scripts in

paleoseismic_caves/scripts/

dem_lineament_analysis.py

san_diego_dem_processing.py

san_diego_lineament_analysis.py

san_diego_integrated_overlay.py

Important Notes

1. Resolution matters - 10-30m optimal for fault detection
2. Multiple methods - Sobel + Canny + drainage gives confidence
3. Ground truth needed - DEM alone is not proof
4. Orphan correlation key - Lineaments matching orphan clusters are strongest evidence
5. Publication quality - 300 DPI, proper color scales, scale bars

Regional DEM Sources