diff --git a/README.ko.md b/README.ko.md index 34aeacb..90e13d9 100644 --- a/README.ko.md +++ b/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` | 도착 영역 접근과 곡선 감속 | diff --git a/README.md b/README.md index 91ca44e..4796c12 100644 --- a/README.md +++ b/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 | diff --git a/pyproject.toml b/pyproject.toml index 5a0d1e2..464d61e 100644 --- a/pyproject.toml +++ b/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"} diff --git a/src/posim_routing/__init__.py b/src/posim_routing/__init__.py index 6ebeeab..930ef3d 100644 --- a/src/posim_routing/__init__.py +++ b/src/posim_routing/__init__.py @@ -1,2 +1,2 @@ """POSIM-ROUTING: reusable ocean simulation components.""" -__version__ = "0.1.1" +__version__ = "0.1.3" diff --git a/src/posim_routing/coastal.py b/src/posim_routing/coastal.py index 644762b..36a0402 100644 --- a/src/posim_routing/coastal.py +++ b/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 diff --git a/src/posim_routing/land.py b/src/posim_routing/land.py new file mode 100644 index 0000000..ffcf392 --- /dev/null +++ b/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)) diff --git a/src/posim_routing/land_guidance.py b/src/posim_routing/land_guidance.py new file mode 100644 index 0000000..d23273d --- /dev/null +++ b/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]=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 indexdeadline: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)] diff --git a/src/posim_routing/port_approach.py b/src/posim_routing/port_approach.py index a467981..57dc8f3 100644 --- a/src/posim_routing/port_approach.py +++ b/src/posim_routing/port_approach.py @@ -48,7 +48,7 @@ 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 @@ -56,6 +56,9 @@ def corridor(tiles, spec, start, target): 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, diff --git a/src/posim_routing/regional.py b/src/posim_routing/regional.py index ab79c42..7d5db43 100644 --- a/src/posim_routing/regional.py +++ b/src/posim_routing/regional.py @@ -6,13 +6,13 @@ import asyncio,math,time from pyproj import Geod GEO=Geod(ellps='WGS84') -REGIONAL_POLICY='buffered-water-v7-coast-250m-priority-2m' +REGIONAL_POLICY='buffered-water-v8-mapped-land-exclusion' REGIONAL_BUFFER_M=40. REGIONAL_DEPTH_M=2. def distance(a,b):return GEO.inv(*a,*b)[2] -def corridor(terrain,spec,start,target): +def corridor(terrain,spec,start,target,guard=None): from posim_terrain.terrain import projection from .navigation import TerrainField,DetourPlanner from .planning import CachedTerrain @@ -20,6 +20,9 @@ def corridor(terrain,spec,start,target): from posim_terrain.canal_override import in_region canal=in_region(start) field=CachedTerrain(TerrainField(terrain['xs'],terrain['ys'],terrain['heights'])) + if guard is not None: + from .land import LandExcludedField + field=LandExcludedField(field,guard) origin=[*forward.transform(*start),0.] heading,_,remaining=GEO.inv(*start,*target) deadline=time.monotonic()+60;errors=[] @@ -151,6 +154,11 @@ def update(self,state,target,needed=False,terminal=False,retain_coast=False): self.discard(state,'suez_fine_grid') point=self.active_point(state,key) plan=state.get('regional_navigation') or {} + # A mapped around-land guide must renew from its checked endpoint. It + # cannot fall back to the original straight line across that land. + guide_exhausted=(plan.get('route_key')==key and plan.get('mapped_guidance') + and plan.get('points') and plan['index']>=len(plan['points'])) + needed=needed or bool(guide_exhausted) prefer_port=terminal and plan.get('mode')!='port_approach' if terminal and point is not None and not prefer_port and not (self.task and self.task.done()):return point # Prepare the next coastal lookahead while the current checked plan @@ -205,10 +213,10 @@ def update(self,state,target,needed=False,terminal=False,retain_coast=False): if point is not None:return point plan=state.get('regional_navigation') or {} if (plan.get('route_key')==key and plan.get('policy')==REGIONAL_POLICY - and plan.get('mode') in ('coastal','coastal_approach') and plan.get('points') + and (plan.get('mode') in ('coastal','coastal_approach') or plan.get('mapped_guidance')) and plan.get('points') and plan['index']>=len(plan['points']) and distance(state['position'],plan['points'][-1])<1500 - and (nearest is not None or self.task or self.index_task)): + and (nearest is not None or self.task or self.index_task or plan.get('mapped_guidance'))): # An exhausted checked path is not permission to steer offshore. # The caller holds thrust, keeping physics and sensors running. # Retry backoff and optional empty results retain this boundary. @@ -216,6 +224,9 @@ def update(self,state,target,needed=False,terminal=False,retain_coast=False): state['regional_waiting']=dict(reason='coastal_renewal',job=plan['job'], endpoint=plan['points'][-1]) return plan['points'][-1] + if (nearest is not None or guide_exhausted or needed) and (self.task or time.monotonic()50000.: + destination=[origin[i]+(destination[i]-origin[i])*50000./length for i in (0,1)] + points=await asyncio.to_thread(detour,guard,origin,destination) + if len(points)>1: + target=list(inverse.transform(*local_target(origin,points))) + guidance=dict(points=[list(inverse.transform(*p)) for p in points], + course_target=list(course_target), + scope='mapped land exclusion only; numeric depth checked locally') value=None;coast_error=None if not canal and (nearest or needed): from .coastal import coastal_corridor - try:value=await asyncio.to_thread(coastal_corridor,grid,job['spec'],current,target) + try:value=await asyncio.to_thread(coastal_corridor,grid,job['spec'],current,target,guard=guard) except RuntimeError as error:coast_error=str(error)[:160] if value is None: current=list(state['position']) @@ -282,7 +317,7 @@ async def prepare(self,state,target,key,nearest=None,needed=True): try: # This small step is checked against real numeric depth, # then checked again against the live Gazebo collision mesh. - value=await asyncio.to_thread(corridor,grid,job['spec'],current,approach) + value=await asyncio.to_thread(corridor,grid,job['spec'],current,approach,guard=guard) endpoint=value['points'][-1] if (distance(current,target)-distance(endpoint,target)<500 or distance(current,nearest[1])-distance(endpoint,nearest[1])<500): @@ -293,12 +328,14 @@ async def prepare(self,state,target,key,nearest=None,needed=True): except RuntimeError as error: coast_error=str(error)[:160] if value is None: - value=await asyncio.to_thread(corridor,grid,job['spec'],current,target) if needed else dict(points=[],index=0,mode='course_retained') + value=await asyncio.to_thread(corridor,grid,job['spec'],current,target,guard=guard) if needed or nearest else dict(points=[],index=0,mode='course_retained') if value['points']:value['mode']='detour' if canal and value['points']:value.update(mode='canal',assumed_depth_m=10,measured_bathymetry=False) if coast_error:value['coastal_preference_error']=coast_error + if guidance:value['mapped_guidance']=guidance if not value['points']:value['coastal_preference_error']=coast_error or 'no_verified_progressive_coastal_corridor' - value.update(route_key=key,job=job['id'],source=job['manifest']['source'],policy=REGIONAL_POLICY) + value.update(route_key=key,job=job['id'],source=job['manifest']['source'],policy=REGIONAL_POLICY, + mapped_land_checked=guard is not None) value['probe_position']=current return value async def close(self): diff --git a/tests/test_coastal_navigation.py b/tests/test_coastal_navigation.py index b653f5b..7584cdd 100644 --- a/tests/test_coastal_navigation.py +++ b/tests/test_coastal_navigation.py @@ -155,6 +155,28 @@ def test_real_progressive_shore_excludes_backward_coast(self): class BackgroundPlanning(unittest.IsolatedAsyncioTestCase): + async def asyncSetUp(self): + from unittest.mock import patch + import shapely + from posim_routing.land import LandGuard + # These RPC fixtures carry synthetic numeric terrain, not map data. + fixture=LandGuard(shapely.GeometryCollection(),shapely.box(-100000,-100000,100000,100000)) + mock=patch('posim_routing.land.guard_from_spec',return_value=fixture) + mock.start();self.addCleanup(mock.stop) + + async def test_exhausted_mapped_detour_cannot_fall_back_across_the_peninsula(self): + async def rpc(*args,**kwargs):await asyncio.Event().wait() + key=[0,0,1];end=geo(2000,0);target=geo(-8000,-8000) + sample=dict(lap=0,leg=0,waypoint=1,position=end,status='sailing',native={}, + regional_navigation=dict(route_key=key,points=[end],index=1,mode='detour', + mapped_guidance=dict(course_target=target),policy=REGIONAL_POLICY, + probe_position=geo(0,0),job='checked')) + navigator=RegionalNavigator(rpc);navigator.key=key + self.assertEqual(navigator.update(sample,target),end) + self.assertIsNotNone(sample.get('regional_waiting')) + self.assertIsNotNone(navigator.task) + await navigator.close() + async def test_coastal_lookahead_is_prefetched_before_its_end(self): async def rpc(*args,**kwargs): await asyncio.Event().wait() @@ -314,7 +336,8 @@ async def rpc(path, body=None, timeout=10): target = geo(6500, 9000) point = geo(800, 2000) with patch('posim_routing.coastal.coastal_corridor', return_value=dict(points=[point], index=0, mode='coastal')): - self.assertEqual(navigator.update(state, target), target) + self.assertEqual(navigator.update(state, target), state['position']) + self.assertEqual(state['regional_waiting']['reason'],'land_checked_coastal_route_pending') self.assertEqual(state['status'], 'sailing') await navigator.task self.assertEqual(navigator.update(state, target), point) @@ -351,7 +374,7 @@ async def rpc(path, body=None, timeout=10): return terrain return dict(id='terrain-fixture', status='ready', spec=SPEC, manifest=dict(source='synthetic fixture')) - def plan(terrain, spec, start, target): + def plan(terrain, spec, start, target, guard=None): starts.append(start) return dict(points=[geo(800, 2000)], index=0, mode='coastal') navigator = RegionalNavigator(rpc) diff --git a/tests/test_land_exclusion.py b/tests/test_land_exclusion.py new file mode 100644 index 0000000..95c8f43 --- /dev/null +++ b/tests/test_land_exclusion.py @@ -0,0 +1,115 @@ +"""Negative bathymetry cannot override mapped land. Synthetic regressions.""" +import unittest +import math +import tempfile +import zipfile +from pathlib import Path +from unittest.mock import patch + +import numpy as np +import shapely + +from posim_routing.land import LandGuard, LandExcludedField, guard_from_spec +from posim_routing.navigation import TerrainField, DetourPlanner +from posim_routing.planning import CachedTerrain +from posim_routing.coastal import coastal_corridor +from posim_routing.land_guidance import detour,local_target +from posim_terrain.terrain import projection + +SPEC=dict(latitude=35.,longitude=129.,width_m=20000.,height_m=20000.,samples=257) +_,FORWARD,INVERSE=projection(SPEC) + + +class LandExclusion(unittest.TestCase): + def test_deep_numeric_cache_cannot_admit_reclaimed_land(self): + axis=np.arange(-1000.,1001.,50.) + numeric=TerrainField(axis,axis,np.full((len(axis),len(axis)),-20.)) + guard=LandGuard(shapely.box(-200,-300,200,300),shapely.box(-1100,-1100,1100,1100)) + field=LandExcludedField(CachedTerrain(numeric),guard) + self.assertIsNone(field.maximum(0,0,8.)) + self.assertEqual(field.maximum(-800,0,8.),-20.) + planner=DetourPlanner(field,[-800,0,0],True,step=50,timeout=5,maximum_radius=2000) + planner.radius=8. + route=planner.route([-800,0,0],[800,0,0]) + self.assertGreater(max(abs(p[1]) for p in route),300.) + self.assertTrue(all(guard.segment(a,b,8.) for a,b in zip([[-800,0,0],*route[:-1]],route))) + np.testing.assert_array_equal(numeric.heights,-20.) + + def test_thin_barrier_between_grid_vertices_is_blocked(self): + axis=np.arange(-200.,201.,100.) + guard=LandGuard(shapely.box(44,-180,46,180),shapely.box(-1000,-1000,1000,1000)) + self.assertFalse(guard.segment([-100,0],[100,0],8.)) + self.assertFalse(guard.grid(axis,axis,80.)[2,2]) + self.assertTrue(guard.segment([-100,-300],[100,-300],8.)) + + def test_land_extends_coastal_search_when_numeric_field_calls_it_water(self): + axis=np.linspace(-10000.,10000.,257);xx,yy=np.meshgrid(axis,axis) + numeric=np.where(xx<=0,20.,-20.) + grid=dict(xs=axis.tolist(),ys=axis.tolist(),heights=numeric.tolist()) + land=shapely.union_all([shapely.box(-11000,-11000,0,11000),shapely.box(0,-1000,1800,1500)]) + guard=LandGuard(land,shapely.box(-11000,-11000,11000,11000)) + start=list(INVERSE.transform(800,-4000));target=list(INVERSE.transform(6500,9000)) + plan=coastal_corridor(grid,SPEC,start,target,guard=guard) + self.assertIsNotNone(plan) + points=[[*FORWARD.transform(*p),0] for p in [start,*plan['points'],*plan['onward_points']]] + self.assertGreater(max(p[0] for p in points),1800.) + self.assertTrue(all(guard.segment(a,b,40.) for a,b in zip(points,points[1:]))) + + def test_unknown_coverage_and_missing_dataset_are_rejected(self): + guard=LandGuard(shapely.GeometryCollection(),shapely.box(-100,-100,100,100)) + self.assertFalse(guard.segment([0,0],[200,0],8.)) + with patch.dict('os.environ',{'WWOS_COAST_OSM_ARCHIVE':'/nonexistent/land.zip'}): + with self.assertRaisesRegex(RuntimeError,'data_unavailable'): + guard_from_spec(SPEC) + + def test_nearer_disconnected_exits_do_not_hide_reachable_water_exit(self): + axis=np.linspace(-10000.,10000.,257) + grid=dict(xs=axis.tolist(),ys=axis.tolist(),heights=np.full((257,257),-20.).tolist()) + guard=LandGuard(shapely.box(-200,-11000,0,11000),shapely.box(-11000,-11000,11000,11000)) + start=list(INVERSE.transform(800,0));target=list(INVERSE.transform(-20000,-20000)) + plan=coastal_corridor(grid,SPEC,start,target,guard=guard) + self.assertIsNotNone(plan) + exit_x,exit_y=FORWARD.transform(*plan['onward_exit']) + self.assertGreater(exit_x,0.) + self.assertLess(exit_y,-9000.) + + def test_complete_water_detour_can_require_initial_distance_loss(self): + # A bay opens east, while the onward course is west of its banks. + land=shapely.union_all([shapely.box(-2500,-6000,-1500,6000), + shapely.box(-2500,5000,4000,6000), + shapely.box(-2500,-6000,4000,-5000)]) + guard=LandGuard(land,shapely.box(-30000,-30000,30000,30000)) + start=[0,0];goal=[-8000,-8000] + points=detour(guard,start,goal) + self.assertGreater(points[0][0],4000.) + self.assertGreater(math.dist(points[0],goal),math.dist(start,goal)) + self.assertEqual(points[-1],goal) + self.assertTrue(all(guard.segment(a,b,40.) for a,b in zip([start,*points[:-1]],points))) + local=local_target(start,points,1000.) + self.assertAlmostEqual(math.dist(start,local),1000.) + self.assertTrue(guard.segment(start,local,40.)) + + def test_mapped_guidance_cannot_place_a_goal_on_land(self): + guard=LandGuard(shapely.box(-1000,-1000,1000,1000),shapely.box(-30000,-30000,30000,30000)) + with self.assertRaisesRegex(RuntimeError,'endpoint_blocked'): + detour(guard,[-5000,0],[0,0]) + + def test_land_dataset_queries_both_sides_of_date_line(self): + spec={**SPEC,'longitude':179.999} + _,forward,inverse=projection(spec) + west=shapely.box(179.997,34.997,180.,35.003) + east=shapely.box(-180.,34.997,-179.997,35.003) + def read(*args,bbox,**kwargs): + geometry=west if bbox[0]>0 else east + return {'crs':'EPSG:4326'},None,np.array([geometry.wkb],dtype=object),None + with tempfile.TemporaryDirectory() as folder: + archive=Path(folder)/'land.zip' + with zipfile.ZipFile(archive,'w') as out:out.writestr('land.shp',b'fixture') + with patch.dict('os.environ',{'WWOS_COAST_OSM_ARCHIVE':str(archive)}),patch('pyogrio.raw.read',side_effect=read) as loader: + guard=guard_from_spec(spec) + self.assertEqual(loader.call_count,2) + self.assertFalse(guard.segment([-1000,0],[1000,0],8.)) + self.assertTrue(guard.segment([-1000,1000],[1000,1000],8.)) + + +if __name__=='__main__':unittest.main()