Exclude mapped land and retain checked around-peninsula guidance

README.ko.md

@@ -2,7 +2,7 @@
 
 [English](README.md)
 
-세계 항로·수치 지형을 검사한 해안 경로·항구 접근 경로를 생성합니다. 전역 그래프로 방문 순서를 정하고, 수치 지형으로 지역·로컬 경로를 확인합니다. WWW-POSIM은 이 계획기를 지형 작업자와 항해 관리기에 연결합니다.
+세계 항로·해안 경로·항구 접근 경로를 생성합니다. 전역 그래프로 방문 순서를 정합니다. 수치 지형으로 수심을 검사하고, 육지 지도와 겹치는 경로를 제외합니다. WWW-POSIM은 이 계획기를 지형 작업자와 항해 관리기에 연결합니다.
 
 ## 설치
 
@@ -50,6 +50,17 @@ path = planner.route(start, goal)
 
 확인된 수역 안에서 출발점과 목표를 선택합니다. 자료가 없는 범위는 통과할 수 없습니다. 차량이 지나가는 전체 폭을 삼각형 충돌 지형과 비교합니다.
 
+`land.LandExcludedField(field, guard)`로 수치 지형에 육지 검사를 추가합니다.
+`land.guard_from_spec(spec)`는 캐시의 OSM 육지 자료를 읽습니다. 다른 로컬
+압축 파일을 쓰려면 `WWOS_COAST_OSM_ARCHIVE`를 지정합니다. WWW-POSIM은 지형
+자료 캐시에 이 파일을 준비합니다. 자료가 없으면 지역 경로 준비를 멈춥니다.
+육지 자료로 경로를 제외할 수 있으며, 수심은 수치 지형으로 확인합니다.
+
+지역 경로는 반도를 돌아가는 전체 지도 경로를 찾은 뒤 수심을 검사한 로컬 목표를
+선택합니다. 우회 초반에는 항로 목표까지의 거리가 늘어날 수 있습니다. 검증된
+구간의 끝에서는 다음 구간이 준비될 때까지 기다립니다. 날짜 변경선에서도
+경도가 이어지도록 처리합니다.
+
 | 기준 | 기본값 |
 |---|---|
 | 수상 차량 로컬 수심 여유 | 1.5 m |
@@ -85,6 +96,7 @@ path = planner.route(start, goal)
 | `connections.py` | 해안을 고려한 항구 연결 |
 | `navigation.py`, `planning.py` | 충돌 지형·차량 폭을 고려한 A*·지형 캐시 |
 | `coastal.py` | 해안 탐색과 수치 지형 경로 |
+| `land.py`, `land_guidance.py` | 육지 경로 제외와 범위를 제한한 우회 안내 |
 | `regional.py` | 비동기 지역 경로 준비 |
 | `arrival_targets.py` | 수치 지형과 선택적인 OSM 육지 검사로 도착 수역 선택 |
 | `port_approach.py` | 도착 영역 접근과 곡선 감속 |

README.md

@@ -2,7 +2,7 @@
 
 [한국어](README.ko.md)
 
-Plan global maritime routes, terrain-checked coastal corridors and port approaches. A global graph supplies the itinerary; numeric terrain approves regional and local paths. WWW-POSIM connects these planners to its terrain workers and voyage coordinator.
+Plan global maritime routes, coastal corridors and port approaches. A global graph supplies the itinerary. Numeric terrain checks depth, and mapped land excludes shore crossings. WWW-POSIM connects these planners to its terrain workers and voyage coordinator.
 
 ## Install
 
@@ -50,6 +50,17 @@ path = planner.route(start, goal)
 
 Choose start and goal inside confirmed water. The planner rejects unknown coverage and continuously checks the swept footprint against the piecewise planar collision field.
 
+`land.LandExcludedField(field, guard)` adds mapped-land exclusion to a numeric
+field. `land.guard_from_spec(spec)` reads the cached OSM land polygon archive.
+Set `WWOS_COAST_OSM_ARCHIVE` to another local archive path. WWW-POSIM prepares
+this dataset in its terrain provider cache. Missing data holds regional planning.
+OSM land can reject a path; numeric terrain must still establish water depth.
+
+Regional lookahead finds a complete mapped detour around a peninsula before
+selecting a local depth-checked goal. The detour can initially increase distance
+to the course target. At its checked endpoint, the caller waits for the next
+checked segment. Longitude wrapping supports the international date line.
+
 Coastal lookahead follows a checked water route to the onward target or a forward
 exit from the regional grid. Its endpoint has an onward route around bays and
 headlands. Ordinary coastal course vertices guide progress; ports and required
