GeoAI Skills

View saved version →

---
name: arcgis-pro-automation
description: >-
  Automate controlled local ArcGIS Pro and ArcPy workflows through
  arcgis-mcp-bridge: inspect .aprx projects and file geodatabases, run
  geoprocessing, projection, raster, network, spatial-statistics, editing,
  symbology, and layout export with path and mutation guards. Use when the
  user explicitly names ArcGIS Pro, ArcPy, .aprx, .gdb, Esri geoprocessing,
  muend/arcgis-mcp-bridge, its health_check, PathGuard, or confirmation
  gates, or sketch-to-GIS extraction. Do not trigger for ArcGIS Online or
  Enterprise administration, QGIS or PyQGIS, generic open-source GIS, or
  live GUI control of an already-open ArcGIS Pro session.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# ArcGIS Pro Automation

Operate saved ArcGIS Pro projects and local GIS datasets through the guarded
`arcgis-mcp-bridge` tool surface. Treat the bridge as an execution backend,
not as permission to mutate data or claim a result that was never observed.

## Establish the boundary

- Require Windows, a licensed ArcGIS Pro installation, and a bridge worker
  interpreter that can import ArcPy for real geoprocessing.
- Use the bridge for local, headless, repeatable ArcPy work. It does not drive
  the visible ArcGIS Pro UI and does not administer ArcGIS Online or Enterprise.
- If the bridge tools are not callable, provide setup or a dry execution plan.
  Never simulate tool output or imply that an `.aprx` or `.gdb` changed.
- Use the schemas exposed by the connected MCP server. Tool names may be host-
  namespaced; match the semantic catalog name and never invent parameters.
- Keep every read and write inside the configured allowed roots. Use absolute
  paths and a dedicated scratch geodatabase for intermediate outputs.

Read [runtime and licensing](references/runtime-and-licensing.md) before the
first execution or whenever setup, extension availability, or worker failures
are in scope.

## Execute the workflow

1. **Define the contract.** Identify inputs, intended output artifacts, CRS and
   units, whether a saved `.aprx` must change, and which operations create,
   overwrite, append, or edit data.
2. **Preflight.** Call `health_check` first. Inspect its worker interpreter,
   allowed roots, scratch GDB, timeout, and concurrency. Because this call does
   not import ArcPy, follow it with a non-mutating ArcPy-backed request such as
   `get_spatial_reference`, `describe_dataset`, or `list_layers` before treating
   the runtime as execution-ready.
3. **Inspect before acting.** Read dataset descriptions, counts, fields, CRS,
   extents, project maps/layers, and required extension licenses. Resolve datum
   transformations, units, and output naming before analysis.
4. **Plan exact tools.** Select the smallest declarative sequence from
   [tool routing](references/tool-routing.md). Separate read-only inspection,
   new-output creation, and in-place mutation. Prefer new outputs and working
   copies over edits to source data.
5. **Authorize mutations.** Read [safety and validation](references/safety-and-validation.md)
   before any write. Set `confirm=true` only when the user's request clearly
   authorizes the exact mutating operation and target. Treat `overwrite=true`
   as a separate explicit decision. Stop on ambiguous target, scope, or intent.
6. **Execute incrementally.** Run one dependent step at a time; preserve returned
   output paths and geoprocessing messages. Parallelize only independent jobs
   when license seats, memory, and `max_workers` are known to support them.
7. **Verify independently.** Re-open outputs and test the relevant invariants:
   existence, count, geometry validity, CRS and units, raster statistics and
   NoData, network solve status, project contents, or exported layout dimensions.
   A success status alone is not evidence of a correct GIS result.

## Route domain method and execution separately

Use this skill for the ArcGIS execution layer. Pair it with the relevant domain
skill when method choice or scientific validity is substantive:

- DEM, slope, drainage, or viewsheds → `terrain-hydrology`
- routing, service areas, OD, or facilities → `network-accessibility-analysis`
- Moran's I, Gi*, kernel density, or spatial inference → `spatial-statistics`
- raster/imagery preprocessing and interpretation → `remote-sensing-analysis`
- map design before `.aprx` styling or layout export → `cartography-geoviz`
- multi-stage cross-domain delivery → `geoai-orchestrator`

Do not activate ArcGIS automation merely because an open format can be read by
ArcGIS. A GeoPackage, GeoJSON, raster, or `.gdb` request that explicitly chooses
GDAL, GeoPandas, QGIS, PostGIS, or another runtime belongs to that specialist.

## Handle failures without bypasses

- `validation`: correct the payload against the live schema; do not relax types.
- `security`: move or copy only with user authorization; never broaden allowed
  roots merely to make a call pass.
- `license`: report the unavailable base or extension license. Use an alternative
  only when it is methodologically equivalent and disclose the change.
- `geoprocessing`: preserve ArcPy messages, inspect inputs/environment, and retry
  only after a concrete correction. Never retry an in-place mutation blindly.
- `internal` or worker crash: preserve the error boundary and stop claiming state.

## Deliver evidence

Report the tool sequence, inspected inputs, exact output paths, mutations and
authorization basis, relevant ArcPy messages, verification results, license and
runtime constraints, and unresolved limitations. For current tool or platform
claims, consult [the authoritative source registry](references/authoritative-sources.md)
and record the checked date.

## Execution contract

- **Workflow:** establish the local ArcGIS boundary; preflight server and ArcPy separately; inspect data and project state; plan exact tools; authorize writes; execute incrementally; reopen and verify outputs.
- **Decision rules:** use the bridge only for explicit local ArcGIS Pro or ArcPy work, pair it with domain-method skills when needed, prefer new outputs, and reserve confirmation and overwrite flags for clearly authorized targets.
- **Verification protocol:** reconcile counts, geometry, CRS and units, raster or network diagnostics, `.aprx` contents, exported artifacts, returned paths, and geoprocessing messages against the stated success criteria.
- **Failure modes:** stop for absent bridge tools, failed ArcPy preflight, unknown CRS or datum transform, paths outside allowed roots, missing licenses, ambiguous mutation scope, worker crashes, or unverified output state.
- **Deliverables:** reproducible tool plan and calls, exact inputs and outputs, mutation record, ArcPy messages, validation evidence, runtime and license provenance, and limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before applying bridge, ArcPy, licensing, or platform rules and record the checked date.

Referenced files: 5

View saved version →

---
name: cartography-geoviz
description: >-
  Always invoke before answering any request to create, compare, design, or
  review a user-facing map, even if the request is terse or underspecified.
  Covers publication maps, choropleths, map series and small multiples,
  comparable multi-date panels, proportional/bivariate/flow maps, raster
  rendering, and interactive web maps. Includes classification, color,
  legends, projections, accessibility, and large-data aggregation. Do not
  trigger for a temporary diagnostic plot inside another analysis.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Cartography & Geovisualization

Purpose: maps that communicate honestly. Cartographic choices (class
breaks, ramps, normalization, projection) can manufacture or hide
patterns; this skill treats them as analytical decisions with stated
rationale, not styling.

## The first three questions

1. **What's the message?** One map = one message. If two variables
   compete, consider small multiples or a bivariate scheme — not twelve
   legend classes.
2. **Normalized?** Choropleths of raw counts are population maps in
   disguise. Rates, densities, or per-capita for area-based color; raw
   magnitudes → proportional symbols instead.
3. **Static or interactive?** Print/PDF/paper → matplotlib/QGIS layout;
   exploration/stakeholders → Folium/MapLibre; big point data →
   Kepler.gl/deck.gl (GPU).

## Thematic map type selection

| Data | Map type |
|---|---|
| Rate/ratio by polygon | Choropleth |
| Count/magnitude by place | Proportional/graduated symbols |
| Two related rates | Bivariate choropleth (3×3 max) |
| Individual-level density | Dot density or KDE surface (label bandwidth) |
| Continuous field (raster) | Classified or stretched render + hillshade context |
| Movement/OD | Flow map (width∝volume), aggregate to avoid hairballs |
| Change over time | Small multiples > animation for analysis; animation for outreach |

## Classification — the honesty lever

- **Natural breaks (Jenks)**: default for skewed data; breaks are
  data-specific, so NOT comparable across maps/dates.
- **Quantiles**: guaranteed color balance; can split near-identical values.
- **Equal interval**: comparable and intuitive; fails on skew.
- **Manual/defined**: the ONLY correct choice for map series (same breaks
  across all dates/regions) and for domain thresholds (WHO limits, slope
  classes).
- 5±2 classes; show the histogram with breaks in the workflow; state the
  scheme in the caption/metadata. Try two schemes — if the story changes
  materially, the story is the classification, and the reader must be told.

## Color

- Ramps from ColorBrewer/`cmcrameri`/viridis family: sequential (ordered),
  diverging (meaningful midpoint — zero, mean, threshold), qualitative
  (categories, ≤ 8).
- Colorblind-safe by default (~8% of male readers); never red-green
  diverging without checking a CVD simulator.
- NoData ≠ zero: render as neutral gray with its own legend entry, never
  the ramp's low end.
- Muted basemaps (CartoDB Positron) under thematic layers — the basemap
  must never win.

## Projection for display

- Web tiles = Web Mercator: fine for city scale; area comparisons at
  continental scale on Mercator are visual lies — use equal-area
  projections (Albers, Mollweide, Equal Earth) for static thematic maps of
  large extents.
- National mapping → the national grid; polar work → polar stereographic.
- Label the projection on publication maps.

## Required furniture (publication static maps)

