From 13ff1373594eee147747cdf3cac279c2c8385b52 Mon Sep 17 00:00:00 2001 From: Hamza-Ayed Date: Tue, 18 Aug 2026 12:44:58 +0300 Subject: [PATCH] fix(tactical): implement scientific D8 hydrological flow routing for natural wadis and isolate from ridges --- apps/web/src/pages/TacticalDefenseView.tsx | 10 +- apps/web/src/utils/elevationService.ts | 290 +++++++++++++++++---- 2 files changed, 242 insertions(+), 58 deletions(-) diff --git a/apps/web/src/pages/TacticalDefenseView.tsx b/apps/web/src/pages/TacticalDefenseView.tsx index 8a9ac58..04385bf 100644 --- a/apps/web/src/pages/TacticalDefenseView.tsx +++ b/apps/web/src/pages/TacticalDefenseView.tsx @@ -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; diff --git a/apps/web/src/utils/elevationService.ts b/apps/web/src/utils/elevationService.ts index ab57ec2..3f7f4c7 100644 --- a/apps/web/src/utils/elevationService.ts +++ b/apps/web/src/utils/elevationService.ts @@ -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