@@ -86,6 +97,7 @@ search retains the supplied route. Numeric terrain checks remain necessary.
 | `connections.py` | Coast-aware port connectors |
 | `navigation.py`, `planning.py` | Collision field, swept-footprint A* and cached terrain |
 | `coastal.py` | Coastal discovery and numeric corridors |
+| `land.py`, `land_guidance.py` | Mapped-land exclusion and bounded around-land guidance |
 | `regional.py` | Asynchronous regional lookahead |
 | `arrival_targets.py` | Buffered numeric arrival water and optional OSM land checks |
 | `port_approach.py` | Arrival-region corridors and bend braking |

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "posim-routing"
-version = "0.1.2"
+version = "0.1.3"
 description = "Global maritime and terrain checked coastal route planning"
 requires-python = ">=3.10"
 license = {text = "Apache-2.0"}

src/posim_routing/__init__.py

@@ -1,2 +1,2 @@
 """POSIM-ROUTING: reusable ocean simulation components."""
-__version__ = "0.1.1"
+__version__ = "0.1.3"

src/posim_routing/coastal.py

@@ -1,7 +1,8 @@
-"""Prefer visible coast without treating a cartographic shoreline as collision data.
+"""Prefer visible coast with numeric depth and optional mapped-land exclusions.
 
 Natural Earth only discovers regions worth querying. Numeric terrain determines
-the shoreline, depth, and every regional corridor. The native Gazebo planner
+depth and every regional corridor. Mapped land can only reject travel.
+The native Gazebo planner
 independently approves each local corridor against its resident collision mesh.
 """
 import gzip
@@ -14,7 +15,7 @@
 
 import numpy as np
 from pyproj import Geod
-from scipy.ndimage import distance_transform_edt, maximum_filter
+from scipy.ndimage import distance_transform_edt, maximum_filter, label
 from scipy.spatial import cKDTree
 
 GEO = Geod(ellps="WGS84")
@@ -110,7 +111,7 @@ def coastal_approach(start, target, shore):
     return candidate
 
 
-def coastal_corridor(terrain, spec, start, target, timeout=25.):
+def coastal_corridor(terrain, spec, start, target, timeout=25., guard=None):
     """Return a progressive coastal path, or None to retain the maritime course.
 
     A bounded weighted A* prefers a 250 m stand-off. It cannot invent land or
@@ -128,6 +129,8 @@ def coastal_corridor(terrain, spec, start, target, timeout=25.):
     xs, ys, heights = field.xs, field.ys, field.heights
     finite = np.isfinite(heights)
     land = finite & (heights >= 0)
+    if guard is not None:
+        land |= ~guard.grid(xs, ys)
     if not land.any():
         return None
     origin = np.array(forward.transform(*start))
@@ -145,6 +148,8 @@ def coastal_corridor(terrain, spec, start, target, timeout=25.):
     bounds = maximum_filter(np.where(finite, heights, math.inf), size=size,
                             mode="constant", cval=math.inf)
     safe = bounds <= -MINIMUM_DEPTH_M
+    if guard is not None:
+        safe &= guard.grid(xs, ys, margin)
     coast_distance = distance_transform_edt(~land, sampling=(dy, dx))
     xx, yy = np.meshgrid(xs, ys)
     travel = np.hypot(xx-origin[0], yy-origin[1])
@@ -159,12 +164,20 @@ def coastal_corridor(terrain, spec, start, target, timeout=25.):
     # a through route to the onward target or a forward water exit, then take
     # its coastal prefix. Distance gain at an isolated endpoint is insufficient.
     goal_node = (int(np.argmin(abs(ys-destination[1]))), int(np.argmin(abs(xs-destination[0]))))