Title (the message, not the filename), legend (units!, sensible number
formatting), scale bar (projected CRS only — degrees have no fixed scale),
north arrow (only when north isn't up or the audience expects it), data
source + date + projection + author, and an inset locator map for
unfamiliar regions.

```python
# GeoPandas static map core
ax = gdf.plot(column="rate_per_1k", scheme="naturalbreaks", k=5,
              cmap="YlGnBu", legend=True, edgecolor="white", linewidth=0.3,
              missing_kwds={"color": "#d9d9d9", "label": "No data"})
ax.set_axis_off()
```

Export: 300 dpi PNG/PDF for print; SVG when editors will touch it; COG +
style for GIS handoff.

## Interactive maps

- Folium/MapLibre: tooltips with formatted values, layer control, sensible
  initial bounds (`fit_bounds`), legend included (Folium needs a manual
  HTML/branca legend — don't ship without one).
- Performance: >~50k vector features → tile it (tippecanoe → PMTiles) or
  switch to deck.gl/Kepler; never dump 500k GeoJSON features into Leaflet.
- Every popup number formatted (thousands separators, units, rounding
  matched to precision honesty).

## Verification protocol

1. Squint test: does the message survive at thumbnail size?
2. CVD simulation pass.
3. Legend audit: units, rounding, class edges non-overlapping.
4. Cross-check 3 features' rendered values against the attribute table
   (classification bugs are silent).
5. For map series: identical breaks, ramp, and extent across panels.

## Pitfalls checklist

- Raw-count choropleth (population in disguise).
- Jenks breaks compared across two dates.
- Red-green diverging ramp, unlabeled midpoint.
- NoData painted as the lowest class.
- Scale bar on an unprojected (degree) map.
- Continental-area comparisons on Web Mercator.
- Interactive map with no legend or units.

## Execution contract

- **Workflow:** inspect audience, data semantics, scale, and output medium; select projection, normalization, classification, and visual hierarchy; render; verify; export.
- **Decision rules:** choose map type from the analytical question, normalize counts when exposure differs, and keep breaks fixed for comparisons.
- **Verification protocol:** run the five checks above and reconcile rendered values, units, class edges, and missing-data treatment against the source.
- **Failure modes:** stop or qualify delivery when denominators, CRS, units, accessibility, or cross-panel comparability are unresolved.
- **Deliverables:** final map, legend and units, data/source note, projection and classification rationale, accessibility note, and reproducible style or code.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before using version-sensitive APIs and record the checked date.

Referenced files: 2

View saved version →

---
name: change-detection
description: >-
  Change analysis, once the observations are comparable. Not for cases whose
  blocker is comparability itself: mixed sensors, product levels or processing
  baselines to remote-sensing-analysis, undocumented vertical datums to
  point-cloud-lidar, multi-decade archive trends over large areas to
  google-earth-engine. Matching product level does not prove comparability.
  Otherwise invoke for what, where or how much changed: two-scene comparison,
  deforestation, urban growth, disaster damage, parcel-change audits,
  bi-temporal differencing, post-classification comparison, adjusted area,
  break detection in a series in hand (BFAST/LandTrendr/CCDC). Seasonal
  mismatch is this skill's own confounder; a documented datum with a stated
  accuracy budget is settled comparability. Keep both.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Change Detection & Spatio-temporal Analysis

Purpose: separate real surface change from the four great impostors —
misregistration, radiometric drift, phenology, and classification error.
Every method below exists to control one of them; skipping the controls
produces confident maps of nothing.

## Preconditions (where change detection is won or lost)

Preconditions 1 and 2 are *checked* here but *established* elsewhere. When
either fails, the owning skill leads and this skill resumes once comparable
observations exist. Precondition 3 is this skill's own problem and is never a
reason to route away.

1. **Co-registration**: sub-pixel alignment between dates (AROSICS or
   manual tie-points). Half a pixel of shift creates edge-shaped phantom
   change everywhere. Verify: flicker-compare crisp features. For elevation
   surfaces or point clouds, vertical datum agreement, co-registration and
   the vertical-accuracy budget belong to `point-cloud-lidar` — a datum
   offset is not subsidence.
2. **Radiometric consistency**: same processing level (surface
   reflectance), same sensor, same processing baseline. If any of the three
   differ, this is a harmonization problem, not a thresholding one: hand it
   to `remote-sensing-analysis` (HLS for Landsat↔Sentinel-2, relative
   normalization with PIFs, `BOA_ADD_OFFSET` across the Sentinel-2 2022
   baseline change).
3. **Same season / phenological stage** for bi-temporal work — a May vs
   September pair "detects" summer. If season can't be matched, use
   composites or time-series methods instead.
4. **Cloud/shadow masks intersected across dates**; analyze only mutually
   valid pixels and report that coverage %.

## Method selection

| Situation | Method |
|---|---|
| Two dates, continuous "how much" | Index differencing (ΔNDVI, ΔNBR...) with statistical thresholding |
| Two dates, categorical "from-what-to-what" | Post-classification comparison (only with strong classifiers) |
| Two dates, multivariate robust | Change vector analysis (CVA); MAD/iMAD for sensor-robust detection |
| Dense stack, gradual + abrupt | Trend + break analysis (BFAST/LandTrendr/CCDC family; at archive scale → `google-earth-engine`) |
| Structure change (buildings) | DL bi-temporal segmentation (siamese U-Net) → `geo-deep-learning` |
| SAR pairs (clouds, disasters) | Log-ratio of calibrated backscatter + speckle handling |
| Vector vintages (parcels, buildings) | Geometry+attribute diff with tolerance (below) |

## Thresholding — never eyeball it

Difference images need a defensible threshold: μ ± k·σ on the difference
histogram (report k), Otsu when bimodal, or supervised thresholds
calibrated on labeled change/no-change samples. Deliver the histogram with
the chosen cut marked. Sensitivity: report changed-area at k-0.5 and
k+0.5; if the story flips, the detection is fragile — say so.

## Post-classification comparison (PCC) — handle with care

PCC error compounds: two 90%-accurate maps yield ≤ ~81% change accuracy,
and biased errors create systematic false transitions. Rules:

- Use ONE classifier trained on both dates' imagery (same legend, same
  features) rather than two independent legacy maps.
- Build the full **transition matrix** (from-class → to-class areas), not
  just a change/no-change binary — impossible transitions (water→forest
  in 1 year) are your error detector.
- Apply a minimum mapping unit consistent across dates before differencing.

## Time-series (dense stack) analysis

- Build a gap-filled, cloud-masked index stack (xarray, time dimension).
- Decompose trend + seasonality + breaks; per-pixel linear trends need
  significance testing (Mann-Kendall + Sen's slope for monotonic trends —
  and FDR correction across millions of pixels, or your "greening map" is
  noise).
- Label break DATES, not just presence — timing is usually the analytic
  payload (when did clearing start?).
- Validate detected breaks against known events (fires, construction
  permits, disaster dates) wherever records exist.

## Vector change audit (two vintages of the same layer)

- Match features by stable ID if it exists; else spatial matching with IoU
  threshold (report it).
- Classify: added / removed / geometry-changed (area delta > tolerance) /
  attribute-changed. Tolerances absorb digitization jitter — 1-2 m for
  cadastre-grade, more for digitized-from-imagery.
- Sum area deltas by class and reconcile totals; unexplained residual =
  matching bugs.

## Accuracy assessment (the deliverable's spine)

Change is rare, so random sampling wastes effort on stable pixels — use
**stratified sampling** (strata: change/no-change or per-transition) with
good-practice area estimation (Olofsson et al. protocol): report
user's/producer's accuracy per stratum AND **area estimates with
confidence intervals** adjusted for map error. A raw pixel count of the
change map is a biased area estimate — always say the adjusted number.

## Reporting template

```
## Change: <phenomenon>, <T1> → <T2 or period>
- Data: <sensor/level>, co-registration RMSE: <px>, valid overlap: <%>
- Method: <...> threshold/params: <...> (sensitivity: <stable/fragile>)
- Transitions: <matrix or top-5 list with areas ± CI>
- Accuracy: stratified n=<>, UA/PA per class, adjusted areas ± CI
- Impostor controls: season <matched?>, radiometry <harmonized?>
```

## Pitfalls checklist

- Phantom edge-change from misregistration.
- Seasonal difference sold as land cover change.
- PCC with two independently produced legacy maps.
- Threshold chosen "because it looked right", no sensitivity.
- Raw changed-pixel counts reported as area (no error-adjusted estimate).
- Trend maps without multiple-testing control.
- SAR change on unfiltered linear-power images.

## Execution contract

- **Workflow:** define the change question; harmonize extent, season, radiometry, resolution, and registration; select method; estimate change; validate; report uncertainty.
- **Decision rules:** use direct differencing only for comparable continuous signals, post-classification comparison for stable class legends, and time-series methods when a dense temporal stack exists.
- **Verification protocol:** quantify co-registration, valid overlap, threshold sensitivity, transition accounting, and accuracy-adjusted area with confidence intervals.
- **Failure modes:** reject causal change claims when season, sensor, clouds, registration, or independent map errors can explain the signal.
- **Deliverables:** change map, transition or trend table, parameter record, validation sample and metrics, adjusted-area estimate, and limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before using version-sensitive products or APIs and record the checked date.

Referenced files: 2

View saved version →

---
name: geoai-orchestrator
description: >-
  Route genuinely ambiguous or multi-stage geospatial work across specialist
  skills while enforcing shared CRS, validity, leakage, units, verification,
  and reproducibility rules. Use for requests spanning multiple stages such as
  acquisition, imagery, modeling, analysis, and map delivery, or for an
  explicit end-to-end pipeline. Never invoke for one domain merely because a
  parameter is unclear. Code implementation/review, backend or platform
  choice, and production-readiness review are direct specialist tasks. Do not
  add this skill as a layer around one specialist.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# GeoAI Orchestrator

The hub of an 18-skill geospatial module. Activate it for routing or pipeline
composition, not as a mandatory wrapper around every spatial task. Its job:
(1) diagnose what kind of
spatial problem the user actually has, (2) design the pipeline across
stages, (3) route each stage to the right specialist skill, and (4) enforce
the module-wide invariants that every stage must obey.

## Routing gate — read before producing any output

This orchestrator routes by **invoking**, never by naming. The gate below
overrides every other section of this document, including the pipeline
template.

1. **Invoke, do not list.** Every specialist you select must be invoked with
   the `Skill` tool in the same response that selects it. Naming a skill in a
   table, plan, or prose sentence is not a handoff. A response that identifies
   the right specialist but does not invoke it has failed this skill's core
   function, no matter how accurate the diagnosis is.
2. **Route every correction, not the first one.** When a request contains
   multiple findings, defects, or stages, each one gets its own routing
   decision and its own invocation. Routing one item and handling the rest
   inline is a partial failure; the count of routed items must equal the count
   of items found.
3. **Never make routing conditional on permission.** Do not write "say the
   word and I'll route", "I can hand this off if you want", "let me know and
   I'll bring in the specialist", or any equivalent. Offering to route later is
   the single most common failure of this skill. If you have identified the
   specialist, invoke it now.
4. **Clarification is not a substitute for routing.** Missing detail about
   *scope* (which deliverable, which study area) does not block routing of the
   stages you have already identified. Ask the scope question and route in the
   same response. Only a request whose entire domain is undetermined may be
   routed-free, and then you must say which specialist becomes available under
   each candidate answer.
5. **Audit requests are `deliver` requests.** "Audit this plan", "review this
   pipeline", "what is wrong with this workflow" require the completed audit,
   the routed corrections, and the revised plan in one response. Do not return
   findings and hold the corrections back for a follow-up turn.

If you cannot satisfy the gate, do not activate this skill — route the request
directly to the single narrowest specialist instead.

## Module map — route by problem type

| Stage / problem | Specialist skill |
|---|---|
| Data acquisition, formats, CRS, tiling, pipelines | `geo-data-engineering` |
| Satellite/aerial imagery, spectral indices, classification | `remote-sensing-analysis` |
| Planetary-scale archives, GEE Python API, cloud compositing | `google-earth-engine` |
| CNN/U-Net/ViT on EO data, segmentation, detection | `geo-deep-learning` |
| Autocorrelation, hotspots, clusters, spatial regression | `spatial-statistics` |
| Site selection, suitability, AHP/weighted overlay | `mcda-suitability-analysis` |
| Interpolation from point samples, kriging, variograms | `geostatistics-interpolation` |
| DEM, slope, watersheds, flow, viewshed | `terrain-hydrology` |
| LiDAR / point clouds, DTM/DSM/CHM, PDAL | `point-cloud-lidar` |
| Routing, service areas, accessibility, OD matrices | `network-accessibility-analysis` |
| GPS tracks, trajectories, stops/trips, map matching | `movement-trajectory` |
| Multi-temporal comparison, land cover change, trends | `change-detection` |
| Map design, choropleths, web maps, publication figures | `cartography-geoviz` |
| Spatial SQL, PostGIS, large-scale spatial joins | `postgis-spatial-sql` |
| Local ArcGIS Pro, ArcPy, `.aprx`, or `.gdb` execution | `arcgis-pro-automation` |

This table selects specialists; it does not hand off to them. Every row you
select must be invoked under the routing gate. For cross-cutting method
standards (leakage, metrics, reproducibility), invoke `ml-experiment-standards`
and `swe-devops-standards` when their rules apply.

## Pipeline design protocol

For any multi-stage request, produce a short pipeline plan BEFORE writing
code, then invoke the specialists that plan names in the same response:

```
## Pipeline: <goal>
1. <stage> → <skill> → output: <artifact> → check: <verification criterion>
2. ...
Success criterion: <what the user can inspect to accept the result>
```

The plan is a routing manifest, not a proposal awaiting approval. Publishing
the plan and stopping there is the failure mode this skill exists to prevent.
Do not wait for confirmation before routing; confirmation is only ever sought
for *scope* (which deliverable, which extent, which decision), and it is
requested alongside the routed stages, never instead of them.

Every stage ends with a verification criterion. Spatial work fails silently
(wrong CRS, empty joins, inverted axes produce plausible-looking garbage),
so a stage without a check is not a stage.

## Module-wide invariants (enforced in every stage)

1. **CRS is explicit, always.** Report the CRS of every input on first
   contact. Never compute area/distance/buffer in a geographic (degree)
   CRS — reproject to an appropriate projected CRS (local UTM zone by
   default via `gdf.estimate_utm_crs()`; equal-area such as EPSG:6933 for
   global area statistics). If a CRS is undefined, stop and resolve it;
   never guess silently.
2. **Axis order discipline.** GeoJSON is lon/lat; many APIs and humans say
   lat/lon. Verify with a known landmark before pipeline-scale processing.
3. **Geometry validity before analysis.** Check `is_valid`; repair with
   `shapely.make_valid` (not `buffer(0)`, which can silently drop parts).
4. **Row-count accounting.** After every join/overlay/filter, report rows
   in vs rows out. Silent duplication or loss is the top geospatial bug.
5. **Spatial autocorrelation awareness.** Random train/test splits on
   spatial data leak. Any ML stage follows the canonical protocol in
   `ml-experiment-standards` → `references/spatial-cv-protocol.md`.
6. **Units in column names.** `area_ha`, `dist_km`, `elev_m` — never bare
   `area`. Unit confusion survives code review; column names don't lie.
7. **Visual + numeric verification.** Every spatial output gets both a
   summary table AND a quick map check (`.explore()`, a PNG, or GIS
   software). A confusion matrix cannot show spatially clustered errors.
8. **Reproducibility.** Pin package versions, seed randomness, log
   parameters. Intermediate artifacts go to GeoPackage or GeoParquet, never
   shapefile (10-char column truncation, 2 GB limit, no proper encoding).

## Internationalization note

Attribute tables in non-ASCII locales break naive string handling.
Canonical example: Turkish dotted/dotless I — `'İ'.lower()` yields a
2-character string in Python. Before any string matching on attributes,
apply a locale-aware normalization step and show `value_counts()` of
cleaned categorical fields. Prefer UTF-8 formats; legacy shapefiles may
carry cp1252/cp125x mojibake silently.

## Choosing the stack

Default to the open Python stack: GeoPandas + Shapely 2 + Rasterio +
xarray/rioxarray + PyProj. Route to PostGIS when data exceeds comfortable
memory (~millions of features) or needs concurrent/repeated querying; to
Earth Engine when the data is a planetary archive rather than local files.
Use GDAL CLI for bulk format conversion. If the user works in ArcGIS Pro or
QGIS, generate headless-runnable scripts (arcpy / PyQGIS) rather than click
instructions, and keep the analysis logic portable.

## Anti-patterns to catch early

- Buffering in degrees ("0.01 degree buffer") — reproject first.
- `EPSG:4326 → Web Mercator` area statistics — Mercator distorts area
  massively away from the equator.
- Joining datasets from different CRS without alignment.
- Treating a DEM's nodata value (-9999, 3.4e38) as real elevation.
- Classifying imagery without checking cloud/shadow masks.
- Reporting model accuracy without a spatially independent test set.

## Execution contract

- **Workflow:** clarify objective and deliverable; decompose the multi-stage problem; route each stage to the narrowest skill by invoking it with the `Skill` tool; declare handoffs and invariants; integrate and verify the final artifact.
- **Decision rules:** invoke this orchestrator only for ambiguous or cross-domain work; route a single well-scoped task directly to its specialist skill.
- **Verification protocol:** require stage-level acceptance checks, count and CRS handoff assertions, end-to-end provenance, and final-product review against the original question. Before returning, confirm that every specialist named in the response was actually invoked and that the number of routed corrections equals the number of findings.
- **Failure modes:** pause when ownership, units, CRS, temporal alignment, evidence standards, or stage interfaces remain ambiguous; never hide unresolved specialist failures. Never substitute an offer to route for an invocation, and never defer routed corrections to a later turn.
- **Deliverables:** pipeline plan, skill-routing table, stage inputs and outputs, verification gates, risk register, and final integration checklist.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) and the selected specialists' registries before fixing interfaces.

Referenced files: 2

View saved version →

---
name: geo-data-engineering
description: >-
  Always invoke when geospatial data must be acquired, prepared, repaired,
  scaled, or moved through a repeatable pipeline. Covers open-data/OSM/STAC
  acquisition, spatial formats, CRS transforms, quality checks, and batch ETL
  architecture for growing or recurring joins. Invoke alongside PostGIS for
  database execution and alongside SWE standards when code is delivered. Do
  not trigger merely because another specialist reads analysis-ready data.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Geospatial Data Engineering

Purpose: get spatial data into a clean, validated, analysis-ready state with
a repeatable pipeline — the stage where most real-world GIS time is spent
and most silent errors are born.

## Format selection

| Format | Use for | Avoid because |
|---|---|---|
| **GeoParquet** | Analysis interchange, big vector, columnar workflows | Not yet readable by some legacy desktop GIS |
| **GeoPackage** | Desktop GIS exchange, multi-layer projects | Slower than Parquet at scale; SQLite locking |
| **FlatGeobuf** | Streaming, HTTP range reads | Single layer |
| **COG** (Cloud-Optimized GeoTIFF) | All raster deliverables | — (make every GeoTIFF a COG) |
| **Zarr/NetCDF** | Multi-dimensional (time × band × y × x) | Overkill for single rasters |
| Shapefile | Only when a legacy tool demands it | 10-char columns, 2 GB cap, encoding chaos, multi-file fragility |
| CSV + WKT/lon-lat | Simple point exchange | No CRS metadata — document it explicitly |

## Acquisition playbook

- **OpenStreetMap**: small areas → `osmnx`; large extracts → Geofabrik PBF +
  `pyrosm`/`osmium`. Respect tag heterogeneity: always inspect tag value
  distributions before filtering.
- **Buildings/places at scale**: Overture Maps (GeoParquet on S3/Azure,
  query with DuckDB spatial — often the fastest path).
- **Satellite/raster**: STAC APIs via `pystac-client` + `odc-stac` — see
  `remote-sensing-analysis`; planetary archives → `google-earth-engine`.
- **Boundaries**: authoritative national source first; Natural Earth / GADM /
  geoBoundaries for global work — record which, versions differ materially.
- Record every acquisition: source URL, query parameters, retrieval date,
  license. Put it in a `DATA_SOURCES.md` next to the data.

## CRS engineering

- Store in EPSG:4326 or source CRS; **analyze** in a projected CRS suited to
  the extent: local UTM zone (`gdf.estimate_utm_crs()`), national grid, or
  equal-area (EPSG:6933/Mollweide) for cross-region area stats.
- Datum shifts matter at sub-meter precision: transformations between datums
  need the right transformation grid (`pyproj.network.set_network_enabled(True)`
  when accuracy matters).
- Never strip or overwrite a CRS to "fix" misaligned layers — diagnose which
  layer is wrong with a known landmark instead.

## Cleaning pipeline

Run `scripts/clean_vector.py` (or import its `clean_vector()` function) as
the standard hygiene pass: drops empty/null geometries, repairs invalid ones
with `make_valid`, de-duplicates, reprojects, and **prints an accounting
report** so silent data loss is impossible.

