fix(tactical): implement scientific D8 hydrological flow routing for natural wadis and isolate from ridges

2026-08-18 12:44:58 +03:00
parent c9b3249aa5
commit 13ff137359
2 changed files with 242 additions and 58 deletions
@@ -585,15 +585,17 @@ export const TacticalDefenseView: React.FC = () => {
            ['get', 'layerType'],
            'cliff', 3.5,
            'road', 4.0,
-            'wadi', 3.0,
+            'wadi', 3.5,
+            'ridge', 2.5,
            2.5
          ],
          'line-dasharray': [
            'match',
            ['get', 'layerType'],
            'cliff', ['literal', [2, 1]],
+            'ridge', ['literal', [3, 1.5]],
            'road', ['literal', [1, 0]],
-            'wadi', ['literal', [2, 1.5]],
+            'wadi', ['literal', [1, 0]],
            ['literal', [1, 0]]
          ]
        }
@@ -835,9 +837,9 @@ export const TacticalDefenseView: React.FC = () => {
    }

    const filtered = terrainResult.spatialGeoJson.features.filter((f: any) => {
-      if (f.properties.category === 'peak' || f.properties.category === 'valley') return showPeaksAndValleys;
+      if (f.properties.category === 'peak' || f.properties.category === 'valley' || f.properties.layerType === 'ridge') return showPeaksAndValleys;
      if (f.properties.layerType === 'slope-sector') return showSlopeSectors;
-      if (f.properties.layerType === 'wadi' || f.properties.layerType === 'ridge' || f.properties.layerType === 'cliff') return showNaturalObstacles;
+      if (f.properties.layerType === 'wadi') return showNaturalObstacles;
      if (f.properties.layerType === 'road') return showRoadCorridors;
      if (f.properties.layerType === 'urban') return showUrbanZones;
      if (f.properties.layerType === 'hazard' || f.properties.layerType === 'cliff') return showHazardousGround;
@@ -733,8 +733,8 @@ export async function calculateTerrainStudy(
  const centerElev = await sampleElevationAt(centerLat, centerLng);
  const areaKm2 = Math.round(Math.PI * Math.pow(radiusMeters / 1000, 2) * 10) / 10;

-  // High-Density Sampling Matrix: 21x21 = 441 nodes (400 terrain cells)
-  const gridSize = 21;
+  // High-Density Sampling Matrix: 25x25 = 625 nodes (576 terrain cells)
+  const gridSize = 25;
  const R_earth = 6371000;
  const dLatTotal = (radiusMeters / R_earth) * (180 / Math.PI);
  const dLngTotal = (radiusMeters / (R_earth * Math.cos((centerLat * Math.PI) / 180))) * (180 / Math.PI);
@@ -963,72 +963,254 @@ export async function calculateTerrainStudy(
  const avgSlope = slopeValues.length > 0 ? Math.round((slopeValues.reduce((a, b) => a + b, 0) / slopeValues.length) * 10) / 10 : 2.5;
  const maxSlope = slopeValues.length > 0 ? Math.round(Math.max(...slopeValues) * 10) / 10 : 5.0;

-  // 6. Trace Drainage Defiles & Natural Wadis (Lowland Thalweg Corridors)
-  const valleySegments: [number, number][] = [];
+  // 6. Scientific Hydrological D8 Flow Routing & Natural Wadi (Thalweg) Extraction
+  // Compute steepest descent neighbor for every grid cell
+  interface FlowTarget {
+    r: number;
+    c: number;
+    slope: number;
+  }
+
+  const flowTarget: (FlowTarget | null)[][] = Array(gridSize).fill(null).map(() => Array(gridSize).fill(null));
+
  for (let r = 0; r < gridSize; r++) {
-    let minRowElev = 9999;
-    let minColIdx = 0;
    for (let c = 0; c < gridSize; c++) {
-      if (elevationGrid[r][c] < minRowElev) {
-        minRowElev = elevationGrid[r][c];
-        minColIdx = c;
+      const zCurr = elevationGrid[r][c];
+      let maxDropSlope = 0;
+      let bestTarget: FlowTarget | null = null;
+
+      for (let dr = -1; dr <= 1; dr++) {
+        for (let dc = -1; dc <= 1; dc++) {
+          if (dr === 0 && dc === 0) continue;
+          const nr = r + dr;
+          const nc = c + dc;
+          if (nr < 0 || nr >= gridSize || nc < 0 || nc >= gridSize) continue;
+
+          const zNeighbor = elevationGrid[nr][nc];
+          const distM = Math.hypot(dr * dyMeters, dc * dxMeters);
+          const drop = zCurr - zNeighbor;
+          const slope = drop / distM;
+
+          if (slope > maxDropSlope) {
+            maxDropSlope = slope;
+            bestTarget = { r: nr, c: nc, slope };
+          }
+        }
      }
-    }
-    const pt = gridCoords[r * gridSize + minColIdx];
-    const dist = calculateDistance(centerLat, centerLng, pt.lat, pt.lng);
-    if (dist <= radiusMeters * 1.02) {
-      valleySegments.push([pt.lng, pt.lat]);
+
+      flowTarget[r][c] = bestTarget;
    }
  }

-  if (valleySegments.length >= 3) {
-    spatialFeatures.push({
-      type: 'Feature',
-      properties: {
-        layerType: 'wadi',
-        category: 'natural-obstacle',
-        title: `🌊 بطن وادٍ ومصرف سيل طبيعي (${minElev}م)`,
-        color: '#06b6d4'
-      },
-      geometry: {
-        type: 'LineString',
-        coordinates: valleySegments
-      }
-    });
-  }
+  // Calculate Flow Accumulation (Upslope contributing area)
+  const flowAcc: number[][] = Array(gridSize).fill(1).map(() => Array(gridSize).fill(1));

-  // 7. Trace Dominant Mountain Ridges (Key Terrain Crest Lines)
-  const ridgeSegments: [number, number][] = [];
+  // Sort cells from highest elevation to lowest
+  const sortedCells: Array<{ r: number; c: number; elev: number }> = [];
  for (let r = 0; r < gridSize; r++) {
-    let maxRowElev = -9999;
-    let maxColIdx = 0;
    for (let c = 0; c < gridSize; c++) {
-      if (elevationGrid[r][c] > maxRowElev) {
-        maxRowElev = elevationGrid[r][c];
-        maxColIdx = c;
-      }
+      sortedCells.push({ r, c, elev: elevationGrid[r][c] });
    }
-    const pt = gridCoords[r * gridSize + maxColIdx];
-    const dist = calculateDistance(centerLat, centerLng, pt.lat, pt.lng);
-    if (dist <= radiusMeters * 1.02 && maxRowElev >= minElev + relief * 0.5) {
-      ridgeSegments.push([pt.lng, pt.lat]);
+  }
+  sortedCells.sort((a, b) => b.elev - a.elev);
+
+  for (const cell of sortedCells) {
+    const target = flowTarget[cell.r][cell.c];
+    if (target) {
+      flowAcc[target.r][target.c] += flowAcc[cell.r][cell.c];
    }
  }

-  if (ridgeSegments.length >= 3) {
-    spatialFeatures.push({
-      type: 'Feature',
-      properties: {
-        layerType: 'ridge',
-        category: 'key-terrain',
-        title: `⛰️ سلسلة تلال حاكمة (Dominant Crest - ${maxElev}م)`,
-        color: '#38bdf8'
-      },
-      geometry: {
-        type: 'LineString',
-        coordinates: ridgeSegments
+  // Extract Thalweg Valley & Stream Channels
+  // Threshold for channel initiation: cells where accumulated flow exceeds minimum threshold
+  const minAccThreshold = Math.max(5, Math.floor(gridSize * 0.3));
+  const visitedStream = Array(gridSize).fill(false).map(() => Array(gridSize).fill(false));
+
+  // Find all channel heads (cells where flowAcc >= threshold, but none of upstream neighbors >= threshold)
+  const channelHeads: Array<{ r: number; c: number; acc: number }> = [];
+  for (let r = 1; r < gridSize - 1; r++) {
+    for (let c = 1; c < gridSize - 1; c++) {
+      if (flowAcc[r][c] >= minAccThreshold) {
+        // Check if any upstream neighbor that flows to (r, c) already had flowAcc >= minAccThreshold
+        let hasUpstreamThreshold = false;
+        for (let dr = -1; dr <= 1; dr++) {
+          for (let dc = -1; dc <= 1; dc++) {
+            if (dr === 0 && dc === 0) continue;
+            const nr = r + dr;
+            const nc = c + dc;
+            if (nr >= 0 && nr < gridSize && nc >= 0 && nc < gridSize) {
+              const target = flowTarget[nr][nc];
+              if (target && target.r === r && target.c === c && flowAcc[nr][nc] >= minAccThreshold) {
+                hasUpstreamThreshold = true;
+                break;
+              }
+            }
+          }
+          if (hasUpstreamThreshold) break;
+        }
+
+        if (!hasUpstreamThreshold) {
+          channelHeads.push({ r, c, acc: flowAcc[r][c] });
+        }
      }
-    });
+    }
+  }
+
+  // Trace streams downstream from each channel head
+  for (const head of channelHeads) {
+    const rawPath: [number, number][] = [];
+    const streamElevations: number[] = [];
+    let currR = head.r;
+    let currC = head.c;
+    let maxStreamAcc = head.acc;
+    let steps = 0;
+
+    while (currR >= 0 && currR < gridSize && currC >= 0 && currC < gridSize && steps < gridSize * 2) {
+      steps++;
+      const pt = gridCoords[currR * gridSize + currC];
+      const distFromCenter = calculateDistance(centerLat, centerLng, pt.lat, pt.lng);
+      
+      if (distFromCenter <= radiusMeters * 1.05) {
+        rawPath.push([pt.lng, pt.lat]);
+        streamElevations.push(elevationGrid[currR][currC]);
+        if (flowAcc[currR][currC] > maxStreamAcc) {
+          maxStreamAcc = flowAcc[currR][currC];
+        }
+      }
+
+      visitedStream[currR][currC] = true;
+      const next = flowTarget[currR][currC];
+      if (!next) break; // Reached sink or local depression
+      if (visitedStream[next.r][next.c]) {
+        // Confluence with existing river channel
+        const confluencePt = gridCoords[next.r * gridSize + next.c];
+        rawPath.push([confluencePt.lng, confluencePt.lat]);
+        streamElevations.push(elevationGrid[next.r][next.c]);
+        break;
+      }
+      currR = next.r;
+      currC = next.c;
+    }
+
+    // Only keep streams that have at least 3 points and actual length
+    if (rawPath.length >= 3) {
+      // Smooth the stream line using moving average
+      const smoothedPath: [number, number][] = [];
+      for (let i = 0; i < rawPath.length; i++) {
+        if (i === 0 || i === rawPath.length - 1) {
+          smoothedPath.push(rawPath[i]);
+        } else {
+          const pPrev = rawPath[i - 1];
+          const pCurr = rawPath[i];
+          const pNext = rawPath[i + 1];
+          smoothedPath.push([
+            Number((pPrev[0] * 0.25 + pCurr[0] * 0.5 + pNext[0] * 0.25).toFixed(6)),
+            Number((pPrev[1] * 0.25 + pCurr[1] * 0.5 + pNext[1] * 0.25).toFixed(6))
+          ]);
+        }
+      }
+
+      const streamMinZ = Math.min(...streamElevations);
+      const streamMaxZ = Math.max(...streamElevations);
+
+      spatialFeatures.push({
+        type: 'Feature',
+        properties: {
+          layerType: 'wadi',
+          category: 'natural-obstacle',
+          title: `🌊 مجرى وادٍ ومصرف سيل طبيعي (${Math.round(streamMinZ)}م - ${Math.round(streamMaxZ)}م)`,
+          flowAcc: maxStreamAcc,
+          color: '#0284c7',
+          width: 3.5
+        },
+        geometry: {
+          type: 'LineString',
+          coordinates: smoothedPath
+        }
+      });
+    }
+  }
+
+  // 7. Trace Natural Ridge Watershed Divides (Dominant High Crests)
+  // Crest cells: high elevation, local peak curvature, divergent flow (low flow accumulation)
+  const ridgeCandidates: Array<{ r: number; c: number; elev: number; tpi: number }> = [];
+  for (let r = 1; r < gridSize - 1; r++) {
+    for (let c = 1; c < gridSize - 1; c++) {
+      const el = elevationGrid[r][c];
+      if (el >= minElev + relief * 0.45 && flowAcc[r][c] <= 2) {
+        // Compute 3x3 local mean
+        let sum = 0;
+        let count = 0;
+        for (let dr = -1; dr <= 1; dr++) {
+          for (let dc = -1; dc <= 1; dc++) {
+            sum += elevationGrid[r + dr][c + dc];
+            count++;
+          }
+        }
+        const localMean = sum / count;
+        const tpi = el - localMean;
+        if (tpi >= 1.5) {
+          ridgeCandidates.push({ r, c, elev: el, tpi });
+        }
+      }
+    }
+  }
+
+  // Group connected ridge candidates into crest lines
+  const visitedRidge = Array(gridSize).fill(false).map(() => Array(gridSize).fill(false));
+  for (const rc of ridgeCandidates) {
+    if (visitedRidge[rc.r][rc.c]) continue;
+    const ridgeLine: [number, number][] = [];
+    let currR = rc.r;
+    let currC = rc.c;
+
+    while (currR >= 0 && currR < gridSize && currC >= 0 && currC < gridSize) {
+      visitedRidge[currR][currC] = true;
+      const pt = gridCoords[currR * gridSize + currC];
+      const distFromCenter = calculateDistance(centerLat, centerLng, pt.lat, pt.lng);
+      if (distFromCenter <= radiusMeters * 1.05) {
+        ridgeLine.push([pt.lng, pt.lat]);
+      }
+
+      // Find best unvisited ridge neighbor
+      let bestNext: { r: number; c: number } | null = null;
+      let highestNeighborElev = -9999;
+      for (let dr = -1; dr <= 1; dr++) {
+        for (let dc = -1; dc <= 1; dc++) {
+          if (dr === 0 && dc === 0) continue;
+          const nr = currR + dr;
+          const nc = currC + dc;
+          if (nr >= 0 && nr < gridSize && nc >= 0 && nc < gridSize && !visitedRidge[nr][nc]) {
+            const isCand = ridgeCandidates.some(c => c.r === nr && c.c === nc);
+            if (isCand && elevationGrid[nr][nc] > highestNeighborElev) {
+              highestNeighborElev = elevationGrid[nr][nc];
+              bestNext = { r: nr, c: nc };
+            }
+          }
+        }
+      }
+
+      if (!bestNext) break;
+      currR = bestNext.r;
+      currC = bestNext.c;
+    }
+
+    if (ridgeLine.length >= 3) {
+      spatialFeatures.push({
+        type: 'Feature',
+        properties: {
+          layerType: 'ridge',
+          category: 'key-terrain',
+          title: `⛰️ سلسلة تلال وقمم حاكمة (خط تقسيم مياه)`,
+          color: '#d97706',
+          width: 2.5
+        },
+        geometry: {
+          type: 'LineString',
+          coordinates: ridgeLine
+        }
+      });
+    }
  }

  // 8. Regional Sector Intelligence Context