+    start_node = (int(np.argmin(abs(ys-origin[1]))), int(np.argmin(abs(xs-origin[0]))))
+    if not safe[start_node]:
+        return None
+    # Rank reachable exits. A nearer exit beyond a reclaimed peninsula can
+    # belong to a disconnected water component. Diagonal corner cutting is
+    # forbidden, so four-neighbour connectivity matches the search contract.
+    components, _ = label(safe)
+    reachable = components == components[start_node]
     target_inside = (xs[0]+margin < destination[0] < xs[-1]-margin
                      and ys[0]+margin < destination[1] < ys[-1]-margin
-                     and safe[goal_node])
+                     and safe[goal_node] and reachable[goal_node])
     edge = ((xx <= xs[0]+margin+2*dx) | (xx >= xs[-1]-margin-2*dx)
             | (yy <= ys[0]+margin+2*dy) | (yy >= ys[-1]-margin-2*dy))
-    exits = safe & edge & (travel >= 1000) & (progress >= min(1000., remaining*.2))
+    exits = safe & reachable & edge & (travel >= 1000) & (progress >= min(1000., remaining*.2))
     scores = np.where(exits, to_target, math.inf)
     goals = []
     if target_inside:
@@ -177,7 +190,11 @@ def coastal_corridor(terrain, spec, start, target, timeout=25.):
         scores[np.hypot(xx-xs[node[1]], yy-ys[node[0]]) < 600] = math.inf
     penalty = (1 + np.minimum(1,np.abs(coast_distance-STANDOFF_M)/STANDOFF_M)
                + 8*np.minimum(6,np.maximum(0,coast_distance-(STANDOFF_M+150))/500))