Then: normalize text attributes (trim, collapse whitespace, locale-aware
casefold — beware Turkish İ/ı, German ß), coerce dtypes explicitly, and show
`value_counts()` of every categorical you will later filter on.

## Scale strategies

- **Fits in RAM**: GeoPandas + Shapely 2 vectorized ops. Ensure the spatial
  index is used (`sjoin`, `query_bulk`) — hand-rolled loops are O(n²).
- **Bigger than RAM, single machine**: DuckDB `spatial` extension over
  GeoParquet (predicate pushdown + spatial SQL), or `dask-geopandas`.
- **Served / concurrent / transactional**: PostGIS — see `postgis-spatial-sql`.
- Rasters: windowed reads (`rasterio.windows`), chunked xarray + dask;
  never `read()` a 50 GB mosaic into memory.

## Pipeline standards

- Idempotent steps with explicit inputs/outputs on disk; re-running never
  corrupts state.
- Checkpoint after expensive stages (download, big join) in GeoParquet/GPKG.
- Log an accounting line per stage: rows/features/pixels in → out.
- Deterministic ordering before writing (sort by stable key) so diffs are
  meaningful.

## Pitfalls checklist

- CSV opened without declaring lon/lat columns' CRS.
- Shapefile column names silently truncated on export.
- Encoding mojibake from legacy files (try `encoding="utf-8"` then cp1252).
- Mixed geometry types in one layer (Polygon + MultiPolygon breaks some
  tools — normalize with `.explode()` or promote to Multi*).
- Antimeridian and pole-crossing geometries after naive reprojection.
- Downloaded "latest" data with no recorded version/date — unreproducible.

## Execution contract

- **Workflow:** inventory sources and contracts; acquire with provenance; inspect CRS, schema, geometry, and scale; clean deterministically; validate; write an analysis-ready artifact.
- **Decision rules:** select formats and engines from size, geometry, concurrency, and downstream access needs; never infer CRS or destructive repairs silently.
- **Verification protocol:** reconcile feature or pixel counts at every stage, assert CRS and geometry invariants, sample outputs spatially, and rerun to confirm idempotence.
- **Failure modes:** quarantine ambiguous CRS, mixed units, invalid encodings, lossy format conversions, or unexplained row loss instead of guessing.
- **Deliverables:** validated dataset, machine-readable schema and CRS, provenance manifest, accounting log, rejected-record report, and reproducible pipeline.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before using version-sensitive formats or APIs and record the checked date.

Referenced files: 3

View saved version →

---
name: geo-deep-learning
description: >-
  Invoke before recommending, training, or auditing a neural method for
  geospatial imagery, including vision transformers, U-Net/DeepLab/SegFormer,
  object detection, pixel classification, building/road extraction, and EO
  foundation-model fine-tuning. Also invoke for neural chip-split validity,
  IoU/accuracy claims, augmentation, imbalanced losses, spatial validation,
  or sliding-window inference. Use remote-sensing-analysis for non-neural
  methods and change-detection when temporal change is the deliverable.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Geospatial Deep Learning

Purpose: deep learning on Earth observation with the two failure modes that
dominate this field designed out from the start: **spatial leakage**
(inflated metrics from nearby train/test pixels) and **georeferencing loss**
(predictions that no longer align with the map).

## Characterise the label set before naming an architecture

Architecture advice given without knowing the label set is guesswork. Before
recommending U-Net versus a foundation model versus a non-deep baseline, state
or ask for:

- **Label count and labelled area** — polygons alone say nothing; 40 polygons
  covering 2 ha and 40 covering 2 000 km² are different problems.
- **Geographic spread** — are the labels clustered in one scene, one season and
  one sensor, or distributed across the deployment domain? Clustered labels cap
  what any model can generalise to, and they decide whether a geographically
  independent validation split is even constructible.
- **Class balance and minority-class pixel fraction**, so loss and sampling
  choices are grounded rather than assumed.
- **Deployment geography** — where predictions will be made, relative to where
  the labels are.

Do not answer "fine-tune a large model or use a simpler approach" before these
are known. When the user has not supplied them, ask and give the provisional
recommendation *conditioned on* the answers ("if the 40 polygons sit in one
scene, then …; if they span the region, then …"), never a single unconditional
recommendation.

## Problem framing first

| Task | Head/architecture default | Metric |
|---|---|---|
| Pixel-wise classes (land cover) | U-Net / DeepLabv3+ (pretrained encoder) | mIoU, per-class IoU |
| Binary extraction (buildings, water, roads) | U-Net + Dice/CE hybrid | IoU, F1; boundary F1 for roads |
| Object detection (vehicles, ships, trees) | YOLO-family / Faster R-CNN, rotated boxes if oriented | mAP@50 |
| Scene classification | Fine-tuned CNN/ViT | F1 (macro) |
| Regression (height, biomass, density) | U-Net with regression head | RMSE/MAE + spatial residual map |

Before any deep model: run a cheap baseline (random forest on bands+indices,
or thresholded index). If the DL model can't beat it clearly, the problem is
data, not architecture. `segmentation-models-pytorch` and `torchgeo` cover
most needs — don't hand-build architectures without a reason.

## Chipping (dataset construction)

- Chip size: 256–512 px; stride < chip size only for training (overlap
  augments), never let overlapping chips straddle the train/val boundary.
- **Preserve georeferencing**: store each chip's transform/bounds (torchgeo
  datasets or a sidecar index in GeoParquet). A prediction you can't put
  back on the map is worthless.
- Keep chips in the native data range; normalize with **dataset-computed**
  per-band statistics (ImageNet stats only for 3-band RGB with a pretrained
  encoder, and say so).
- Class imbalance is the norm (buildings ≈ 2-5% of pixels). Log per-chip
  class fractions; oversample positive-containing chips rather than
  distorting the loss beyond recognition.

## Split policy — the non-negotiable

Split by **geographic block or scene**, never by random chip. Adjacent
chips are near-duplicates; random splits produce beautiful, fake validation
curves. Follow the canonical protocol:
`ml-experiment-standards` → `references/spatial-cv-protocol.md`.
For generalization claims across regions, hold out an entire region.

## Training defaults

- Loss: Dice + CE (segmentation, imbalanced); plain CE when balanced; Focal
  only after comparing — it's not a free win.
- Augmentation: flips/rot90 are safe for nadir imagery; be careful with
  color jitter on multispectral (it breaks radiometric meaning — prefer
  band dropout or slight scaling); never augment in ways that violate the
  physics.
- Encoder pretrained; multispectral input → inflate/replace first conv, or
  use an EO foundation model checkpoint (Prithvi, SatMAE, Clay) when bands
  match.
- Early stopping on val mIoU (patience 10-15); cosine or plateau LR
  schedule; AMP on by default.
- Log config + metrics + git hash per run — see `ml-experiment-standards`.

## Inference on large scenes

Sliding window with overlap (25-50%) and blending (feather/gaussian or
center-crop stitching) to kill tile-edge artifacts. Then:

```python
import rasterio

with rasterio.open(scene_path) as src:
    profile = src.profile
profile.update(count=1, dtype="uint8", nodata=255, compress="deflate")
with rasterio.open(out_path, "w", **profile) as dst:
    dst.write(mask.astype("uint8"), 1)  # same transform/CRS as the scene
```

Post-process: sieve tiny blobs (min mapping unit), optionally regularize
building polygons, and vectorize (`rasterio.features.shapes`) for GIS
delivery. Report metrics AFTER post-processing too — that's what the user
ships.

## Verification protocol

1. Metrics table: per-class IoU/F1 with CI across seeds or folds.
2. **Error map**: prediction vs reference overlaid on imagery for 3+
   representative areas including a known-hard one.
3. Sanity inference on an out-of-distribution patch (different season/
   region) with an honest note on degradation.
4. Alignment check: overlay predictions on the source scene in a GIS at
   two zoom levels — catches transform bugs instantly.

## Pitfalls checklist

- Random chip split → leaked, unreproducible "SOTA".
- Normalizing test data with train-time stats not saved → skewed inference.
- Losing the geotransform in NumPy-land; writing predictions with default
  north-up transform.
- Tile-edge seams from no-overlap inference.
- uint16 imagery fed to a float pipeline without scaling → dead gradients.
- Accuracy reported on chip level while the product is a stitched map.

## Execution contract

- **Workflow:** frame target and unit of prediction; build chips and labels; create spatial splits; train against a baseline; run overlap-aware inference; validate the stitched product.
- **Decision rules:** use deep learning only when label volume, spatial texture, compute, and expected uplift justify it; otherwise prefer a simpler remote-sensing or ML workflow.
- **Verification protocol:** report spatial holdout metrics across seeds or folds, inspect error maps and hard areas, test geographic transfer, and check output georeferencing.
- **Failure modes:** invalidate results for leaked chips, label misalignment, train/inference normalization drift, tile seams, or metrics computed at the wrong product unit.
- **Deliverables:** model and configuration, split manifest, preprocessing contract, metrics with uncertainty, georeferenced predictions, error maps, and model card limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before selecting framework APIs, datasets, or weights and record the checked date.

Referenced files: 2

View saved version →

---
name: geostatistics-interpolation
description: >-
  Turn scattered point measurements into continuous surfaces with quantified
  uncertainty: variogram modeling, ordinary/universal/regression kriging,
  IDW, and spatially honest cross-validation. Use when unobserved values must
  be estimated from sparse samples such as stations, wells, or soundings.
  Trigger on "interpolate", "kriging", "variogram", or "IDW"; a named
  interpolation method is sufficient even when the requested surface is
  described informally as a heatmap. Do not use for point-density heatmaps,
  zonal aggregation, or raster resampling without value interpolation.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Geostatistics & Interpolation

Purpose: interpolation that reports what it doesn't know. The difference
between a professional product and a pretty raster is the uncertainty
surface and an honest cross-validation — both are non-optional here.

## Method selection

| Situation | Method |
|---|---|
| Dense, smooth phenomenon, quick look | IDW (report power parameter; test 1-3) |
| Physical phenomenon with spatial structure, need uncertainty | **Ordinary kriging** (default professional choice) |
| Clear trend (elevation gradient, coastal effect) | Universal kriging or regression kriging on covariates |
| Strong covariates available (DEM, land cover, distances) | Regression kriging / random-forest residual kriging |
| Categorical target | Indicator kriging |
| Honeycomb-free tessellation, no extrapolation wanted | Natural neighbor |

IDW is a reasonable baseline but has no error model and bullseyes around
extremes; say so when delivering IDW-only products.

## Exploratory phase (before any interpolation)

- Map the points with values; look for duplicates at identical coordinates
  (average or offset them — kriging matrices go singular otherwise).
- Histogram + skew: strongly skewed variables (rainfall, concentrations)
  usually want a log/normal-score transform; back-transform predictions
  properly (bias correction for lognormal kriging).
- Trend check: regress value on x, y, and candidate covariates; visible
  trend → universal/regression kriging path.
- Declustering if sampling is preferential (dense where values are high).

## Variogram discipline

The variogram is a MODELING decision, not an auto-fit output:

```python
import gstools as gs

bin_center, gamma = gs.vario_estimate((x, y), values, max_dist=dmax)  # dmax ≈ half extent
model = gs.Exponential(dim=2)
model.fit_variogram(bin_center, gamma, nugget=True)
print(model)   # report: nugget, sill, range — plus the fitted plot
```

- Max lag ≈ half the domain diameter; ≥30 pairs per bin.
- Check anisotropy with directional variograms (0/45/90/135°); geological
  and meteorological fields are often anisotropic — fit an anisotropic
  model rather than ignoring it.
- Interpret and report the parameters in words: nugget (measurement error +
  micro-scale variance), range (correlation distance), sill. A nugget near
  the sill means the data barely support interpolation — say that honestly.
- Never interpolate meaningfully beyond the variogram range from the
  nearest sample; mask or flag those cells.

## Kriging execution

```python
from pykrige.ok import OrdinaryKriging

ok = OrdinaryKriging(x, y, values, variogram_model="exponential",
                     variogram_parameters={"sill": s, "range": r, "nugget": n},
                     coordinates_type="euclidean")
z, ss = ok.execute("grid", gridx, gridy)   # ss = kriging VARIANCE — keep it!
```

- Work in a projected CRS (kriging distances in degrees are wrong except
  with explicitly geographic-aware models).
- Deliver TWO rasters: prediction AND kriging standard deviation
  (`np.sqrt(ss)`). The uncertainty map drives where to sample next and
  where not to trust the map.
- Search neighborhood: 16-32 nearest points typical; document it.

## Validation — leave-one-out and beyond

- **LOOCV** for small n: report ME (bias ≈ 0?), RMSE, and standardized
  RMSE (should be ≈ 1 if kriging variances are honest).
- For clustered samples, LOOCV flatters; also run spatial block CV
  (canonical recipe: `ml-experiment-standards` →
  `references/spatial-cv-protocol.md`) and report both.
- Scatter plot observed vs predicted with 1:1 line; map the CV residuals —
  spatially clustered residuals mean missing trend/covariate.
- Compare against the dumb baseline (global mean, IDW): kriging must earn
  its complexity.

## Reporting template

```
## Surface: <variable>
- n = <>, extent, CRS, transform applied: <log/none>
- Variogram: <model>, nugget/sill/range = <...>, anisotropy: <...>
- Method: <OK/UK/RK + covariates>
- LOOCV: ME <>, RMSE <>, RMSE_std <>; block-CV RMSE <>
- Uncertainty: kriging SD raster delivered; area beyond reliable range masked
```

## Pitfalls checklist

- Auto-fitted variogram accepted without looking at the plot.
- Kriging in EPSG:4326 over large extents.
- Prediction raster delivered without its variance/SD companion.
- Log-kriging back-transformed by naive `exp()` (bias!).
- Extrapolation far beyond the data hull presented at equal confidence.
- Duplicated station coordinates crashing or silently distorting the fit.

## Execution contract

- **Workflow:** inspect sampling and transform needs; model spatial dependence; fit candidate surfaces; validate against simple baselines; quantify uncertainty; mask unsupported extrapolation.
- **Decision rules:** choose IDW or trend methods for transparent baselines, kriging when a defensible variogram exists, and regression kriging when covariates improve blocked validation.
- **Verification protocol:** examine variogram diagnostics, LOOCV and spatial block-CV residuals, baseline uplift, uncertainty calibration, and behavior beyond the sample hull.
- **Failure modes:** withhold a confident surface when sampling is too sparse, dependence is absent, duplicates dominate, distance units are invalid, or uncertainty is uncalibrated.
- **Deliverables:** prediction surface, uncertainty surface, variogram and parameters, validation table, sampling-support mask, and method limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before relying on package APIs or defaults and record the checked date.

Referenced files: 2

View saved version →

---
name: google-earth-engine
description: >-
  Invoke when Earth Engine, GEE, ee., or geemap is named; when work needs its
  server-side catalog; or when choosing Earth Engine versus local xarray or
  desktop processing for a large area or long archive. Covers image
  collections, masking, compositing, reducers, zonal statistics, time series,
  classification, quota-aware batching, and exports. This is an execution
  platform skill; combine it with remote-sensing-analysis or change-detection
  when those skills own the scientific method.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Google Earth Engine

