version: "1.0.1"

Stage 1: Mesh Generation — Skill Document

Purpose

Generate the hexagonal (or other polygonal) mesh with DEM-derived elevation values in the JSON format consumed by HexWatershed. This stage bridges raw geospatial data (GeoTIFF DEMs) to the model's cell-based representation.

The mesh JSON is the single most critical input to HexWatershed. Every subsequent computation depends on correct cell geometry, neighbor relationships, and elevation values.

Inputs

Outputs

Output JSON Structure

{
  "sMesh_type": "hexagon",
  "nCell": 12345,
  "dResolution": 5000.0,
  "aCells": [
    {
      "lCellID": 1,
      "dLongitude_center_degree": -76.5,
      "dLatitude_center_degree": 39.2,
      "dElevation_mean": 152.3,
      "dElevation_raw": 152.3,
      "dArea": 21650000.0,
      "nVertex": 6,
      "vVertex": [{"dLongitude_degree": -76.48, "dLatitude_degree": 39.22}, ...],
      "aNeighbor": [2, 3, 4, 5, 6, 7],
      "aNeighbor_distance": [5000.0, 5000.0, 5000.0, 5000.0, 5000.0, 5000.0]
    }
  ]
}

Procedure

1. Obtain DEM in GCS (EPSG:4326) with meters elevation:

``bash # Reproject if needed gdalwarp -t_srs EPSG:4326 -r bilinear input_utm.tif output_gcs.tif ``

1. Verify DEM CRS and units:

``bash gdalinfo dem.tif | grep -E "AUTHORITY|Unit|NoData" # Must show: EPSG:4326, degree units, elevation in meters ``

1. Generate mesh using PyFlowline (recommended):

``python from pyflowline.formats.convert_mesh import pyflowline_generate_mesh pyflowline_generate_mesh( sDEM='/data/dem.tif', sMesh_type='hexagon', dResolution=5000, sOutput='/data/mesh/' ) ``

Or use the KI mesh_converter tool: ``bash python mesh_converter.py --dem /data/dem.tif --resolution 5000 \ --mesh-type hexagon --output /data/mesh/hexagon_mesh_info.json ``

1. Validate mesh:

• Check all cells have valid elevation (not nodata)
• Check neighbor relationships are symmetric (if A neighbors B, B neighbors A)
• Check cell areas are physically reasonable for the resolution

Verification

# Check mesh JSON is valid and has cells
python3 -c "
import json
mesh = json.load(open('hexagon_mesh_info.json'))
print(f'Cells: {mesh[\"nCell\"]}')
cells = mesh['aCells']
elev = [c['dElevation_mean'] for c in cells]
print(f'Elevation range: {min(elev):.1f} to {max(elev):.1f} m')
nodata = [c for c in cells if c['dElevation_mean'] <= -9999]
print(f'Nodata cells: {len(nodata)}')"
 
# Check coordinate system (values should be in degree range)
python3 -c "
import json
mesh = json.load(open('hexagon_mesh_info.json'))
lons = [c['dLongitude_center_degree'] for c in mesh['aCells']]
lats = [c['dLatitude_center_degree'] for c in mesh['aCells']]
print(f'Lon range: {min(lons):.4f} to {max(lons):.4f}')
print(f'Lat range: {min(lats):.4f} to {max(lats):.4f}')
assert max(abs(l) for l in lons) <= 180, 'PROJECTED CRS DETECTED'
assert max(abs(l) for l in lats) <= 90, 'PROJECTED CRS DETECTED'"

Traps

TRAP: DEM in projected CRS (dt_006)

If the DEM is in UTM or State Plane (coordinates in meters, not degrees), the mesh will have coordinates like (500000, 4000000) instead of (-76.5, 39.2). The model will compute nonsensical distances because it expects geographic degrees for its internal spherical geometry calculations. Always check bounds before proceeding.

TRAP: DEM elevation in feet (dt_013)

USGS NED tiles may have elevation in feet. Since all HexWatershed internal calculations assume meters, slopes will be inflated by a factor of 3.28, and depression filling will behave incorrectly. Verify: gdalinfo dem.tif should show meters or check max values.

TRAP: Cell area units (dt_011)

Cell area (dArea) must be in m² (square meters). If a GIS tool exports area in km², the value will be 10⁶ times too small, causing drainage density and all area-based metrics to be wildly wrong. A 5km hex cell should have area ~21.65 × 10⁶ m², not 21.65.

TRAP: Neighbor distance units (dt_012)

aNeighbor_distance must be in meters. If distances are in km or degrees, slope calculations (Δz / distance) will be off by orders of magnitude.

TRAP: Mesh type mismatch (dt_005)

The sMesh_type in the config must match the structure of the mesh JSON. A hexagon mesh has 6 neighbors per cell; an MPAS mesh has variable connectivity. Using the wrong parser silently produces garbage topology.

Example

For a 5km hexagonal mesh over the Chesapeake Bay watershed:

# DEM: MERIT Hydro, already in EPSG:4326, meters
python mesh_converter.py \
    --dem /data/merit_hydro_chesapeake.tif \
    --resolution 5000 \
    --mesh-type hexagon \
    --output /data/mesh/hexagon_mesh_info.json
 
# Expected: ~2000-5000 cells depending on watershed size