-    validator = DetourPlanner(CachedTerrain(field), [*origin, 0.], True,
+    checked = CachedTerrain(field)
+    if guard is not None:
+        from .land import LandExcludedField
+        checked = LandExcludedField(checked, guard)
+    validator = DetourPlanner(checked, [*origin, 0.], True,
                              timeout=timeout, maximum_radius=18000, maximum_length=40000)
     validator.radius = REGIONAL_BUFFER_M
     validator.clearance = MINIMUM_DEPTH_M
@@ -187,7 +204,6 @@ def coastal_corridor(terrain, spec, start, target, timeout=25.):
     def xy(node):
         return [float(xs[node[1]]), float(ys[node[0]]), 0.]
 
-    start_node = (int(np.argmin(abs(ys-origin[1]))), int(np.argmin(abs(xs-origin[0]))))
     if not safe[start_node] or not validator.segment([*origin, 0.], xy(start_node)):
         return None
     expanded = 0

src/posim_routing/land.py

@@ -0,0 +1,100 @@
+"""Mapped land is an exclusion, never evidence of water depth."""
+import json
+import math
+import os
+import zipfile
+from functools import lru_cache
+from pathlib import Path
+
+import numpy as np
+import shapely
+from shapely.ops import transform
+
+
+class LandGuard:
+    def __init__(self, land, coverage):
+        if not land.is_valid or not coverage.is_valid:
+            raise RuntimeError('invalid_land_exclusion')
+        self.land, self.coverage = land, coverage
+        shapely.prepare(self.land)
+        shapely.prepare(self.coverage)
+
+    def free(self, x, y, radius):
+        footprint = shapely.box(x-radius, y-radius, x+radius, y+radius)
+        return self.coverage.covers(footprint) and not self.land.intersects(footprint)
+
+    def segment(self, a, b, radius):
+        # Enclose the swept square used by numeric triangle clearance.
+        footprint = shapely.MultiPoint([(point[0]+dx, point[1]+dy)
+                    for point in (a,b) for dx,dy in
+                    ((-radius,-radius),(-radius,radius),(radius,radius),(radius,-radius))]).convex_hull
+        return self.coverage.covers(footprint) and not self.land.intersects(footprint)
+
+    def grid(self, xs, ys, radius=0.):
+        xx, yy = np.meshgrid(xs, ys)
+        margin = radius*math.sqrt(2)
+        return (shapely.contains_xy(self.coverage.buffer(-margin), xx, yy)
+                & ~shapely.intersects_xy(self.land.buffer(margin), xx, yy))
+
+
+class LandExcludedField:
+    def __init__(self, field, guard):
+        self.field, self.guard = field, guard
+        base = getattr(field, 'field', field)
+        self.xs, self.ys, self.heights = base.xs, base.ys, base.heights
+
+    def maximum(self, x, y, radius):
+        if not self.guard.free(x, y, radius):
+            return None
+        return self.field.maximum(x, y, radius)
+
+
+def guard_from_spec(spec, extent=None):
+    from posim_terrain.terrain_sources import cache_root
+    archive = Path(os.getenv('WWOS_COAST_OSM_ARCHIVE') or cache_root()/'coastal-land-polygons.zip')
+    if not archive.is_file():
+        raise RuntimeError('land_exclusion_data_unavailable')
+    signature = (archive.stat().st_size, archive.stat().st_mtime_ns)
+    return _load(json.dumps(spec, sort_keys=True), str(archive), signature, extent)
+
+
+def geographic_windows(inverse,half):
+    # Sample the perimeter in the actual floating frame. Unwrap around its
+    # centre before splitting queries at the international date line.
+    edges=np.linspace(-half,half,33)
+    x=np.concatenate((edges,edges,np.full(33,-half),np.full(33,half)))
+    y=np.concatenate((np.full(33,-half),np.full(33,half),edges,edges))
+    lons,lats=inverse.transform(x,y);centre=inverse.transform(0.,0.)[0]
+    lons=(np.asarray(lons)-centre+180)%360-180+centre
+    lo,hi=float(lons.min())-.01,float(lons.max())+.01
+    south,north=float(np.min(lats))-.01,float(np.max(lats))+.01
+    if hi-lo>=180 or south < -85 or north>85:
+        raise RuntimeError('land_exclusion_horizon_unavailable')
+    return [(max(-180.,lo+shift),south,min(180.,hi+shift),north)
+            for shift in (-360.,0.,360.) if lo+shift<180 and hi+shift>-180]
+
+
+@lru_cache(maxsize=12)
+def _load(encoded, archive, signature, extent):
+    import pyogrio.raw
+    from posim_terrain.terrain import projection
+    spec = json.loads(encoded)
+    _, forward, inverse = projection(spec)
+    half = max(spec['width_m'], spec['height_m'])/2 + 500 if extent is None else extent
+    windows=geographic_windows(inverse,half)
+    with zipfile.ZipFile(archive) as source:
+        member = next(name for name in source.namelist() if name.endswith('.shp'))
+    parts=[]
+    for bounds in windows:
+        metadata,_,geometries,_=pyogrio.raw.read('/vsizip/'+archive+'/'+member,
+                                               bbox=bounds,columns=[])
+        if metadata['crs']!='EPSG:4326':
+            raise RuntimeError('land_exclusion_projection_mismatch')
+        try:
+            land=shapely.union_all(shapely.intersection(shapely.from_wkb(geometries),shapely.box(*bounds)))
+            parts.append(transform(forward.transform,land))
+        except shapely.GEOSException as error:
+            raise RuntimeError('invalid_land_exclusion_source') from error
+    # Both geographic windows cover this local square; unioning their
+    # projected geographic boundaries can leave a sub-nanometre wrap seam.
+    return LandGuard(shapely.union_all(parts),shapely.box(-half,-half,half,half))

src/posim_routing/land_guidance.py

@@ -0,0 +1,75 @@
+"""Bounded mapped-water guidance; numeric depth must approve every local leg."""
+import heapq
+import math
+import time
+
+import numpy as np
+
+
+def detour(guard, start, target, step=150., radius=40., timeout=12.):
+    """Return a complete mapped-water route, including necessary initial retreat.
+
+    This exclusion search neither creates water nor approves bathymetry. Its
+    endpoints and shortcuts enclose the same swept square as the local planner.
+    """
+    deadline=time.monotonic()+timeout
+    if not guard.free(*start,radius) or not guard.free(*target,radius):
+        raise RuntimeError('mapped_guidance_endpoint_blocked')
+    if guard.segment(start,target,radius):return [list(target)]
+    pad=24000.
+    lo=np.minimum(start,target)-pad;hi=np.maximum(start,target)+pad
+    if np.max(hi-lo)>128000.:
+        raise RuntimeError('mapped_guidance_horizon_exceeded')
+    xs=np.arange(lo[0],hi[0]+step/2,step);ys=np.arange(lo[1],hi[1]+step/2,step)
+    safe=guard.grid(xs,ys,radius+step/math.sqrt(2))
+    def xy(node):return [float(xs[node[1]]),float(ys[node[0]])]
+    def anchor(point):
+        j=int(np.argmin(abs(ys-point[1])));i=int(np.argmin(abs(xs-point[0])))
+        candidates=[(jj,ii) for jj in range(max(0,j-2),min(len(ys),j+3))
+                    for ii in range(max(0,i-2),min(len(xs),i+3))]
+        candidates.sort(key=lambda n:math.dist(xy(n),point))
+        for node in candidates:
+            if safe[node] and guard.segment(point,xy(node),radius):return node
+        raise RuntimeError('mapped_guidance_anchor_blocked')
+    origin,goal=anchor(start),anchor(target)
+    costs={origin:0.};parents={};heap=[(math.dist(xy(origin),xy(goal)),0.,origin)]
+    closed=set();found=False
+    offsets=((0,1),(0,-1),(1,0),(-1,0),(1,1),(1,-1),(-1,1),(-1,-1))
+    while heap:
+        if time.monotonic()>deadline:raise RuntimeError('mapped_guidance_budget_exceeded')
+        _,cost,node=heapq.heappop(heap)
+        if node in closed:continue
+        closed.add(node)
+        if node==goal:found=True;break
+        for j,i in offsets:
+            nxt=(node[0]+j,node[1]+i)
+            if not (0<=nxt[0]<len(ys) and 0<=nxt[1]<len(xs)) or not safe[nxt] or nxt in closed:continue
+            if i and j and (not safe[node[0]+j,node[1]] or not safe[node[0],node[1]+i]):continue
+            value=cost+step*math.hypot(i,j)
+            if value>=costs.get(nxt,math.inf):continue
+            costs[nxt],parents[nxt]=value,node
+            heapq.heappush(heap,(value+math.dist(xy(nxt),xy(goal)),value,nxt))
+    if not found:raise RuntimeError('mapped_guidance_no_water_detour')
+    nodes=[goal]
+    while nodes[-1]!=origin:nodes.append(parents[nodes[-1]])
+    path=[list(start),*[xy(n) for n in reversed(nodes)],list(target)]
+    # Greedy checked shortcuts preserve the complete around-land connection.
+    result=[];index=0
+    while index<len(path)-1:
+        if time.monotonic()>deadline:raise RuntimeError('mapped_guidance_budget_exceeded')
+        chosen=index+1
+        for end in range(len(path)-1,index,-1):
+            if guard.segment(path[index],path[end],radius):chosen=end;break
+        result.append(path[chosen]);index=chosen
+    if not all(guard.segment(a,b,radius) for a,b in zip([start,*result[:-1]],result)):
+        raise RuntimeError('mapped_guidance_validation_failed')
+    if sum(math.dist(a,b) for a,b in zip([start,*result[:-1]],result))>160000.:
+        raise RuntimeError('mapped_guidance_length_exceeded')
+    return result
+
+
+def local_target(start, points, horizon=6500.):
+    """Take a checked prefix, preserving each bend before the numeric horizon."""
+    point=points[0];length=math.dist(start,point)
+    if length<=horizon:return list(point)
+    return [start[i]+(point[i]-start[i])*horizon/length for i in (0,1)]

src/posim_routing/port_approach.py

@@ -48,14 +48,17 @@ def shortest_arrival(planner, field, origin, candidates):
     return best, best_floor
 
 
-def corridor(tiles, spec, start, target):
+def corridor(tiles, spec, start, target, guard=None):
     from posim_terrain.terrain import projection
     from .navigation import TerrainField, DetourPlanner
     from .planning import CachedTerrain
     _, forward, inverse = projection(spec)
     origin = [*forward.transform(*start), 0.]
     centre = [*forward.transform(*target), 0.]
     field = CachedTerrain(TerrainField.tiles(tiles))
+    if guard is not None:
+        from .land import LandExcludedField
+        field = LandExcludedField(field, guard)
     planner = DetourPlanner(field, origin, True, step=8., timeout=20.,
                            maximum_radius=6000., maximum_length=16000.)
     # Preserve the catalog's 3 m / 25 m port endpoint screening. En route,