Purpose: use GEE's server-side model correctly. The recurring failure
modes are **client/server confusion** (calling `.getInfo()` in loops,
Python `if` on server objects), **unbounded computation** (timeouts from
unscaled reductions), and **silent default scales** (statistics computed
at the wrong resolution).

## Should this run here at all? — Earth Engine versus local

Answer this before writing any `ee.` code. The decision turns on six things, and
you cannot make it without them, so establish them first — asking alongside a
provisional recommendation, never instead of one:

1. **Archive extent and duration** — area, and how many years at what revisit.
   This is what makes server-side worth its constraints; a single scene does not.
2. **Algorithm expressibility** — can the work be written as masks, reducers and
   band math? Anything needing arbitrary per-pixel iteration, a custom solver, or
   a Python library GEE does not host belongs local.
3. **Data locality and sensitivity** — restricted or offline data cannot be
   uploaded, and that ends the discussion regardless of scale.
4. **Interactive limits versus batch** — see [Quotas and etiquette](#quotas-and-etiquette).
   Anything beyond a ~5 minute interactive request has to be designed as a batch
   export from the start, not retrofitted when `getInfo` times out.
5. **Export volume** — what actually comes back: a few reduced statistics, or
   full-resolution per-pixel stacks you will store and reprocess locally.
6. **Reproducibility cost** — the real price of moving server-side. The catalog
   version can shift under you and the computation leaves no local trace, so
   choosing GEE obliges you to ship the [provenance record](#provenance-record).
   State this cost when you recommend GEE; a recommendation that omits it is
   incomplete.

Recommend Earth Engine only when 1 and 2 favour it and 3 permits it. When the
answer is genuinely balanced, say so and name the deciding question rather than
defaulting to the platform this skill is about. `xee` and STAC + `stackstac` /
`odc-stac` are the middle paths worth naming: catalog access with local compute.

## Mental model — everything is deferred

`ee.Image`, `ee.ImageCollection`, `ee.FeatureCollection` are **server-side
descriptions**, not data. Nothing computes until an output is requested
(`getInfo`, export, map tile). Consequences:

- Never use Python `if`/`for` on server values — use `ee.Algorithms.If`
  sparingly, prefer `.map()` + filters. A Python loop that calls
  `.getInfo()` per element is the #1 GEE performance bug.
- `.getInfo()` blocks and transfers; use it for tiny scalars only.
  Anything sized → **Export** (to Drive/GCS/Asset).
- Debug with `.aggregate_array()`, `.first()`, `.limit(3)` probes — not by
  printing whole collections.

## Canonical pipeline (Sentinel-2 cloud-free composite)

```python
import ee
ee.Initialize(project="my-project")

aoi = ee.Geometry.Rectangle([27.0, 38.3, 27.4, 38.6])

def mask_s2(img):
    # Cloud Score+ is the current best practice (threshold ~0.5-0.65)
    cs = img.linkCollection(csplus, ["cs_cdf"]).select("cs_cdf")
    return img.updateMask(cs.gte(0.6))

csplus = ee.ImageCollection("GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED")
s2 = (ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
      .filterBounds(aoi)
      .filterDate("2025-05-01", "2025-09-30")
      .map(mask_s2))
composite = s2.median().clip(aoi)
ndvi = composite.normalizedDifference(["B8", "B4"]).rename("ndvi")
```

Collection choices: `S2_SR_HARMONIZED` (post-2022 offset harmonized),
`LANDSAT/LC08/C02/T1_L2` + friends (apply scale factors: optical
`*0.0000275 - 0.2`), `MODIS/061/...` for daily/coarse, ERA5-Land for
climate. Record collection IDs + date filters in the deliverable.

## Reducers and zonal statistics — scale is not optional

```python
stats = ndvi.reduceRegions(
    collection=districts,
    reducer=ee.Reducer.mean().combine(ee.Reducer.stdDev(), sharedInputs=True),
    scale=10,                    # ALWAYS explicit — native resolution
    tileScale=4,                 # raise when "computation timed out"
)
```

- `scale` defaults to the map zoom level in some paths — silently coarse
  statistics. Always set it to the data's native resolution (or state the
  deliberate coarsening).
- `bestEffort=True` silently degrades scale to fit limits — avoid in
  analysis; prefer `tileScale` + exports.
- Large reductions → `Export.table.toDrive`, not `.getInfo()`.
- Weighted vs unweighted reducers differ at polygon edges
  (`.unweighted()` for counts of whole pixels); state which you used.

## Time series

- Build per-period composites with a mapped function over
  `ee.List.sequence` of dates (monthly/seasonal medians), then reduce —
  don't export daily stacks you'll aggregate anyway.
- For per-pixel trends: `ee.Reducer.sensSlope()` (robust) or
  `linearFit`; harmonic regression (`.addBands` of sin/cos terms) for
  phenology. Mask by count of valid observations — trends from 4 pixels
  of 200 possible are noise; report the count band.
- For break detection at archive scale (LandTrendr/CCDC available in GEE),
  method selection follows `change-detection`.

## Classification in GEE

`ee.Classifier.smileRandomForest` covers most cases. Training samples via
`image.sampleRegions`; split train/test **spatially** (add a grid-cell
attribute and filter — random `randomColumn` splits leak; see
`ml-experiment-standards` → `references/spatial-cv-protocol.md`). Report
per-class accuracy from `errorMatrix`; area estimates from a classified
map still need design-based adjustment (`change-detection` / Olofsson).

## Exports and hand-off

- `Export.image.toDrive/toCloudStorage` with explicit `region`, `scale`,
  `crs`, `maxPixels`; use `crsTransform` when pixel alignment with an
  existing raster matters.
- Export > ~10⁸ pixels: shard by tiles or use `toAsset` intermediate.
- Hand off to the local Python stack (rasterio/xarray) via COG exports, or
  `xee` for xarray-native access; visualize interactively with `geemap`.

## Quotas and etiquette

Batch tasks queue (check task status; don't fire hundreds blindly).
Interactive requests time out at ~5 min — long jobs go to batch export.
Cache intermediate products as assets when a pipeline reuses them.

## Provenance record

Server-side computation is invisible after the fact: the catalog moves under
you, a reducer default changes the number, and nothing in the exported file
says which archive produced it. Every Earth Engine deliverable ships with a
provenance record, emitted as a sidecar JSON next to the export — not left
in the notebook:

- **Catalog asset IDs with their version suffix** (`COPERNICUS/S2_SR_HARMONIZED`
  and the specific collection version), plus the date range and filters applied.
- **Mask method and thresholds** — cloud probability source, threshold value,
  and any morphological buffer.
- **Reducers and their arguments**, including `tileScale`, `bestEffort`, and
  any `crsTransform`.
- **Export parameters**: `region`, `scale`, `crs`, `maxPixels`, and the task ID.
- **Run date and the `ee.__version__` / API client version**, because
  server-side defaults change without notice.

Recommending Earth Engine over a local workflow is incomplete without this:
the reproducibility cost is the main thing the user trades away by moving
server-side, so state how it is recovered.

## Verification protocol

1. Probe: `composite.select("B4").projection().nominalScale().getInfo()`
   and band names — confirms scale/CRS assumptions before reductions.
2. Visual check in geemap at 2 zoom levels vs a basemap.
3. Cross-check one zonal statistic against a local computation on an
   exported clip (catches scale/masking discrepancies).
4. Report: collection IDs, date ranges, mask method + threshold, scale,
   reducer types.

## Pitfalls checklist

- `.getInfo()` inside a loop (move logic server-side).
- Missing `scale` in reduceRegion(s) → zoom-dependent statistics.
- Landsat C2 used without scale factors → reflectance > 1.
- `bestEffort=True` hiding resolution degradation.
- Median composite including cloudy pixels (mask BEFORE reduce).
- Python conditionals on server-side objects (always false-y).
- Trend maps without valid-observation-count masking.

## Execution contract

- **Workflow:** define collection and period; build a server-side mask and transform pipeline; test on a small region; compute; verify scale and projection; export reproducibly.
- **Decision rules:** use Earth Engine for planetary archives and scalable aggregation, local tools for sensitive or offline data, and batch exports for work beyond interactive limits.
- **Verification protocol:** probe bands, projection, scale, masks, and observation counts; inspect spatial samples; cross-check one exported statistic locally; record collection versions and parameters.
- **Failure modes:** stop for client-side loops, implicit scale, masked-pixel bias, quota-driven silent degradation, expired assets, or unbounded region operations.
- **Deliverables:** runnable script, collection and date manifest, mask and reducer parameters, task/export settings, verification evidence, and exported asset inventory.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) at execution time for catalog, API, quota, and policy changes.

Referenced files: 2

View saved version →

---
name: mcda-suitability-analysis
description: >-
  Always invoke for spatial suitability, site selection, AHP, criteria
  weights, or weighted-overlay work, including audits of inconsistent
  pairwise judgments and requests for only a final map. Covers consistency,
  standardization, constraints, ranked surfaces, shortlists, and sensitivity.
  Route travel-time placement and location-allocation to
  network-accessibility-analysis.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# MCDA & Suitability Analysis

Purpose: produce suitability maps whose weights, scales, and assumptions are
explicit, consistent, and stress-tested. A suitability map without a
sensitivity analysis is an opinion with a legend.

## Workflow

1. **Structure**: goal → criteria (factors) → constraints. Constraints are
   binary masks (legal exclusions, water bodies, slope > threshold) applied
   at the END by multiplication; factors are continuous and weighted.
   Keep them apart — encoding a constraint as a heavily-weighted factor is
   a classic error that lets forbidden areas score "acceptable".
2. **Criteria layers**: each factor as a raster on a COMMON grid (same CRS,
   extent, cell size, snap). Resample categorical layers with nearest,
   continuous with bilinear; document each.
3. **Standardization** to a common suitability scale (0-1 or 0-255):
   - Linear min-max for monotonic "more is better/worse".
   - Fuzzy membership (sigmoid/linear with control points) when suitability
     saturates — justify control points from domain knowledge.
   - Categorical layers: explicit reclass table, shown to the user.
   Direction check: confirm for EVERY layer whether high raw value means
   high or low suitability (slope: low=good; distance-to-road: usually
   low=good). Direction bugs survive to the final map invisibly.
4. **Weights** (AHP below, or direct/ranked methods with rationale).
5. **Aggregation**: weighted linear combination (WLC) default; OWA when
   the decision-maker's risk attitude (AND-like vs OR-like) matters.
6. **Constraint mask** multiply; classify the result (equal interval or
   quantiles — say which and why); **sensitivity analysis**; validate
   against known good/bad sites if any exist.

## AHP with consistency enforcement

Pairwise comparisons on Saaty's 1-9 scale; weights from the principal
eigenvector; consistency ratio (CR) must be < 0.10 or the matrix goes back
for revision. Run `scripts/ahp_weights.py` to compute weights + CR from a
reciprocal comparison matrix (it validates reciprocity and reports λ_max).

Practices: elicit comparisons pair by pair with verbal anchors ("moderately
more important" = 3); with multiple experts, aggregate judgments by
geometric mean BEFORE computing weights; report the full matrix, weights,
λ_max and CR in the deliverable. If CR ≥ 0.10, identify the most
inconsistent triad and ask the expert to revisit it — do not silently
massage numbers.

## Aggregation

```python
suit = np.zeros_like(factors[0], dtype="float32")
for w_i, f in zip(weights, factors):   # factors already standardized 0-1
    suit += w_i * f
suit *= constraint_mask                 # binary 0/1, applied last
```

OWA variant: sort factor values per cell and apply order weights — full
AND (min) to full OR (max) continuum; use when stakeholders disagree on
risk tolerance and show 2-3 scenarios.

## Sensitivity analysis — mandatory

A result that flips with a small weight change is not a result:

- **One-at-a-time**: perturb each weight ±20% (renormalize), recompute,
  report % of area changing suitability class and a stability map (cells
  that never change class across perturbations).
- **Scenario**: 2-3 alternative weight sets from different stakeholder
  priorities; present side-by-side.
- If a Monte Carlo budget exists: sample weights from Dirichlet around the
  AHP vector; per-cell probability of "highly suitable" is a far stronger
  product than a single map.

## Deliverable standard

Suitability map (classified + continuous), constraint mask map, weights
table with CR, standardization functions per criterion (with direction),
sensitivity/stability summary, and limitations paragraph (data currency,
resolution, criteria omitted). Route cartography to `cartography-geoviz`;
network-access criteria come from `network-accessibility-analysis`.

## Pitfalls checklist

- Direction inversion on a criterion (the silent killer — double-check
  distance-based factors).
- Mixing resolutions without declaring the resampling rule.
- CR ignored or unreported.
- Constraints blended as weights → forbidden zones scored medium.
- Classifying with quantiles then reading them as absolute suitability.
- No sensitivity analysis; single map presented as truth.

## Execution contract

- **Workflow:** define decision and stakeholders; separate constraints from factors; standardize criteria; elicit and validate weights; aggregate; test sensitivity; communicate uncertainty.
- **Decision rules:** use MCDA for transparent criteria-ranked surfaces, network analysis for route-constrained access, and optimization when discrete placement or capacity decisions dominate.
- **Verification protocol:** check criterion direction and alignment, AHP consistency, constraint enforcement, weight and threshold perturbations, and stable-versus-fragile areas.
- **Failure modes:** reject the model when criteria double-count the same construct, weights lack provenance, constraints leak into compensation, or rankings collapse under plausible perturbations.
- **Deliverables:** continuous and classified suitability maps, constraints, criteria transformations, weights and consistency ratio, sensitivity results, and limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before applying methods or implementation APIs and record the checked date.

Referenced files: 3

View saved version →

---
name: ml-experiment-standards
description: >-
  Always invoke for training, validating, tuning, benchmarking, or claiming
  readiness of a predictive model. Covers leakage audits, spatial and grouped
  splits, metrics, reproducibility, and honest reporting. Invoke especially
  when spatial dependence, split design, or deployment geography is unknown;
  uncertainty is a reason to use this skill. Do not trigger for descriptive
  EDA or non-predictive statistical inference.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# ML Experiment Standards

Purpose: every ML job (quick prototypes included) is reproducible,
leakage-free, and metric-justified. These are not optional polish; every
skipped item typically returns as "the model collapsed in production" or
"the result didn't replicate".

## 1. EDA comes first

Before any model, produce and show: distributions, missingness rates,
outliers, target balance, salient correlations. Metric and loss choice
depend on this information; a model recommendation without EDA is a guess.

## 2. Leakage audit

At every split decision, answer explicitly (and write the answer as a code
comment): "Does the training set contain indirect information about any
test sample?"

| Data type | Correct split | Why |
|---|---|---|
| Independent samples | Stratified k-fold | Preserves class ratios |
| Time series | TimeSeriesSplit / walk-forward | Future must not leak into past |
| **Spatial data** | Spatial block CV — see `references/spatial-cv-protocol.md` | Neighbors are near-duplicates |
| Grouped data (patient, parcel, scene) | GroupKFold | A group must not straddle the split |

- Scalers/encoders/imputers are **fit on train only**; the clean path is
  `sklearn.pipeline.Pipeline` — CV then fits correctly by construction.
- Target-derived features (target encoding etc.) must be computed
  out-of-fold, and shown to be.

The spatial protocol in `references/spatial-cv-protocol.md` is the single
canonical source for this repo — other skills link here; do not restate it.

## 3. Metric selection — justified

Never choose a metric by default; write a one-sentence rationale:

- Imbalanced classes → **F1 / AUC-PR**, not accuracy (accuracy rewards
  majority-class memorization).
- Segmentation → **IoU/Dice** (pixel accuracy is inflated by background).
- Regression → RMSE (sensitive to large errors) vs MAE (robust) vs R²
  (variance explained) — justify from the use case.
- Every point estimate gets uncertainty: bootstrap CI or mean ± std across
  CV folds. A single number hides whether a difference is signal or noise.

## 4. Reproducibility skeleton

Every training script follows this shape (script-first; no notebook magic):

```python
"""Experiment: <name>. Goal and success criterion: <one sentence>."""
from dataclasses import dataclass, asdict
import json, random
import numpy as np

@dataclass
class Config:
    seed: int = 42
    lr: float = 1e-3
    batch_size: int = 32
    epochs: int = 100
    patience: int = 10  # early stopping

def set_seed(seed: int) -> None:
    random.seed(seed)
    np.random.seed(seed)
    # if torch: torch.manual_seed(seed); torch.cuda.manual_seed_all(seed)

def main(cfg: Config) -> None:
    set_seed(cfg.seed)
    ...  # data -> split -> pipeline -> train -> evaluate
    with open("runs/run_meta.json", "w", encoding="utf-8") as f:
        json.dump({"config": asdict(cfg), "metrics": metrics}, f, indent=2)

if __name__ == "__main__":
    main(Config())
```

- Config lives in a dataclass/YAML, never hardcoded — sweeps and run
  comparison depend on it.
- Pin library versions (`pip freeze > requirements.txt`).
- Use MLflow/W&B when available; the JSON log above is the minimum.

## 5. Deep learning extras

- **Loss rationale**: Dice/Dice+CE for imbalanced segmentation; write why.
  Focal only after comparison — not a free win.
- **Augmentation rationale**: state which transforms respect the physics
  of the problem (orientation-dependent tasks forbid some rotations;
  multispectral forbids naive color jitter).
- **Overfitting control**: early stopping with patience + a train/val
  curve in the report; no curve, no "the model is good".
- **Capacity order**: small model + simple baseline first (logistic
  regression, RF); a deep model that can't beat the baseline is a data
  problem, not an architecture problem.
- EO-specific chipping/inference details → `geo-deep-learning`.

## 6. System context (MLOps)

Position every model in its chain in one paragraph: data source →
cleaning → features/versioning → training → evaluation → deployment
(batch/real-time) → monitoring (data/model drift). Even for a prototype,
note "what this step becomes in production".

## 7. Report format

```
## Experiment: <name>
- Data: n=<>, split: <strategy + rationale>
- Baseline: <model> → <metric ± CI>
- Model: <model> → <metric ± CI>
- Leakage audit: <what was checked>
- Next step: <single recommendation>
```

When reporting differences, respect statistical honesty: if the gap
doesn't exceed the across-fold std, say "no clear difference" — no
p-hacking, no selective reporting.

## Execution contract

- **Workflow:** define prediction target and decision use; establish a baseline; audit leakage; create spatially valid splits; train reproducibly; quantify uncertainty; inspect errors and deployment fit.
- **Decision rules:** apply this skill only to predictive model experiments; use spatial statistics for inference, geostatistics for sampled-surface estimation, and descriptive analysis without forcing a model.
- **Verification protocol:** reproduce from a clean environment, compare against baseline across folds or seeds, inspect spatial residuals, verify split independence, and test the final decision threshold.
- **Failure modes:** invalidate uplift claims for leakage, post-split preprocessing, inappropriate metrics, non-independent test units, selective runs, or train-serving skew.
- **Deliverables:** experiment configuration, split and seed manifest, baseline and model metrics with uncertainty, leakage audit, error analysis, artifacts, and deployment caveats.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before using version-sensitive split, metric, or reproducibility APIs.

Referenced files: 3

View saved version →

---
name: movement-trajectory
description: >-
  Movement and trajectory analytics from GPS/GNSS tracks: cleaning, stop/trip
  detection, road-network map matching, speed/direction, flow aggregation,
  and origin-destination construction. Use for fleets, human mobility, animal
  tracking, AIS, or sports tracks. Trigger on GPS points, GPX, trajectories,
  stop detection, map matching, or timestamped positions per moving object.
  Also invoke for privacy, aggregation, de-identification, or release of
  individual trajectories. Use network-accessibility-analysis for
  hypothetical routes, isochrones, or static OD costs without observed tracks.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Movement & Trajectory Analytics

Purpose: turn noisy timestamped points into defensible movement facts. The
recurring failure modes: **speed computed through GPS noise** (teleporting
points → 400 km/h pedestrians), **stops invented by signal drift**, and
**privacy-blind delivery** of individual-level traces.

## Data model first

A trajectory = ordered fixes per object: `(object_id, timestamp, x, y,
[accuracy, ...])`. Before analysis, report per object: fix count, time
span, median sampling interval, and interval distribution — **sampling
rate drives every method choice** (1 s vehicle traces and 1 fix/hour
animal tags are different problems wearing the same schema).

```python
import movingpandas as mpd
import geopandas as gpd

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat),
                       crs=4326).to_crs(gdf_utm_epsg)
tc = mpd.TrajectoryCollection(gdf, "object_id", t="timestamp")
```

### Declare CRS and time base before any threshold

Every distance radius, speed limit and dwell duration in this skill is
meaningless until two things are stated **in the answer, before the number is
used**:

1. **The projected CRS** all distance and speed computation runs in — a "50 m
   stop radius" applied to raw lon/lat degrees is not 50 m anywhere, and the
   error scales with latitude. Name the CRS (`estimate_utm_crs()` for a local
   fleet, an equal-distance projection for continental extents).
2. **The timestamp base**, normalised to timezone-aware UTC. Fleet logs mix
   local times, DST shifts and naive strings; a dwell that straddles a DST
   boundary gains or loses an hour, and stop durations silently corrupt.

State both before proposing a radius or duration, not afterwards as a caveat.

**Declaring is not withholding.** An unknown CRS or timezone is never grounds to
stop and ask instead of answering. State it as an explicit, named assumption and
deliver the method anyway:

> Assuming a local UTM zone for distance and that timestamps are naive local time
> needing UTC normalisation — confirm both, since they change dwell durations.

Then give the cleaning steps, the parameters, and the sensitivity check. A
response that asks for the CRS, the timezone, or the file *in place of* the
method has failed this skill even if the question is a good one. Ask alongside
the answer, never instead of it. Never silently treat naive timestamps as UTC —
but "silently" is the operative word: an assumption you have labelled and
surfaced is exactly what is wanted.

## Cleaning pipeline (in order)

1. **Deduplicate** identical (object, timestamp) fixes.
2. **Accuracy filter**: drop fixes above an HDOP/accuracy threshold if
   the column exists (report the threshold and % dropped).
3. **Speed filter**: drop fixes implying impossible speed for the mode
   (walk > 15 km/h sustained, car > 200 km/h...); iterate — one bad fix
   creates two bad segments (`mpd.OutlierCleaner`).
4. **Gap splitting**: split trajectories at temporal gaps (e.g., > 5×
   median interval) — interpolating across a tunnel/power-off invents
   movement.
5. Optional smoothing (Kalman/rolling median) for jittery urban-canyon
   data — AFTER outlier removal, and never before stop detection tuning.

Accounting line per step: fixes in → out.

## Stops and trips

Stop = spatial dwell: fixes within a distance radius for a minimum
duration (`mpd.TrajectoryStopDetector(max_diameter=50,
min_duration=timedelta(minutes=5))`). The two parameters ARE the result —
report them and run a ±50% sensitivity check; urban-canyon drift mimics
movement, so diameter < GPS noise level yields zero stops.

Trips = segments between stops. Deliver per trip: origin, destination,
start/end time, duration, length, main mode guess if applicable. OD
matrices: aggregate trip endpoints to zones (see privacy below);
accessibility questions on the resulting flows → `network-accessibility-analysis`.

## Map matching

Raw GPS does not sit on the road. For any road-referenced claim (distance
driven, street-level flows, speeding), match to the network first:

- Tools: Valhalla (Meili), OSRM `match`, or `mappymatch`; HMM-based
  matchers are the standard.
- Sampling interval > ~30 s degrades matching sharply — report match
  confidence and the % of unmatched points; don't silently keep unmatched
  geometry.
- Never map-match animal tracks or off-road movement (obviously) — and
  don't compute "distance traveled" from raw noisy fixes either
  (noise inflates path length ~5-20%); smooth first, state the method.

## Aggregate analytics

- **Flow maps / desire lines**: aggregate OD pairs before plotting
  (`cartography-geoviz` for delivery); hairball avoidance = zone-level
  aggregation + minimum-flow threshold.
- **Density**: KDE or hex-bin of fixes vs of trips — fixes overweight slow
  movement (dwell = many fixes); use trip-based or time-weighted density
  and say which.
- **Space-time clustering** (co-location, convoys): ST-DBSCAN family;
  cluster parameters in both space and time reported together.
- Sequence/periodicity: hour-of-day × day-of-week activity matrices per
  object class before any behavioral claims.

## Privacy — non-optional

Individual trajectories are personal data almost everywhere (GDPR etc.)
and are notoriously re-identifiable (home/work anchor pairs identify most
people). Defaults: aggregate before sharing (zones ≥ k objects,
suppress cells < k, typical k=5-10), truncate trip ends near homes,
and never publish raw individual traces without explicit clearance.
State the anonymization applied in every deliverable.

## Verification protocol

1. Speed histogram per mode after cleaning — tail must be physically
   plausible.
2. Map 3 sample trajectories (raw vs cleaned vs matched) over a basemap.
3. Stop-detection sensitivity: parameters ±50%, report stop-count change.
4. OD totals reconcile with trip counts (accounting).

## Pitfalls checklist

- Speeds computed across gaps or through outlier fixes.
- Distance traveled from raw (unsmoothed, unmatched) fixes.
- Stops detected with radius below GPS noise, or drift counted as trips.
- Mixed timezones / DST jumps creating phantom teleports.
- Fix-density maps read as movement-density maps.
- Individual traces shipped without aggregation/suppression.
- Trajectories split by object but not by temporal gap.

## Execution contract

- **Workflow:** validate identifiers, time, and CRS; segment tracks; remove impossible fixes; infer stops or trips; optionally map-match; aggregate; apply privacy controls; verify.
- **Decision rules:** use trajectory methods for observed timestamped movement, network analysis for possible routes or access, and point-pattern methods when sequence and identity are absent.
- **Verification protocol:** inspect speed and gap distributions, map raw-versus-cleaned samples, perturb stop parameters, reconcile trip and OD counts, and audit disclosure risk.
- **Failure modes:** suppress or qualify results for timezone ambiguity, long gaps, implausible speeds, poor network matching, sparse sampling, or re-identification risk.
- **Deliverables:** cleaned trajectories or approved aggregates, segmentation rules, quality report, derived stop/trip tables, privacy treatment, maps, and limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before using format, library, or privacy guidance and record the checked date.

Referenced files: 2

View saved version →

---
name: network-accessibility-analysis
description: >-
  Always invoke for access to facilities or opportunities by walking,
  driving, cycling, or public transport, even for a conceptual question with
  no routing terms or data yet. Covers hospital and service access,
  transit/GTFS, routes, isochrones, OD matrices, closest facility, 2SFCA,
  walkability, coverage, and equity. Invoke when Euclidean buffers proxy for
  network access. Use movement-trajectory for observed tracks and MCDA for
  suitability without network costs.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Network & Accessibility Analysis

Purpose: replace as-the-crow-flies guesswork with network-true travel
costs, at the right scale and with honest assumptions about speeds and
modes. First decision on every task: Euclidean distance is only acceptable
as a declared approximation — flag it whenever you see it standing in for
access.

## Tool selection by scale

| Scale | Tool |
|---|---|
| Neighborhood-city, research flexibility | **OSMnx + NetworkX** |
| City-region, many-to-many OD (>10⁴×10⁴) | **r5py** (multimodal + transit w/ GTFS) or **pandana** (contraction-hierarchy speed) |
| Production routing service | Valhalla / OSRM / OpenRouteService API |
| Proprietary stacks | ArcGIS Network Analyst (script it headlessly) |

NetworkX chokes on metro-scale many-to-many — don't loop `shortest_path`
over thousands of origins; switch tools instead.

## Graph construction (OSMnx)

```python
import osmnx as ox

G = ox.graph_from_place("City, Country", network_type="drive")  # walk/bike/all
G = ox.add_edge_speeds(G)          # imputes from highway tags where maxspeed missing
G = ox.add_edge_travel_times(G)    # edge attr: travel_time (s)
G = ox.project_graph(G)            # metric CRS before any distance work
```

- **network_type matters**: pedestrian analysis on a `drive` graph misses
  paths, stairs, plazas; driving on `all` uses footpaths. Match mode.
- Imputed speeds are averages by road class — a systematic bias, not
  noise. State it; calibrate against known trips when stakes are high.
- Keep the strongly connected component for routing
  (`ox.truncate.largest_component(G, strongly=True)`); orphan islands
  cause spurious infinities.
- **Snapping**: origins/destinations map to nearest nodes/edges
  (`ox.distance.nearest_nodes`). Report the snap-distance distribution;
  a facility snapped 2 km away (riverside, gated area) silently corrupts
  results.

## Core products

- **Isochrones / service areas**: ego-graph by travel_time cutoff → alpha
  shape or buffered edge union around reached edges. Node-based convex
  hulls overstate coverage across rivers/highways — prefer edge-based
  polygons. Always label the assumptions: mode, speed model, cutoff.
- **OD matrix**: many-to-many travel costs; the substrate for
  accessibility and location-allocation. For big matrices use
  pandana/r5py; store as Parquet with origin/destination IDs.
- **Closest facility**: k-nearest by network cost (not Euclidean); report
  both the assigned facility and the cost.
- **Centrality**: betweenness on travel_time (sampled `k` for big
  graphs — exact is O(nm)); edge betweenness ≈ through-traffic potential.
  Interpret as network structure, not observed traffic.

## Accessibility metrics — pick deliberately

| Metric | Question it answers | Weakness |
|---|---|---|
| Cumulative opportunities (# jobs/POIs within T min) | Simple, communicable | Cliff at T; all-or-nothing |
| Gravity-based (distance-decayed sum) | Smooth access | Decay parameter must be justified |
| **2SFCA / E2SFCA** | Supply-demand ratio access (health care standard) | Catchment size choice drives results |
| Closest-facility time | Worst-case need | Ignores capacity/congestion |

For equity analyses, join metrics to population/demographic polygons
(area-weighted or dasymetric — see `geo-data-engineering`) and report
distributions per group, not just city means. Route statistical testing of
disparities to `spatial-statistics`.

## Location-allocation

Optimal siting (p-median, max-coverage) on the OD matrix: formulate with
PuLP/OR-Tools; inputs are the OD matrix + demand weights + candidate
sites. State the objective explicitly — minimize mean travel time
(p-median) vs maximize covered demand within T (max-coverage) give
different answers, and stakeholders rarely know which they asked for.
Feed results back to `mcda-suitability-analysis` when siting mixes network
access with other criteria.

## Transit (GTFS)

Use r5py with OSM + GTFS feeds; results are departure-time sensitive —
compute over a time window (e.g., 07:00-09:00 percentiles), never a single
departure. Validate the feed (calendar coverage on your analysis date!) —
an expired GTFS calendar yields walking-only times that look plausible.

## Verification protocol

1. Spot-check 3 routes against an external router (Google/OSRM) — within
   ~20% or explain why.
2. Map unreachable/infinite-cost pairs — usually snapping or connectivity
   artifacts, not real inaccessibility.
3. Isochrone eyeball: does it respect rivers, highways, one-ways?

## Pitfalls checklist

- Euclidean buffers presented as "service areas".
- Wrong network_type for the mode.
- Convex-hull isochrones bridging barriers.
- Snap distances unchecked.
- One departure time for transit accessibility.
- Betweenness sold as traffic volume.
- OD matrix in degrees-CRS travel "distances".

## Execution contract

- **Workflow:** define mode, time, impedance, origins, destinations, and equity question; build and validate the network; snap inputs; compute routes or matrices; summarize access; verify.
- **Decision rules:** use network costs for constrained travel, movement analytics for observed tracks, and MCDA only when accessibility becomes one criterion in a broader preference model.
- **Verification protocol:** audit connectivity and snapping, spot-check routes, map unreachable pairs, test departure-time or impedance sensitivity, and reconcile OD dimensions and units.
- **Failure modes:** withhold access claims for disconnected graphs, wrong mode or turn rules, expired GTFS service, excessive snapping, Euclidean substitution, or unstable departure-time results.
- **Deliverables:** network provenance, assumptions and cost function, routes or OD matrix, isochrones or access metrics, unreachable-case report, validation evidence, and equity caveats.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before using network, GTFS, or routing APIs and archive source dates.

Referenced files: 2

View saved version →

---
name: point-cloud-lidar
description: >-
  LiDAR and point cloud processing: PDAL pipelines, LAS/LAZ/COPC handling,
  ground classification, DTM/DSM/CHM generation, canopy and building
  metrics, and photogrammetric (SfM) point clouds. Use when the primary input
  is LAS, LAZ, COPC, LiDAR, or an unstructured 3D point cloud. This skill owns
  vertical datum agreement, co-registration and the vertical-accuracy budget
  when two acquisitions are differenced with comparability not yet established;
  once datum, geoid and accuracy are documented, a subsidence or
  elevation-change question is change-detection's. Route a derived DEM, DTM,
  DSM or CHM to terrain-hydrology unless point-level classification,
  comparability, or metrics remain in scope.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Point Clouds & LiDAR

Purpose: from raw returns to defensible elevation and structure products.
The recurring failure modes: **trusting vendor classification blindly**,
**mixing return types in surfaces** (DSM from last returns, DTM with
vegetation), and **ignoring point density** when choosing output
resolution.

## First contact with any cloud

```bash
pdal info input.laz --summary        # counts, bounds, CRS, classes, returns
```

Report before touching anything: point count, density (pts/m² — decides
achievable raster resolution), CRS (horizontal AND vertical datum —
ellipsoidal vs orthometric heights differ by the geoid undulation, tens of
meters in places), classification present?, return numbers present?,
flight-line overlap artifacts. A cloud without CRS metadata: resolve from
the provider, never assume.

## Format and scale

| Format | Use |
|---|---|
| **LAZ** | Compressed interchange/archive — default |
| **COPC** (cloud-optimized LAZ) | Streaming/HTTP range access, web viewers |
| LAS | Only when a tool can't read LAZ |
| Entwine/EPT | Massive multi-tile collections, indexed |

Tile large collections; process per-tile with buffered edges (~2× search
radius) to avoid seam artifacts in filters and surfaces; drop the buffer
on write.

## PDAL pipeline pattern

```json
{
  "pipeline": [
    "input.laz",
    {"type": "filters.reprojection", "out_srs": "EPSG:32636"},
    {"type": "filters.outlier", "method": "statistical",
     "mean_k": 8, "multiplier": 2.5},
    {"type": "filters.smrf", "slope": 0.15, "window": 18.0,
     "threshold": 0.5, "scalar": 1.2},
    {"type": "writers.las", "filename": "classified.laz",
     "extra_dims": "all"}
  ]
}
```

Run: `pdal pipeline pipeline.json`. Denoise BEFORE ground classification
(low outliers below ground destroy SMRF/CSF); tune `slope` up for steep
terrain, `window` to the largest non-ground object (big buildings need
bigger windows).

## Ground classification & DTM

- If vendor class 2 (ground) exists: **audit it** on 2-3 cross-sections
  (bridges, dense canopy, steep slopes) before trusting; reclassify where
  it fails.
- Algorithms: SMRF (PDAL default, robust), CSF (cloth simulation, good in
  steep forest). Parameters are terrain-dependent — show a cross-section
  plot as evidence, not just the parameter list.
- DTM from ground-only points; interpolation: TIN → raster (standard for
  DTM) or IDW for dense clouds. Output resolution ≥ ~1/√density; a 0.5 m
  DTM from 1 pt/m² data is invented detail.
- DSM from **first returns / highest-point** binning. CHM = DSM − DTM,
  clamp negatives to 0, and use a pit-free algorithm for forestry (naive
  CHMs are pocked by within-crown pits).

## Structure metrics

- **Forestry**: height percentiles (p95 ≈ canopy height), canopy cover
  (first returns > 2 m / all first returns), density metrics per grid cell
  or plot; normalize heights against the DTM first (`filters.hag_dem` or
  `filters.hag_nn`). Individual tree detection: local maxima on pit-free
  CHM + watershed segmentation — validate count against field plots or
  manual photo-interpretation samples.
- **Buildings**: class 6 or planar-patch extraction; building height =
  p90(roof points HAG); footprint fusion with cadastre/OSM polygons via
  zonal statistics on HAG.
- Downstream terrain analysis (slope, watersheds) → `terrain-hydrology`;
  DL on point clouds or derived rasters → `geo-deep-learning`.

## SfM/photogrammetric clouds — not LiDAR

Drone photogrammetry clouds have no returns, no canopy penetration
(ground under vegetation is guessed), correlated noise, and possible doming
from poor camera calibration. A "DTM" from SfM over forest is a canopy
model. State the sensor type in every deliverable; use LiDAR-specific
claims (penetration, return metrics) only for LiDAR.

## Vertical datum: resolve, transform, record

Every elevation product carries three obligations, and the third is the one
that gets skipped. **Stating the datum in your answer is not recording it.**
A height product whose vertical datum lives only in a chat reply is
indistinguishable from one with no datum at all the moment the file is
handed to anyone else.

1. **Resolve.** Read the vertical CRS from the header/VLR. Where it is
   missing or contradicted, resolve it against acquisition metadata (vendor
   flight report, project spec) or diagnose it: a tile-wide constant offset
   matching the local geoid undulation is the ellipsoidal-vs-orthometric
   fingerprint. Never infer a datum from elevation magnitude alone.
2. **Transform.** Apply an explicit, named transformation — a compound CRS
   plus geoid model through `filters.reprojection`, or a fitted per-tile
   offset when no geoid grid is available. Re-difference the overlaps
   afterwards and confirm strips agree within noise.
3. **Record it into the output, not just the reply.** Every delivered
   product must carry, in machine-readable form:
   - the compound or vertical CRS written into the file itself (LAS/LAZ
     header VLR, GeoTIFF CRS, or PROJJSON in the sidecar);
   - the **geoid model name and version** actually applied (e.g. EGM2008,
     GEOID18) and the transformation pipeline or EPSG operation code;
   - the **per-tile offsets applied**, where correction was per tile, with
     the control or reference each was fitted against;
   - the source of truth used to resolve an originally missing datum;
   - the residual strip-edge disagreement after correction.

   Emit this as a sidecar (`*.prj`/PROJJSON, a metadata JSON, or embedded
   raster tags) alongside the product, and never publish a height product
   whose vertical datum is unresolved. If the datum cannot be resolved,
   deliver the product labelled provisional with the unresolved datum
   recorded in the same metadata block — silence is not an option.

## Verification protocol

1. Cross-sections (2-3, including a building edge and a vegetated slope):
   ground class hugs terrain, DSM caps surface.
2. DTM minus known control points / national DEM: report RMSE and check
   for a constant offset = vertical datum mismatch.
3. Hillshade the DTM — classification artifacts (pits, pimples,
   flight-line stripes) are instantly visible.
4. Report density, CRS + vertical datum, classifier + parameters, and
   output resolution rationale in the answer — **and** confirm the verified
   vertical datum, geoid model, and transformation were written into the
   output metadata before the product is considered delivered.

## Pitfalls checklist

- Ellipsoidal heights delivered as orthometric (whole product offset by
  the geoid).
- Vertical datum resolved during the audit but never written into the
  delivered product's metadata — the next consumer inherits the same
  ambiguity you just spent the analysis removing.
- DTM resolution finer than point density supports.
- Vendor ground class trusted under dense canopy.
- CHM with negative values or crown pits (no pit-free processing).
- Per-tile processing without buffers → seam lines in derivatives.
- Outlier filter run AFTER ground classification.
- SfM cloud treated as canopy-penetrating LiDAR.

## Execution contract

- **Workflow:** inspect header, CRS, vertical datum, density, classes, and returns; tile with buffers; filter noise; classify; derive products; mosaic; validate in 3D and cross-section.
- **Decision rules:** use point-cloud workflows when return-level 3D evidence matters, terrain workflows after a validated DEM exists, and separate assumptions for LiDAR versus SfM clouds.
- **Verification protocol:** reconcile point counts and classes, inspect buffered seams and cross-sections, compare elevations to control, hillshade derived terrain, and report density-supported resolution.
- **Failure modes:** stop for unknown vertical datum, insufficient density, corrupt classification, tile seams, unbounded outliers, or product resolution finer than sampling supports. Never resolve a vertical datum and then ship the product without that datum and its transformation recorded in the output metadata.
- **Deliverables:** validated cloud or derived DTM/DSM/CHM, pipeline parameters, CRS and vertical datum, density and class report, QA graphics, accuracy metrics, and limitations. The verified vertical datum, the geoid model and transformation applied, and any per-tile offsets are written into the output metadata or a sidecar, not only into the answer text.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before applying format, quality, or processing rules and record the checked date.

Referenced files: 2

View saved version →

---
name: postgis-spatial-sql
description: >-
  Invoke whenever spatial SQL or its execution backend is the decision:
  PostGIS, DuckDB Spatial, SpatiaLite, ST_* functions, recurring spatial
  joins, concurrent/growing workloads, or large GeoParquet queries. Covers
  backend selection, schemas, GiST/BRIN indexes, KNN, geometry versus
  geography, correctness benchmarks, and EXPLAIN optimization. Use PostGIS
  for managed concurrent services and embedded engines for bounded local
  analytics when evidence supports that choice. Use geo-data-engineering for
  acquisition, conversion, and file-based ETL without spatial SQL.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# PostGIS & Spatial SQL

Purpose: correct-and-fast spatial SQL. The two recurring failure modes are
semantic (geometry vs geography, SRID mismatches → wrong answers) and
performance (missing index usage → hour-long joins); this skill guards
both.

## When the database is the right tool

Move from files/GeoPandas to PostGIS when any of: features > a few
million, concurrent readers/writers, repeated ad-hoc querying, a serving
API on top, or transactional integrity needs. For single-shot analytical
scans over GeoParquet, **DuckDB Spatial** is often the fastest
zero-install path — same SQL mindset, no server.

When requirements are incomplete, do not turn this heuristic into a final
recommendation. First obtain current and forecast data volume, concurrency,
delivery and mutation pattern, latency/SLA, serving needs, and operational
ownership (including backup and recovery). Define representative ingestion,
join, and read queries for both viable backends; compare runtime and resource
use only after row counts, join cardinality, SRID, geometry validity, and sample
outputs agree. Include this benchmark and correctness plan in the current
response; do not merely offer to draft it later.

## Schema fundamentals

This runnable example assumes the data is contained in UTM zone 33N. Replace
EPSG:32633 with a projected CRS verified for the actual area of interest.

```sql
CREATE TABLE parcels (
  id          bigint GENERATED ALWAYS AS IDENTITY PRIMARY KEY,
  parcel_no   text NOT NULL,
  landuse     text,
  area_m2     double precision,          -- unit in the name, always
  geom        geometry(MultiPolygon, 32633) NOT NULL
);
CREATE INDEX parcels_geom_gix ON parcels USING gist (geom);
ANALYZE parcels;
```

- **Type the geometry column fully**: `geometry(MultiPolygon, SRID)` — an
  untyped `geometry` column happily accepts mixed garbage.
- Promote to Multi* on load (`ST_Multi`) so Polygon/MultiPolygon mixing
  never bites.
- **geometry vs geography**: geometry in a projected SRID for regional
  analysis (fast, full function set); geography (SRID 4326) when the
  extent is global/cross-zone and you want meters without picking a
  projection (slower, smaller function set). Never store in 4326 geometry
  and call `ST_Area` expecting m² — that's square degrees.
- Never use EPSG:3857/Web Mercator for area or length measurement. When the
  analysis CRS is not yet known, either use 4326 geography for a geodesic
  result or stop and select a verified local/equal-area CRS; do not present a
  known-distorting CRS as a runnable measurement alternative.
- **Any stored geometry column you recommend must be typed with its SRID.**
  Advising a "second projected geometry column" for repeated measurement is
  incomplete until it is written as `geometry(<Type>, <SRID>)` with the index
  and the populating `ST_Transform`. An untyped column recommended as a fix
  reintroduces the mixed-SRID problem it was meant to solve:

  ```sql
  ALTER TABLE parcels ADD COLUMN geom_32633 geometry(MultiPolygon, 32633);
  UPDATE parcels SET geom_32633 = ST_Transform(geom, 32633);
  CREATE INDEX parcels_geom_32633_gix ON parcels USING gist (geom_32633);
  ```
- GiST index on every geometry column, `ANALYZE` after bulk loads; BRIN
  only for huge, spatially-ordered, append-only tables.
- Load paths: `ogr2ogr -f PostgreSQL`, `shp2pgsql`, or GeoPandas
  `to_postgis` (small/medium). `COPY` beats INSERT by orders of magnitude.

## Correct spatial predicates

- `ST_Intersects` for "touches at all", `ST_Contains`/`ST_Within` for
  containment, `ST_DWithin(a, b, dist)` for proximity — **never**
  `ST_Distance(a,b) < dist` (that form can't use the index).
- The classic point-in-polygon join:

```sql
SELECT p.id, a.district
FROM points p
JOIN admin a ON ST_Intersects(a.geom, p.geom);   -- GiST on both sides
```

- KNN nearest-neighbor with the distance operator (index-assisted):

```sql
SELECT h.id, h.name
FROM hospitals h
ORDER BY h.geom <-> (SELECT geom FROM incident WHERE id = 42)
LIMIT 3;
```

`<->` gives true-distance ordering on modern PostGIS for geometry; wrap
with `ST_DWithin` to bound the search when tables are huge.

## Performance playbook

1. `EXPLAIN (ANALYZE, BUFFERS)` first — confirm the GiST index is used
   (look for "Index Scan ... _gix"); a Seq Scan on a big spatial join
   means a rewrite, not a bigger server.
2. Same SRID on both sides of every predicate — `ST_Transform` inside a
   join predicate kills index use; store a transformed, indexed copy
   instead.
3. Big-polygon problem: country/basin-sized geometries make index bboxes
   useless → `ST_Subdivide` into a work table (typical 10-100× speedup on
   joins against them).

The following example assumes `countries(country_id, geom)`.

```sql
CREATE TABLE country_parts AS
SELECT c.country_id, part.geom
FROM countries AS c
CROSS JOIN LATERAL ST_Subdivide(c.geom, 256) AS part(geom);

CREATE INDEX country_parts_geom_gix ON country_parts USING gist (geom);
ANALYZE country_parts;
```

`ST_Subdivide` is a set-returning function; do not access its result as
`(ST_Subdivide(...)).geom`.

4. Validity in-database: `ST_IsValid` audit, `ST_MakeValid` repair, add a
   `CHECK (ST_IsValid(geom))` if writers are untrusted.
5. Simplify for serving, not for analysis: keep full-resolution geometry;
   generate `ST_SimplifyPreserveTopology` copies or vector tiles
   (`ST_AsMVT`) for the web tier.
6. Batch updates in transactions; `VACUUM ANALYZE` after churn.

## Common analytical patterns

```sql
-- Area-weighted aggregation (e.g., population into custom zones)
SELECT z.zone_id,
       SUM(b.pop * ST_Area(ST_Intersection(z.geom, b.geom)) / ST_Area(b.geom)) AS pop_est
FROM zones z JOIN blocks b ON ST_Intersects(z.geom, b.geom)
GROUP BY z.zone_id;

-- Dissolve with attribute
SELECT landuse, ST_Multi(ST_Union(geom))::geometry(MultiPolygon, 32633) AS geom
FROM parcels GROUP BY landuse;
```

Area-weighted interpolation assumes uniform density within source units —
state that assumption when reporting. Validity repair is `ST_MakeValid`,
never `ST_Buffer(geom, 0)`.

## DuckDB Spatial quick path

```sql
INSTALL spatial; LOAD spatial;
SELECT a.name, count(*)
FROM 'admin.parquet' a, 'points.parquet' p
WHERE ST_Intersects(a.geom, p.geom)
GROUP BY a.name;
```

Reads GeoParquet/Shapefile/GPKG directly, parallel by default — ideal for
one-off large joins and pipeline steps without a server. No GiST; it plans
its own joins — benchmark, don't assume.

## Verification protocol

1. Row-count accounting query after each join/overlay CTE.
2. `SELECT DISTINCT ST_SRID(geom), GeometryType(geom)` on every table
   touched — one query kills two classic bug families.
3. Sample 5 output features rendered over a basemap (QGIS connects
   directly) — numbers can pass while geometries are garbage.
4. Treat every `sql` fence presented as runnable as a syntax and alias
   boundary: it must execute top-to-bottom after stated schema assumptions.
   Never put angle-bracket placeholders, ellipses, pseudocode, abandoned joins,
   or incomplete aliases inside it. If a schema value such as an SRID is
   unknown, ask for it or keep the template in a labeled `text` block.

## Pitfalls checklist

- `ST_Area`/`ST_Length` on 4326 geometry (square degrees).
- EPSG:3857/Web Mercator for area or length measurement (systematic distortion).
- `ST_Distance < x` instead of `ST_DWithin` (no index).
- `ST_Transform` in join predicates.
- Untyped geometry columns with mixed SRIDs.
- Country-sized polygons joined without `ST_Subdivide`.
- `buffer(0)` as validity repair (silent part loss) — `ST_MakeValid`.
- Serving full-resolution geometries to web clients.

## Execution contract

- **Workflow:** inspect schema, SRID, geometry type, size, and query goal; choose predicates and indexes; write auditable CTEs; inspect the plan; reconcile results; operationalize safely.
- **Decision rules:** use PostGIS for concurrent, repeated, or transactional spatial workloads; use file pipelines or DuckDB Spatial for bounded one-off transformations when a server adds no value.
- **Verification protocol:** assert SRID and geometry invariants, account for rows at each join, compare indexed plans and timings, sample geometries on a map, and test boundary semantics.
- **Failure modes:** block release for mixed SRIDs, accidental many-to-many explosion, invalid geometries, non-indexable predicates, geography/geometry unit confusion, or unexplained plan regressions.
- **Deliverables:** self-contained parameterized SQL or migration with consistent CTE/table aliases, indexes and rationale, query plan evidence, row accounting, sample validation, expected schema, performance notes, and rollback guidance.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) for the deployed database and extension versions before selecting functions or plans.

Referenced files: 2

View saved version →

---
name: remote-sensing-analysis
description: >-
  Always invoke for classical analysis, classification, validation, or
  comparability of satellite, aerial, or drone imagery. This skill owns
  sensor, product, processing-level and processing-baseline harmonization,
  including multi-date inputs; add change-detection only after comparable
  observations exist. Two scenes of the same product level are not
  automatically comparable: Sentinel-2 L2A crossed a reflectance offset at
  Processing Baseline 04.00 in January 2022, so any pair spanning that date
  starts here. Covers spectral indices, masking, compositing, SAR, land cover,
  and accuracy assessment. Route neural methods to geo-deep-learning and
  planetary server-side execution to google-earth-engine.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Remote Sensing Analysis

Purpose: turn raw Earth observation imagery into defensible analytical
products. The failure modes here are subtle — uncorrected DNs treated as
reflectance, clouds counted as land cover change, indices computed on the
wrong bands — so this skill front-loads the checks.

## Data access (STAC-first)

Search via STAC APIs rather than per-provider portals; the workflow is
uniform and scriptable:

```python
import pystac_client
import odc.stac

catalog = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[27.0, 38.3, 27.4, 38.6],
    datetime="2025-05-01/2025-09-30",
    query={"eo:cloud_cover": {"lt": 20}},
).item_collection()
ds = odc.stac.load(items, bands=["red", "nir", "scl"], resolution=10, chunks={})
```

Key collections: `sentinel-2-l2a` (10 m optical, surface reflectance),
`landsat-c2-l2` (30 m, 1982→), `sentinel-1-grd` (SAR, weather-independent).
Microsoft Planetary Computer mirrors most (needs `planetary_computer`
signing). For continental/global extents or decades-long stacks, route to
`google-earth-engine` instead of downloading. Record collection + item IDs +
search parameters for reproducibility.

## Processing-level discipline

| Level | Meaning | Analysis-ready? |
|---|---|---|
| L1C / L1TP | Top-of-atmosphere (TOA) | Indices OK-ish; cross-date comparison risky |
| **L2A / L2SP** | Surface reflectance (BOA) | Yes — default choice |
| GRD (SAR) | Detected amplitude | Needs terrain correction + speckle filter |

Always state which level you used. Never mix TOA and BOA scenes in one
composite or time series. Landsat Collection 2 L2 needs its scale factors
applied (`reflectance = DN * 0.0000275 - 0.2`).

### The Sentinel-2 baseline discontinuity — passes the level check above

Processing Baseline 04.00, applied from **25 January 2022**, added a constant
`BOA_ADD_OFFSET` (currently −1000) to L2A digital numbers so that negative
surface reflectance can be encoded. Two scenes on opposite sides of that date
are **both L2A**: the level check above sees nothing wrong while their DNs sit
1000 apart. Differencing them yields a systematic reflectance shift that reads
as real change and survives every mask, threshold and accuracy report you
apply afterwards.

- Read `BOA_ADD_OFFSET` and `QUANTIFICATION_VALUE` from each product's
  metadata rather than hardcoding −1000 and 10000; both are per-band and the
  baseline has changed before.
- Convert with `reflectance = (DN + BOA_ADD_OFFSET) / QUANTIFICATION_VALUE`.
- Record the **processing baseline of every scene** in the manifest, not just
  the product level. Two L2A scenes is not a sufficient statement.
- **Do not correct twice.** Harmonised collections — Earth Engine's
  `COPERNICUS/S2_SR_HARMONIZED` and several commercial mirrors — have already
  shifted post-baseline data back to the pre-2022 range. Applying the offset
  again inverts the error rather than removing it.
- If the baseline is undocumented for either scene, the comparison is not
  defensible. Say that instead of assuming pre- or post-2022.

## Cloud and quality masking — before anything else

- Sentinel-2: mask with SCL band (drop classes 3 cloud shadow, 8-9 clouds,
  10 cirrus, 11 snow — keep 4 vegetation, 5 bare, 6 water, 7 unclassified
  with care).
- Landsat C2: decode `QA_PIXEL` bitfields (cloud, shadow, cirrus bits).
- Report the % of valid pixels after masking per scene; scenes below ~60%
  valid usually deserve exclusion.
- For gap-free products, build median composites over a season rather than
  cherry-picking single scenes.

## Spectral indices

Compute on surface reflectance, guard against division by zero, and name
bands explicitly — band **numbers differ across sensors** (NIR is B8 on
Sentinel-2, B5 on Landsat 8/9):

```python
import numpy as np
import xarray as xr

def normalized_diff(a: xr.DataArray, b: xr.DataArray) -> xr.DataArray:
    """(a - b) / (a + b) with zero-denominator protection."""
    return xr.where(a + b == 0, np.nan, (a - b) / (a + b))

ndvi = normalized_diff(ds.nir, ds.red)     # vegetation
ndwi = normalized_diff(ds.green, ds.nir)   # open water (McFeeters)
ndbi = normalized_diff(ds.swir16, ds.nir)  # built-up
```

Interpretation guardrails: NDVI thresholds are scene- and season-dependent;
never hardcode "NDVI > 0.3 = vegetation" without checking the histogram.
Water confuses NDBI; shadows mimic water in NDWI — cross-check indices
against each other and against true-color.

## Classification workflow

1. Define a legend with mutually exclusive, imagery-separable classes.
2. Collect training samples spatially spread across the scene; record them
   as a versioned vector file.
3. Features: bands + indices + texture (GLCM) + temporal statistics if
   multi-date. For deep learning routes, hand off to `geo-deep-learning`.
4. Validate with a **spatially independent** test set (see
   `ml-experiment-standards` → `references/spatial-cv-protocol.md`) and
   report per-class F1/IoU plus a confusion matrix — overall accuracy alone
   hides rare-class failure.
5. Map the errors: a spatial plot of misclassifications reveals systematic
   problems (terrain shadow, urban/bare confusion) that global metrics hide.

## SAR notes (Sentinel-1)

Preprocess: orbit file → thermal noise removal → calibration (σ⁰) →
terrain correction (Range-Doppler with a DEM) → speckle filter (Lee/Refined
Lee) → dB conversion. Work in dB for statistics; VV/VH ratio is a strong
water/vegetation discriminator. SAR sees through clouds — prefer it for
flood mapping and continuous monitoring.

## Pitfalls checklist

- Comparing scenes across dates without consistent atmospheric correction.
- Mixing Sentinel-2 scenes across the 2022-01-25 baseline change without
  applying `BOA_ADD_OFFSET` — or applying it a second time on a collection
  that is already harmonised.
- Ignoring 20 m→10 m band mixing on Sentinel-2 (B11/B12 are natively 20 m).
- Computing indices on integer DNs without scale factors → nonsense ranges.
- Median composites of SAR in linear units (do statistics in dB).
- Training and test pixels from the same field/polygon → leaked accuracy.
- Forgetting nodata masks after reprojection (edges become zeros → fake
  land cover).

## Execution contract

- **Workflow:** define phenomenon and scale; select sensor, product level, and dates; harmonize calibration, masks, CRS, and resolution; derive features; analyze; validate spatially; publish provenance.
- **Decision rules:** use this skill for imagery preparation and classical analysis, change detection for explicit temporal differencing, deep learning for neural training, and Earth Engine for archive-scale execution.
- **Verification protocol:** inspect masks and valid counts, confirm scale factors, offsets and processing baseline per scene, confirm band resolution, overlay outputs, use spatially independent validation, map errors, and test seasonal or sensor sensitivity.
- **Failure modes:** reject results for cloud or shadow leakage, incomparable processing levels, undocumented or mismatched processing baselines, resampling artifacts, label leakage, nodata contamination, or claims beyond sensor resolution.
- **Deliverables:** analysis-ready imagery or features, processing manifest, masks, derived products, validation metrics and error map, reproducible code, and limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) at execution time for product, calibration, and catalog changes.

Referenced files: 2

View saved version →

---
name: spatial-statistics
description: >-
  Always invoke before testing a geographic pattern for clustering, hotspots,
  dependence, or explanatory regression, even when aggregation or ordinary
  OLS is proposed as routine. Covers Moran's I, LISA, Getis-Ord Gi*, weights,
  MAUP and scale sensitivity for areas/grids, residual dependence, and
  spatial lag/error/GWR/MGWR models. Use ML standards for predictive
  evaluation and geostatistics for continuous surfaces from sparse samples.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Spatial Statistics

Purpose: answer "is it clustered, where, and why" with defensible inference.
The core discipline: spatial data violates independence assumptions, so
standard statistics silently overstate significance — every analysis here
starts with weights design and ends with residual diagnostics.

## Spatial weights (W) — the analysis IS the weights

Every result downstream depends on W; choose it for substantive reasons and
run a sensitivity check with one alternative:

| Weights | Use when |
|---|---|
| Queen/Rook contiguity | Irregular polygons (admin units, parcels) |
| K-nearest neighbors | Points; islands present (contiguity leaves them unconnected) |
| Distance band | Physical process with known range |
| Kernel (distance-decayed) | Smooth influence, GWR-style local models |

```python
from libpysal.weights import Queen

w = Queen.from_dataframe(gdf, use_index=True)
print(f"islands: {w.islands}")   # unconnected units break stats — fix or document
w.transform = "r"                # row-standardize (default for Moran/lag models)
```

Always report: weights type, parameters, number of islands, and whether
results survive an alternative W.

## Global → local workflow

1. **Global Moran's I** (`esda.Moran`, permutation inference ≥999) —
   answers "any clustering at all?" Report I, p_sim, and the permutation
   distribution, not the analytical p.
2. **LISA / local Moran** (`esda.Moran_Local`) — maps WHERE: High-High,
   Low-Low clusters, High-Low/Low-High outliers. Correct for multiple
   testing (FDR at minimum) before coloring a map — uncorrected LISA maps
   overstate clusters and this is the field's most common abuse.
3. **Getis-Ord Gi\*** (`esda.G_Local`, star=True) — hot/cold spots of
   intensity (a distinct question from Moran clusters — Gi* finds
   concentrations of high values, LISA finds similarity structure).
4. Rates, not counts, for population-based phenomena; use Empirical Bayes
   smoothing (`esda.smoothing`) for small-population units before any of
   the above — raw rates in sparse units are noise.

## Point patterns

- Separate first-order intensity (density varies) from second-order
  interaction (points attract/repel) — KDE describes the former, Ripley's
  K/L (`pointpats`) tests the latter.
- Always test against an inhomogeneous null when the study area has obvious
  density gradients (population, roads); CSR against a city is a strawman.
- KDE bandwidth drives the story: report it, justify it (Silverman/CV), and
  show one alternative.

## Spatial regression decision path

Run OLS first, then diagnose — never start with a spatial model:

```python
from spreg import OLS
ols = OLS(y, X, w=w, spat_diag=True, moran=True, name_y="price", name_x=xnames)
```

Decision (Anselin's rule via LM tests): LM-Lag significant & LM-Error not →
**spatial lag (SAR)**; reverse → **spatial error (SEM)**; both → compare
robust LM versions; neither → OLS stands (report that as a finding).
Interpretation caveats: in SAR, coefficients are NOT marginal effects —
report direct/indirect (spillover) effects. In SEM, spatial structure is
nuisance correlation, no spillover story allowed.

**GWR/MGWR** (`mgwr`): when relationships plausibly vary over space.
Bandwidth by AICc search; map local coefficients WITH local t-values masked
for insignificance; MGWR when predictors operate at different scales.
GWR is exploratory — resist causal language on local coefficients.

## Inference honesty

- Permutation p-values over analytical ones wherever available.
- Multiple testing: n local tests = n units; FDR-correct.
- MAUP (modifiable areal unit problem): results can flip with unit
  aggregation — if the aggregation level is a choice, test one alternative
  and disclose.
- Spatial autocorrelation in residuals after modeling = model still wrong;
  report residual Moran's I for every final model.
- Correlation ≠ causation applies doubly here: spatially confounded
  variables (everything correlates with "distance to coast") demand
  explicit identification strategies before causal claims.

## Reporting template

```
## Spatial analysis: <question>
- Units & n, variable(s), rate smoothing: <...>
- W: <type/params>, islands: <n>, sensitivity W: <type>
- Global: Moran's I = <> (p_perm = <>)
- Local: <k> significant clusters after FDR; map attached
- Model: <OLS/SAR/SEM/GWR> chosen because <LM diagnostics>
- Residual Moran's I: <> — <interpretation>
- Caveats: MAUP, W-sensitivity, causal limits
```

## Execution contract

- **Workflow:** define inferential question and unit; inspect distributions and rates; construct and justify spatial weights; run global before local tests; fit models if needed; diagnose residual dependence; report uncertainty.
- **Decision rules:** use spatial statistics for dependence and inference, geostatistics for interpolating sampled continuous surfaces, and predictive ML when out-of-sample prediction is the primary goal.
- **Verification protocol:** test alternative weights and aggregation, use valid permutation or model inference, correct local multiplicity, inspect residual Moran's I, and distinguish association from causation.
- **Failure modes:** withhold inferential claims for arbitrary weights, islands ignored, unstable MAUP results, uncorrected multiple tests, residual autocorrelation, or unsupported causal language.
- **Deliverables:** analysis-ready variables, weights specification, global and local results, corrected significance, diagnostic maps, model and residual checks, sensitivity analysis, and caveats.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before applying version-sensitive statistical APIs or defaults.

Referenced files: 2

View saved version →

---
name: swe-devops-standards
description: >-
  Always invoke to review, repair, or deliver geospatial or GeoAI code,
  including contract compliance, security, error handling, transactions,
  tests, scripts, functions, notebooks, packages, CI/CD, and repository
  changes, even when deployment is not requested. Pair with the domain skill
  for ETL and other production code. Covers CRS/data invariants, dependencies,
  cross-platform reproducibility, automation, and shipping. Do not trigger for
  unrelated software or analysis requesting no code or repository artifact.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Geospatial SWE & DevOps Standards

Purpose: code produced as part of geospatial work should run in the user's
real environment and meet peer-level engineering quality. Apply these rules
only when code or repository artifacts are in scope.

## 1. Environment realities (the top error source)

- **Script-first by default**: no `%matplotlib inline`, `!pip install`, or
  `display()` unless the user is explicitly in a notebook. Every file runs
  from a terminal via `python script.py` behind an
  `if __name__ == "__main__":` block. (Cell markers like `# %%` are fine
  as an addition — the script must also work without them.)
- **Cross-platform paths**: always `pathlib.Path`; never string-concatenate
  or hardcode `/` or `\\`. Ask or detect the user's OS before giving shell
  commands; give CMD/PowerShell syntax on Windows, POSIX elsewhere —
  don't mix (`export` vs `set`, `venv/bin/activate` vs
  `venv\Scripts\activate`).
- **Encodings**: explicit `encoding="utf-8"` on every text file open —
  Windows still defaults to legacy code pages, and non-ASCII content
  corrupts silently.
- **Modern Python (3.11+)**: `X | None` unions, `type` aliases, structural
  pattern matching where they clarify; state the minimum version if a
  feature requires it.

## 2. Code quality defaults

Applied to every generated function/module, even when not asked:

```python
def compute_share(values: list[float], total: float) -> list[float]:
    """Return each value's share of the total.

    Args:
        values: Values to compute shares for.
        total: Denominator; must be non-zero.

    Returns:
        Shares in the same order as values.

    Raises:
        ValueError: If total is zero.
    """
    if total == 0:
        raise ValueError("total must be non-zero — share is undefined.")
    return [v / total for v in values]
```

- Type hints on every signature; `dataclass`/`TypeAlias` for complex types.
- Google-style docstrings; one-liners suffice for trivial functions.
- **Never bare `except:`**; catch specific exceptions, handle or re-raise
  with `raise ... from e`. A silent `pass` costs a week of debugging.
- `logging` over `print` (leveled, formatted), except user-facing CLI
  output.
- Note algorithmic complexity where it matters ("this is O(n log n), safe
  at n>10⁶") — especially around nested loops and pandas `apply`.
- Magic numbers → named module-level constants.

## 3. Testing and verification

- Offer at least a skeleton pytest for every function carrying real logic:

```python
# test_compute.py — run: python -m pytest -q
import pytest
from compute import compute_share

def test_basic() -> None:
    assert compute_share([1, 1], 2) == [0.5, 0.5]

def test_zero_total_raises() -> None:
    with pytest.raises(ValueError):
        compute_share([1.0], 0)
```

- Numerical code: test edge cases — empty input, NaN, negatives, single
  element.
- Run generated code yourself when an execution environment exists;
  otherwise mark it explicitly "not executed" — no silent assumptions.

## 4. Dependencies and reproducibility

- New project → virtual environment + pinned `requirements.txt`
  (`package==version`); never "install the latest".
- Seed randomness and put the seed in config (details in
  `ml-experiment-standards`).
- Note environment-difference risks where relevant (BLAS, CUDA, locale).

## 5. Git practices

- Conventional Commits: `feat(scope): ...`, `fix: ...`, `refactor: ...`;
  the body explains *why* — the diff already shows *what*.
- Commit in meaningful units; warn against 500-line single commits.
- Default `.gitignore`: `venv/`, `__pycache__/`, `*.pyc`, large data files
  (suggest DVC/LFS), IDE folders.

## 6. Automation / DevOps

- **CI**: minimal GitHub Actions for test + lint (ruff); note OS-runner
  differences if jobs must run on Windows too.
- **Docker**: start from `python:3.12-slim`, simple single-stage until
  size/caching demands more; note image size and build-cache implications.
- **Monitoring**: any long-lived service/pipeline ships three signals
  minimum: structured logs, failure alerting, basic metrics (duration,
  volume). ML services add drift checks (see `ml-experiment-standards`).
- **Scheduled jobs**: match the user's platform — cron on POSIX,
  Task Scheduler (`schtasks`) on Windows.

## 7. Code review mode

Review in this order and report findings by severity: correctness (edge
cases, silent failures) → security (injection, secrets, path traversal) →
performance (N+1, needless copies, O(n²)) → readability. Every finding
ships with the suggested fix as code — never "this is bad" and nothing
else.

## Execution contract

- **Workflow:** clarify the geospatial code's contract; reproduce the environment; inspect correctness and data invariants; implement the smallest safe change; test; package; document operations and rollback.
- **Decision rules:** apply this skill to geospatial software and pipeline delivery, not generic non-spatial coding; scale CI, containers, and observability to the actual deployment risk.
- **Verification protocol:** run focused and regression tests, lint and type checks where configured, exercise CRS/nodata/geometry edge cases, verify clean installation, and review CI artifacts.
- **Failure modes:** block release for silent data loss, nondeterminism, mutable hidden state, unpinned critical dependencies, secrets, platform assumptions, missing rollback, or unhandled spatial edge cases.
- **Deliverables:** reviewed code, tests, reproducible environment and lock data, CI configuration, operational notes, risk-ranked findings, observability plan, and rollback instructions.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before applying packaging, CI, testing, or supply-chain guidance.

Referenced files: 2

View saved version →

---
name: terrain-hydrology
description: >-
  Always invoke for terrain, drainage, viewshed, or visibility analysis from
  elevation, even before the DEM or correct surface is chosen. Covers
  DTM-versus-DSM selection, slope, aspect, curvature, hillshade, conditioning,
  flow direction/accumulation, streams, watersheds, and catchments. Use
  point-cloud-lidar first only when an elevation surface must be created from
  LiDAR or photogrammetric points.
license: MIT
metadata:
  author: Muhammed Enes Duran
---

# Terrain & Hydrology

Purpose: terrain products whose numbers are physically meaningful. The two
recurring failure modes: **unit mismatch** (degree coordinates with meter
elevations silently corrupts every derivative) and **unconditioned DEMs**
(flow routed into spurious pits produces fragmented, fictional streams).

## DEM hygiene first

| Check | Rule |
|---|---|
| Surface type | **DTM** (bare earth) for hydrology/slope; **DSM** (with canopy/buildings) for viewshed/solar. Using a DSM for watersheds routes rivers over treetops. |
| Source | Copernicus GLO-30 > SRTM for most global work; national LiDAR DTMs when available (see `point-cloud-lidar` to make your own). Record source + acquisition date. |
| Nodata | Identify the nodata value (-9999, -32768, 3.4e38) and mask it — never let it enter statistics or fill algorithms as "very deep hole". |
| Voids | Fill data voids (interpolation from edges) BEFORE hydrological conditioning; document filled areas. |
| **CRS + units** | Reproject to a projected CRS so horizontal units = vertical units (meters). Slope from a 4326 DEM without z-factor correction is the classic silent error. If staying geographic, apply a latitude-dependent z-factor — better: don't. |

## Derivatives

```python
import whitebox
wbt = whitebox.WhiteboxTools()
wbt.slope("dem.tif", "slope_deg.tif", units="degrees")
wbt.aspect("dem.tif", "aspect_deg.tif")
wbt.plan_curvature("dem.tif", "plan_curv.tif")
```

- Slope: state units (degrees vs percent — 45° = 100%); Horn's method
  (3×3) is the standard; steeper terrain → consider resolution effects
  (slope flattens as cell size grows — report cell size with every slope
  statistic).
- Aspect: circular variable — never average it arithmetically; use vector
  (sin/cos) averaging; flat cells have undefined aspect (mask, don't zero).
- Curvature: plan (flow convergence) vs profile (flow acceleration) —
  pick per question.
- Hillshade is for cartography (see `cartography-geoviz`), never analysis
  input.
- Ruggedness/position: TRI, TPI (radius-dependent — report the radius),
  geomorphons for landform classification.

## Hydrological conditioning — order matters

```
voids filled → breach depressions (preferred) → fill remaining pits
→ flow direction → flow accumulation → streams → watersheds
```

- **Breaching before filling** (WhiteboxTools
  `BreachDepressionsLeastCost`): carves through barriers (road embankments
  over culverts) instead of flooding upstream areas flat. Pure fill on
  flat/embanked terrain creates large artificial lakes with arbitrary flow
  paths.
- Real depressions exist (karst, prairie potholes, reservoirs). If the
  landscape genuinely holds water, don't condition it away — model with
  explicit sink handling and say so.
- Flow direction: **D8** for stream networks/watersheds (discrete,
  standard); **D-infinity/MFD** for dispersal quantities (wetness index,
  erosion) on hillslopes.

## Streams and watersheds

- Stream extraction threshold (min. accumulation) is a MODELING choice:
  derive from a mapped reference network (match total stream length) or
  report the threshold and show two alternatives — never present one
  threshold's network as "the" rivers.
- **Pour point snapping**: outlet coordinates rarely fall on the modeled
  stream cell. Snap to the highest-accumulation cell within a search
  radius (`wbt.jenson_snap_pour_points`) — an unsnapped pour point yields
  a tiny, wrong watershed silently.
- Verify delineation: watershed area vs authoritative basin data (±5-10%),
  and the modeled network overlaid on imagery/topo maps at 3 locations.
- Wetness index (TWI), stream power (SPA): compute from conditioned DEM +
  MFD accumulation; they are relative indices — don't read absolute
  thresholds across regions.

## Viewshed

- Use a **DSM** (or DTM + feature heights) — bare-earth viewsheds
  overstate visibility wherever trees/buildings exist; state which surface
  was used.
- Set observer height (~1.7 m person, tower height for infrastructure) and
  target height explicitly; defaults differ across tools.
- Account for earth curvature + refraction beyond ~5 km
  (`wbt.viewshed` handles it; verify the flag).
- Deliver binary visible/not plus the observer point(s) and parameters in
  the metadata; for siting problems, cumulative viewsheds from candidate
  sets feed `mcda-suitability-analysis`.

## Tooling

WhiteboxTools (conditioning, full hydrology suite, fast) · `pysheds`
(lightweight Python watersheds) · `richdem` (derivatives) · GDAL
(`gdaldem`) for quick slope/hillshade · GRASS (`r.watershed`) for very
large DEMs (no explicit fill needed — least-cost routing).

## Verification protocol

1. Derivative histograms: slope > 60° over large areas or negative
   accumulation = unit/nodata bug.
2. Stream network overlay on imagery at 3 locations, including one flat
   area (where artifacts concentrate).
3. Watershed area cross-check vs authoritative basin polygons.
4. Report: DEM source/date/resolution, conditioning method, flow
   algorithm, stream threshold, all in the deliverable.

## Pitfalls checklist

- Slope from a geographic-CRS DEM without z-factor (values ~100× off).
- DSM used for watershed delineation (rivers over treetops).
- Fill-only conditioning across road embankments → phantom lakes.
- Unsnapped pour point → 3-cell "watershed".
- Arithmetic mean of aspect (350° and 10° average to south, not north).
- Nodata treated as elevation in fill/statistics.
- One arbitrary stream threshold presented as the drainage network.

## Execution contract

- **Workflow:** inspect DEM source, CRS, vertical units, datum, resolution, and nodata; condition terrain; derive gradients and flow; delineate products; test thresholds; validate against imagery and controls.
- **Decision rules:** use terrain workflows on raster elevation products, point-cloud workflows before DEM generation, and choose conditioning and flow algorithms from landscape and scale.
- **Verification protocol:** inspect derivative distributions, hillshade artifacts, stream overlays, watershed area, pour-point snapping, threshold sensitivity, and elevation-control residuals.
- **Failure modes:** reject products from DSM misuse, geographic-unit slope, vertical datum mismatch, unconditioned barriers, nodata contamination, unsnapped outlets, or resolution unsupported by source data.
- **Deliverables:** conditioned DEM, derivatives and hydrologic products, parameter and threshold record, CRS and vertical datum, QA maps, validation metrics, and limitations.
- **Source freshness:** consult [the authoritative source registry](references/authoritative-sources.md) before applying tool algorithms or product rules and record the checked date.

Referenced files: 2