""" S-57 ENC chart manager — parse, store and serve chart features as GeoJSON. """ import json import logging import math import shutil import tempfile import zipfile from pathlib import Path log = logging.getLogger(__name__) # Optional heavy deps — only needed for parsing .000 ENC files. # GPS Navigator may not have them installed; CSV-based charts still work. try: import geopandas as gpd import pandas as pd _HAS_GPD = True except ImportError: gpd = None pd = None _HAS_GPD = False log.warning("geopandas/pandas not installed — .000 ENC parsing disabled. " "Install with: pip install geopandas pyogrio") CHARTS_DIR = Path(__file__).parent.parent / "charts" CHARTS_DIR.mkdir(exist_ok=True) # S-57 layers we care about, mapped to a friendly category S57_LAYERS = { "LIGHTS": "light", "BOYLAT": "buoy", "BOYCAR": "buoy", "BOYISD": "buoy", "BOYSPP": "buoy", "BOYSAW": "buoy", "BCNLAT": "beacon", "BCNCAR": "beacon", "BCNISD": "beacon", "BCNWTW": "beacon", "LNDMRK": "landmark", "RTPBCN": "beacon", "LAKARE": "area", # Lake / lagoon / ciénaga area # BUISGL is rendered as polygon footprint via _parse_land → landLayer. # Removed from S57_LAYERS to avoid a duplicate centroid-point beacon dot. } # Human-readable light character codes LITCHR = { 1:"F", 2:"Fl", 3:"LFl", 4:"Q", 5:"VQ", 6:"UQ", 7:"Iso", 8:"Oc", 9:"IQ", 10:"IVQ", 11:"IUQ", 12:"Mo", 13:"FFl", 14:"FlLFl", 15:"OcFl", 16:"FLFl", 25:"Al.Oc", 26:"Al.LFl", 27:"Al.Fl", 28:"Al.Grp", } COLOUR = { 1:"W", 2:"K", 3:"R", 4:"G", 5:"B", 6:"Y", 7:"Gy", 8:"Br", 9:"Amb", 10:"Vi", 11:"Or", 12:"Mg", } BOYSHP = { 1:"Conical(N)",2:"Can(C)",3:"Spherical",4:"Pillar", 5:"Spar",6:"Barrel",7:"Super-buoy",8:"Ice buoy", } # Cardinal category (CATCAM) per IHO S-57 Object Catalogue CATCAM = {1:"N", 2:"E", 3:"S", 4:"W"} # IALA region defaults by chart producer prefix (cell ID first 2 chars) REGION_A_PREFIXES = { "GB","FR","DE","NL","BE","DK","NO","SE","FI","IS","IE","ES","PT","IT", "GR","TR","RU","UA","HR","CY","MT","EE","LV","LT","PL","RO","BG","SI", "MA","DZ","TN","LY","EG","ZA","NG","SN","CI","KE","TZ","MZ", "IN","PK","BD","LK","MM","TH","ID","MY","SG","VN","CN","HK","TW", "AE","SA","OM","IR","IQ","KW","QA","BH","YE", "AU","NZ","PG", } REGION_B_PREFIXES = { "US","CA","MX","BS","CU","JM","HT","DO","PR","VI", "BZ","GT","SV","HN","NI","CR","PA", "CO","VE","EC","PE","CL","AR","UY","BR","GY","SR","PY","BO", "JP","KR","KP","PH", } def auto_region(cell_id: str) -> str: prefix = (cell_id or "").upper()[:2] if prefix in REGION_A_PREFIXES: return "A" if prefix in REGION_B_PREFIXES: return "B" return "B" # ── CSV → features (DIMAR / custom CSV charts) ──────────────────────────────── # Mapping from DIMAR CSV feat_type to canonical S-57 layer name. _CSV_FEAT_TYPE_MAP = { "BOYSPEC": "BOYSPP", # DIMAR uses BOYSPEC; IHO S-57 is BOYSPP } # Feature categories per layer _CSV_CATEGORY_MAP = { "BOYLAT": "buoy", "BOYCAR": "buoy", "BOYISD": "buoy", "BOYSPP": "buoy", "BOYSAW": "buoy", "BCNLAT": "beacon","BCNCAR": "beacon", "LIGHTS": "light", "LNDMRK": "landmark", } def _csv_num(v): if v is None or str(v).strip() == "": return None try: return float(v) except Exception: return None def _csv_int(v): n = _csv_num(v) return int(n) if n is not None else None def _csv_fmt(v): """Format float: drop trailing .0 for whole numbers.""" try: f = float(v) return str(int(f)) if f == int(f) else str(f) except Exception: return str(v) def _csv_light_desc(litchr, siggrp, colour, sigper, valnmr) -> str: """Build compact light description from CSV fields (e.g. 'Fl G 3s 3M').""" parts = [] lc_int = _csv_int(litchr) lc = LITCHR.get(lc_int, str(lc_int)) if lc_int is not None else "" sg = _csv_int(siggrp) if sg is not None: lc = f"{lc}({sg})" col_int = _csv_int(colour) col_str = COLOUR.get(col_int, "") if col_int is not None else "" if lc: parts.append(f"{lc} {col_str}".strip()) sp = _csv_num(sigper) if sp is not None: parts.append(f"{_csv_fmt(sp)}s") rng = _csv_num(valnmr) if rng is not None: parts.append(f"{_csv_fmt(rng)}M") return " ".join(parts) def _csv_infer_catcam(siggrp: str, name: str) -> str | None: """ Cardinal buoy quadrant from SIGGRP (reliable) or DIMAR naming convention. Q(9)+LFl=W Q(6)+LFl=S Q(3)=E Q=N DIMAR name suffixes: SS/VS→S SN/VN→N SE→E SO→W """ sg = _csv_int(siggrp) if sg is not None: if sg == 9: return "W" if sg == 6: return "S" if sg == 3: return "E" return "N" n = (name or "").upper() if any(x in n for x in (" SS", "(SS)", "VS")): return "S" if any(x in n for x in (" SE", "(SE)")): return "E" if any(x in n for x in (" SO", "(SO)")): return "W" return "N" def _csv_infer_catlam(colour: str) -> int | None: """IALA B (Americas): green(4)=port=1, red(3)=stbd=2.""" c = _csv_int(colour) if c == 4: return 1 if c == 3: return 2 return None def _parse_csvs_to_features(csv_dir: Path) -> list[dict]: """ Read navigation-aid CSV files from csv_dir and return GeoJSON features. Each CSV must have columns: OBJNAM, lon, lat, feat_type, LITCHR, SIGGRP, COLOUR, SIGPER, VALNMR, HEIGHT, ORIENT, INFORM. feat_type values: BOYLAT, BOYCAR, BOYISD, BOYSPEC/BOYSPP, BCNLAT, LIGHTS, LNDMRK. This function is the primary path for DIMAR/custom charts created with QGISS57Converter; it preserves all light attributes (LITCHR, SIGPER, etc.) that the GDAL S-57 driver may drop during the .000 round-trip. """ import csv as _csv_mod features: list[dict] = [] for csv_file in sorted(csv_dir.glob("*.csv")): with open(csv_file, newline="", encoding="utf-8") as fh: reader = _csv_mod.DictReader(fh) for row in reader: feat_type = (row.get("feat_type") or csv_file.stem).strip() layer = _CSV_FEAT_TYPE_MAP.get(feat_type, feat_type) if not layer: continue lon = _csv_num(row.get("lon")) lat = _csv_num(row.get("lat")) if lon is None or lat is None: continue category = _CSV_CATEGORY_MAP.get(layer, "buoy") litchr = row.get("LITCHR", "").strip() siggrp = row.get("SIGGRP", "").strip() colour = row.get("COLOUR", "").strip() sigper = row.get("SIGPER", "").strip() valnmr = row.get("VALNMR", "").strip() height = row.get("HEIGHT", "").strip() orient_r = row.get("ORIENT", "").strip() name = row.get("OBJNAM", "").strip() inform = row.get("INFORM", "").strip() col_int = _csv_int(colour) colours = [col_int] if col_int is not None else [] light_d = _csv_light_desc(litchr, siggrp, colour, sigper, valnmr) props: dict = { "layer": layer, "category": category, "name": name or None, "info": inform or None, "light_desc": light_d or None, "range_nm": _csv_num(valnmr), "height_m": _csv_num(height), "colours": colours, "colour_code": col_int, } if layer in ("BOYCAR", "BCNCAR"): props["catcam"] = _csv_infer_catcam(siggrp, name) if layer in ("BOYLAT", "BCNLAT"): props["catlam"] = _csv_infer_catlam(colour) orient_val = _csv_num(orient_r) if orient_val is not None: props["orient"] = orient_val props["aid_type"] = classify(layer, props) # Remove None values props = {k: v for k, v in props.items() if v is not None} features.append({ "type": "Feature", "geometry": {"type": "Point", "coordinates": [lon, lat]}, "properties": props, }) return features def install_from_csv_dir(csv_dir: Path, cell_id: str) -> str: """ Create or update an AidsMonitoring chart cell from a directory of CSV files. This is the preferred install pathway for custom (DIMAR) charts because it preserves all S-57 attribute codes (LITCHR, SIGGRP, etc.) without loss. csv_dir — directory containing *.csv files (BOYLAT.csv, BOYCAR.csv, etc.) cell_id — chart cell identifier (e.g. 'BAHIA_DE_CARTAGENA') """ cell_id = cell_id.upper() cell_dir = CHARTS_DIR / cell_id cell_dir.mkdir(exist_ok=True) features = _parse_csvs_to_features(csv_dir) with open(cell_dir / "features.geojson", "w", encoding="utf-8") as f: json.dump({"type": "FeatureCollection", "features": features}, f, ensure_ascii=False) log.info("CSV chart %s → %d features", cell_id, len(features)) # Empty auxiliary caches so the frontend doesn't show stale data for fname in ("depths.geojson", "land.geojson", "hazards.geojson", "zones.geojson"): p = cell_dir / fname if not p.exists(): with open(p, "w") as f: json.dump({"type": "FeatureCollection", "features": []}, f) # Update meta.json bbox = None if features: lons = [ft["geometry"]["coordinates"][0] for ft in features] lats = [ft["geometry"]["coordinates"][1] for ft in features] bbox = [min(lons), min(lats), max(lons), max(lats)] meta = get_meta(cell_id) meta["feature_count"] = len(features) meta["bbox"] = bbox meta.setdefault("region", auto_region(cell_id)) _meta_path(cell_dir).write_text(json.dumps(meta)) return cell_id def install_from_csv_zip(zip_path: Path) -> list[str]: """ Install one or more CSV chart cells from a ZIP archive. Expected ZIP layout (any of these is accepted): Option A — single cell (cell_id inferred from folder name or ZIP name): BOYLAT.csv BOYCAR.csv ... Option B — cell_id in subfolder: BAHIA_DE_CARTAGENA/BOYLAT.csv BAHIA_DE_CARTAGENA/BOYCAR.csv Option C — meta.json declares cell_id: meta.json → {"cell_id": "BAHIA_DE_CARTAGENA"} BOYLAT.csv ... """ import csv as _csv_mod installed: list[str] = [] with zipfile.ZipFile(zip_path) as z: namelist = z.namelist() # Collect all CSV files grouped by subfolder csv_files = [n for n in namelist if n.lower().endswith(".csv")] if not csv_files: raise ValueError("No CSV files found in ZIP") # Check for meta.json meta_cell_id: str | None = None if "meta.json" in namelist: try: meta_cell_id = json.loads(z.read("meta.json")).get("cell_id") except Exception: pass # Group by directory prefix import collections groups: dict[str, list[str]] = collections.defaultdict(list) for name in csv_files: prefix = str(Path(name).parent) groups[prefix].append(name) with tempfile.TemporaryDirectory() as tmp_root: tmp_root_p = Path(tmp_root) for prefix, members in groups.items(): # Determine cell_id if meta_cell_id: cid = meta_cell_id elif prefix and prefix != ".": cid = Path(prefix).name else: cid = Path(zip_path).stem # e.g. BAHIA_DE_CARTAGENA # Extract CSVs to temp dir tmp_csv = tmp_root_p / cid tmp_csv.mkdir(exist_ok=True) for member in members: data = z.read(member) (tmp_csv / Path(member).name).write_bytes(data) # Install result_id = install_from_csv_dir(tmp_csv, cid) installed.append(result_id) return installed def _meta_path(cell_dir: Path) -> Path: return cell_dir / "meta.json" def get_meta(cell_id: str) -> dict: cell_dir = CHARTS_DIR / cell_id.upper() p = _meta_path(cell_dir) if p.exists(): try: return json.loads(p.read_text()) except Exception: pass return {} def set_meta(cell_id: str, **patch) -> dict: cell_dir = CHARTS_DIR / cell_id.upper() if not cell_dir.exists(): raise FileNotFoundError(cell_id) meta = get_meta(cell_id) meta.update(patch) _meta_path(cell_dir).write_text(json.dumps(meta)) return meta def get_region(cell_id: str) -> str: return get_meta(cell_id).get("region") or auto_region(cell_id) # ── Aid classification ──────────────────────────────────────────────────── # Canonical aid types — derived from S-57 layer + attributes. # Frontend dispatches symbol drawing on this; renders are region-aware # only for LATERAL_* (which depend on IALA A vs B for color). def classify(layer: str, props: dict) -> str: catlam = props.get("catlam") catcam = props.get("catcam") if layer in ("BOYLAT", "BCNLAT"): if catlam == 1: return "LATERAL_PORT" if catlam == 2: return "LATERAL_STBD" if catlam == 3: return "LATERAL_PREF_STBD" # preferred channel to STBD (acts port) if catlam == 4: return "LATERAL_PREF_PORT" # preferred channel to PORT (acts stbd) # Fallback: infer from primary colour (IALA B — Americas: green=port, red=stbd) colour = props.get("colour_code") if colour == 4: return "LATERAL_PORT" # green if colour == 3: return "LATERAL_STBD" # red return "LATERAL_UNKNOWN" if layer in ("BOYCAR", "BCNCAR"): if catcam == "N": return "CARDINAL_N" if catcam == "E": return "CARDINAL_E" if catcam == "S": return "CARDINAL_S" if catcam == "W": return "CARDINAL_W" return "CARDINAL_UNKNOWN" if layer in ("BOYISD", "BCNISD"): return "ISOLATED_DANGER" if layer == "BOYSAW": return "SAFE_WATER" if layer == "BOYSPP": return "SPECIAL" if layer == "LIGHTS": return "LIGHT_POINT" if layer == "LNDMRK": return "LANDMARK" if layer == "RTPBCN": return "RACON" if layer.startswith("BCN"): return "BEACON_GENERIC" if layer.startswith("BOY"): return "BUOY_GENERIC" return "UNKNOWN" # Topmark shape (TOPSHP) TOPSHP = { 1:"cone_up",2:"cone_down",3:"sphere",4:"X",5:"cross", 6:"can",7:"T",8:"triangle_up",9:"triangle_down", 10:"circle",11:"rectangle",12:"board",13:"diamond", 14:"diamond_2",15:"square",16:"rectangle_2",17:"pennant", 18:"pennant_2",19:"pennant_3",20:"pennant_4", 21:"2sphere",22:"2cone_up",23:"2cone_down", 24:"cone_up_sphere",25:"2cone_base_up",26:"2cone_base_down", 27:"2cone_point",28:"besom_up",29:"besom_down", 30:"flag",31:"sphere_2",32:"2cone_base_apart", 33:"board_2", } def _safe_list(val) -> list: """Return a plain Python list from a pandas/numpy array-like value. GDAL/pyogrio surfaces multi-valued S-57 attributes (like COLOUR) as comma-separated STRINGS (e.g. "4" or "1,4"), not as Python lists. Without the string branch we silently drop colour info.""" if val is None: return [] try: if hasattr(val, 'tolist'): val = val.tolist() if isinstance(val, str): return [int(x) for x in val.split(",") if x.strip()] if isinstance(val, (list, tuple)): return [int(v) for v in val if v is not None] return [int(val)] except Exception: return [] def _safe(val): """Return a plain Python scalar from any pandas/numpy value.""" if val is None: return None # numpy array or list with multiple elements — take first if hasattr(val, '__len__') and not isinstance(val, str): if len(val) == 0: return None val = val[0] # numpy scalar if hasattr(val, 'item'): try: val = val.item() except Exception: val = val.tolist() # NaN check after extraction try: if isinstance(val, float) and _HAS_GPD and pd.isna(val): return None except Exception: pass return val def _light_desc(row) -> str: parts = [] chr_code = _safe(row.get("LITCHR")) grp = _safe(row.get("SIGGRP")) col_list = _safe(row.get("COLOUR")) period = _safe(row.get("SIGPER")) rng = _safe(row.get("VALNMR")) lc = LITCHR.get(int(chr_code), str(chr_code)) if chr_code else "" if grp: lc = f"{lc}({grp})" colours = "" if col_list: try: nums = [int(x) for x in str(col_list).split(",") if x.strip()] colours = ".".join(COLOUR.get(n, "?") for n in nums) except Exception: pass if lc: parts.append(f"{lc} {colours}".strip()) if period: parts.append(f"{period}s") if rng: parts.append(f"{rng}M") return " ".join(parts) def _parse_cell(enc_path: Path) -> list[dict]: if not _HAS_GPD: raise RuntimeError("geopandas not installed — cannot parse .000 ENC files. " "Run: pip install geopandas pyogrio") features = [] # ── Build land union for LDLINE clipping ───────────────────────────────── # Leading lines must stop at the riverbank / land boundary, not extend # indefinitely. We union all land polygons in the cell and use them to # clip each LDLINE so it terminates where it hits the shore. from shapely.ops import unary_union from shapely.geometry import LineString as SLineString, Point as SPoint _land_clip = None for _ll in ("LNDARE", "BUAARE"): try: _gdf_l = gpd.read_file(str(enc_path), layer=_ll, engine="pyogrio") if _gdf_l.empty: continue _gdf_l = _gdf_l.to_crs(epsg=4326) _polys = [g for g in _gdf_l.geometry if g is not None and not g.is_empty and g.geom_type in ("Polygon", "MultiPolygon")] if _polys: _u = unary_union(_polys).buffer(0) _land_clip = _u if _land_clip is None else _land_clip.union(_u) except Exception: pass for layer_name, category in S57_LAYERS.items(): try: gdf = gpd.read_file(str(enc_path), layer=layer_name, engine="pyogrio") except Exception: continue if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue # Use centroid for all geometry types pt = geom.centroid if geom.geom_type != "Point" else geom colours = _safe_list(row.get("COLOUR")) # e.g. [3,4] = Red+Green props = { "layer": layer_name, "category": category, "name": _safe(row.get("OBJNAM")), "info": _safe(row.get("INFORM")), "height_m": _safe(row.get("HEIGHT")), "range_nm": _safe(row.get("VALNMR")), "status": _safe(row.get("STATUS")), "colours": colours, # list of S-57 colour codes "colour_code": colours[0] if colours else None, # primary colour } # Extract light description for ALL navigational aids and lights. # Previously only BOY* layers were included, so BCNLAT/BCNWTW/LNDMRK # silently lost their LITCHR/SIGPER/VALNMR — fixed here. if category in ("light", "buoy", "beacon", "landmark"): props["light_desc"] = _light_desc(row) if layer_name.startswith("BOY") or layer_name.startswith("BCN"): shape = _safe(row.get("BOYSHP")) props["boyshp"] = int(shape) if shape else None # raw code for JS props["shape"] = BOYSHP.get(int(shape), None) if shape else None topshp = _safe(row.get("TOPSHP")) props["topshp"] = int(topshp) if topshp else None if layer_name in ("BOYCAR", "BCNCAR"): catcam = _safe(row.get("CATCAM")) props["catcam"] = CATCAM.get(int(catcam), None) if catcam else None if layer_name in ("BOYLAT", "BCNLAT"): catlam = _safe(row.get("CATLAM")) # 1=port, 2=stbd, 3=pref-stbd, 4=pref-port props["catlam"] = int(catlam) if catlam else None # ORIENT — bearing of the leading line in True degrees orient_raw = _safe(row.get("ORIENT")) if orient_raw: try: props["orient"] = float(orient_raw) except ValueError: pass props["aid_type"] = classify(layer_name, props) features.append({ "type": "Feature", "geometry": {"type": "Point", "coordinates": [pt.x, pt.y]}, "properties": {k: v for k, v in props.items() if v is not None}, }) # ── Leading-line bearing projection for lights with ORIENT ──────── # Project a line from the mark in the ORIENT direction, then clip it # against the land union so it terminates at the riverbank/shoreline. # A long initial line (5 NM) is used so the clip always finds the shore. if "orient" in props: bearing_deg = props["orient"] bearing_rad = math.radians(bearing_deg) lat_rad = math.radians(pt.y) cos_lat = math.cos(lat_rad) # Initial generous length (5 NM) — will be clipped to land boundary PROBE_NM = 5.0 d_lat = (PROBE_NM / 60.0) * math.cos(bearing_rad) d_lon = (PROBE_NM / 60.0) * math.sin(bearing_rad) / max(cos_lat, 1e-6) end = [pt.x + d_lon, pt.y + d_lat] # Densify the probe line (10 intermediate points) so the clip is accurate n_pts = 12 ldcoords = [ [pt.x + d_lon * i / n_pts, pt.y + d_lat * i / n_pts] for i in range(n_pts + 1) ] # Clip against land — the leading line must not cross land if _land_clip is not None: try: probe = SLineString(ldcoords) in_water = probe.difference(_land_clip) if in_water.is_empty: # Mark is on land — use minimum stub ldcoords = [[pt.x, pt.y], end] else: # Keep only the segment that touches the mark start_pt = SPoint(pt.x, pt.y) geoms = (list(in_water.geoms) if in_water.geom_type.startswith("Multi") else [in_water]) # Sort by distance of each segment's first point to mark geoms.sort(key=lambda g: start_pt.distance(SPoint(g.coords[0]))) # Take the first (nearest to mark) segment seg = geoms[0] ldcoords = [list(c[:2]) for c in seg.coords] except Exception: pass # fall through to unclipped coords features.append({ "type": "Feature", "geometry": { "type": "LineString", "coordinates": ldcoords, }, "properties": { "layer": "LDLINE", "category": "ldline", "aid_type": "LEADING_LINE", "name": props.get("name", ""), "orient": bearing_deg, "colours": colours, }, }) # ── Proximity merge: copy light attrs from LIGHTS → co-located buoys/beacons ── # In S-57 a LIGHTS object and a BOYLAT object may share the same coordinates. # The LIGHTS carries LITCHR/VALNMR/HEIGHT/COLOUR; the BOYLAT carries BOYSHP/CATLAM. # We find the nearest LIGHTS within 50 m and backfill missing attrs on the host. MERGE_DEG = 0.00045 # ≈ 50 m at equator (1° lat ≈ 111 km) light_feats = [f for f in features if f["properties"].get("layer") == "LIGHTS" and f["geometry"]["type"] == "Point"] struct_feats = [f for f in features if (f["properties"].get("layer","").startswith("BOY") or f["properties"].get("layer","").startswith("BCN")) and f["geometry"]["type"] == "Point"] for lf in light_feats: lc = lf["geometry"]["coordinates"] lp = lf["properties"] best_host, best_dist = None, MERGE_DEG for bf in struct_feats: bc = bf["geometry"]["coordinates"] d = ((lc[0] - bc[0])**2 + (lc[1] - bc[1])**2) ** 0.5 if d < best_dist: best_dist, best_host = d, bf if best_host: hp = best_host["properties"] if lp.get("light_desc") and not hp.get("light_desc"): hp["light_desc"] = lp["light_desc"] if lp.get("range_nm") and not hp.get("range_nm"): hp["range_nm"] = lp["range_nm"] if lp.get("height_m") and not hp.get("height_m"): hp["height_m"] = lp["height_m"] if lp.get("colours") and not hp.get("colours"): hp["colours"] = lp["colours"] hp["colour_code"] = lp.get("colour_code") # Also sync orient for leading lines if the buoy didn't have it if lp.get("orient") is not None and hp.get("orient") is None: hp["orient"] = lp["orient"] return features def _read_depth_unit(enc_path: Path) -> str: """Read DSPM_DUNI from the S-57 dataset metadata. 1=METERS, 2=FEET, 3=FATHOMS.""" try: from osgeo import gdal ds = gdal.OpenEx(str(enc_path)) if ds: md = ds.GetMetadata() or {} duni = md.get("DSPM_DUNI") or md.get("S57_DUNI") if duni: return {"1": "METERS", "2": "FEET", "3": "FATHOMS"}.get(str(duni), "METERS") except Exception as e: log.debug("Could not read DUNI from %s: %s", enc_path, e) return "METERS" # NOAA / IHO default def _bbox_of_coords(coords) -> list[float]: xs, ys = [], [] def walk(c): if isinstance(c[0], (int, float)): xs.append(c[0]); ys.append(c[1]) else: for sub in c: walk(sub) walk(coords) return [min(xs), min(ys), max(xs), max(ys)] def _parse_depths(enc_path: Path) -> list[dict]: """Extract bathymetry from S-57 cell. SOUNDG: MultiPoint Z (sounding points; depth in geometry's Z value) DEPCNT: LineString (contour at VALDCO meters) DEPARE: Polygon (depth area between DRVAL1 and DRVAL2 meters) Each feature gets a `bbox` property [minx, miny, maxx, maxy] so the /charts/depths endpoint can do fast viewport filtering without touching the geometries. The depth unit (METERS/FEET/FATHOMS) read from the cell's DSPM record is attached to every feature so the frontend legend can display it without a separate request. """ if not _HAS_GPD: return [] out = [] unit = _read_depth_unit(enc_path) # METERS / FEET / FATHOMS # Build a UNIFIED land mask from every S-57 layer that represents # something a depth feature shouldn't cover. NOAA puts urban islands / # reclaimed land in BUAARE, dykes in DAMCON, dry docks in DRYDOC, etc. # — using only LNDARE leaves visible bleed. from shapely.ops import unary_union from shapely.prepared import prep LAND_LAYERS = ("LNDARE", "BUAARE", "DAMCON", "DRYDOC", "PONTON", "LNDRGN") land_polys = [] for ll in LAND_LAYERS: try: gdf = gpd.read_file(str(enc_path), layer=ll, engine="pyogrio") if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) land_polys.extend(g for g in gdf.geometry if g is not None and not g.is_empty and g.geom_type in ("Polygon", "MultiPolygon")) except Exception: pass land_union = None land_prepared = None if land_polys: try: land_union = unary_union(land_polys).buffer(0) land_prepared = prep(land_union) # Also emit the union as a LANDMASK feature so the frontend can # paint land in the chart's own colour when OSM is hidden. geoms = list(land_union.geoms) if land_union.geom_type.startswith("Multi") else [land_union] for g in geoms: if g.geom_type != "Polygon" or g.is_empty: continue rings = [[list(c[:2]) for c in g.exterior.coords]] rings.extend([list(c[:2]) for c in r.coords] for r in g.interiors) out.append({ "type": "Feature", "geometry": {"type": "Polygon", "coordinates": rings}, "properties": {"layer": "LANDMASK", "bbox": _bbox_of_coords(rings)}, }) except Exception as e: log.warning("Land union failed for %s: %s", enc_path, e) # ── SOUNDG (drop sounding points that fall on land) ── try: gdf = gpd.read_file(str(enc_path), layer="SOUNDG", engine="pyogrio") if not gdf.empty: gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue geoms = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] for g in geoms: if not getattr(g, "has_z", False): continue if land_prepared is not None and land_prepared.contains(g): continue # sounding lies on land — skip out.append({ "type": "Feature", "geometry": {"type": "Point", "coordinates": [g.x, g.y]}, "properties": {"layer": "SOUNDG", "depth": round(g.z, 1), "unit": unit, "bbox": [g.x, g.y, g.x, g.y]}, }) except Exception: pass # ── DEPCNT (clip contour lines so they don't cross land) ── try: gdf = gpd.read_file(str(enc_path), layer="DEPCNT", engine="pyogrio") if not gdf.empty: gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue depth = _safe(row.get("VALDCO")) if land_union is not None and geom.intersects(land_union): try: geom = geom.difference(land_union) except Exception: pass if geom.is_empty: continue geoms = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] for g in geoms: if g.geom_type != "LineString" or g.is_empty: continue coords = [list(c[:2]) for c in g.coords] if len(coords) < 2: continue out.append({ "type": "Feature", "geometry": {"type": "LineString", "coordinates": coords}, "properties": {"layer": "DEPCNT", "depth": round(float(depth), 1) if depth is not None else None, "unit": unit, "bbox": _bbox_of_coords(coords)}, }) except Exception: pass def _emit_polygon_feature(g, props): """Write one Polygon (single-part) feature, with bbox + props.""" if g.geom_type != "Polygon" or g.is_empty: return rings = [[list(c[:2]) for c in g.exterior.coords]] rings.extend([list(c[:2]) for c in r.coords] for r in g.interiors) out.append({ "type": "Feature", "geometry": {"type": "Polygon", "coordinates": rings}, "properties": {**props, "bbox": _bbox_of_coords(rings)}, }) # ── DEPARE (clipped to remove land overlap) ── try: gdf = gpd.read_file(str(enc_path), layer="DEPARE", engine="pyogrio") if not gdf.empty: gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue dmin = _safe(row.get("DRVAL1")) dmax = _safe(row.get("DRVAL2")) if land_union is not None and geom.intersects(land_union): try: geom = geom.difference(land_union) except Exception: pass if geom.is_empty: continue props = { "layer": "DEPARE", "depth_min": round(float(dmin), 1) if dmin is not None else None, "depth_max": round(float(dmax), 1) if dmax is not None else None, "unit": unit, } geoms = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] for g in geoms: _emit_polygon_feature(g, props) except Exception: pass return out def _emit_geometry(geom, props: dict) -> list[dict]: """Flatten a shapely geometry into GeoJSON Feature dicts with bbox. Handles Point / LineString / Polygon and their Multi* variants.""" out = [] geoms = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] for g in geoms: if g.is_empty: continue if g.geom_type == "Point": coords = [g.x, g.y] gtype = "Point" bbox = [g.x, g.y, g.x, g.y] elif g.geom_type == "LineString": coords = [list(c[:2]) for c in g.coords] gtype = "LineString" bbox = _bbox_of_coords(coords) elif g.geom_type == "Polygon": rings = [[list(c[:2]) for c in g.exterior.coords]] rings += [[list(c[:2]) for c in r.coords] for r in g.interiors] coords = rings gtype = "Polygon" bbox = _bbox_of_coords(coords) else: continue out.append({ "type": "Feature", "geometry": {"type": gtype, "coordinates": coords}, "properties": {**{k: v for k, v in props.items() if v is not None}, "bbox": bbox}, }) return out # ── Terrain / land features ─────────────────────────────────────────────────── _LAND_POLY_LAYERS = {"LNDARE": "land_area", "BUAARE": "built_up_area", "BUISGL": "building"} # building footprints as polygons _LAND_LINE_LAYERS = {"COALNE": "coastline"} _LAND_POINT_LAYERS = {"BRIDGE": "bridge", "HRBFAC": "harbour", "SILTNK": "silo_tank"} def _split_coalne(geom, props: dict) -> list[dict]: """Parte un COALNE LineString/MultiLineString en segmentos descartando saltos artificiales (bordes de celda de carta que cruzan el océano). Un vértice-a-vértice de una costa real rara vez supera ~0.5° de arco (~55 km). Si el salto supera MAX_DEG se considera borde artificial y se corta ahí, emitiendo el segmento acumulado hasta ese punto. """ MAX_DEG = 0.5 # ~55 km — ajustar si aparecen cortes en costas rectas largas def _emit_seg(coords): if len(coords) < 2: return [] bbox = _bbox_of_coords(coords) return [{"type": "Feature", "geometry": {"type": "LineString", "coordinates": coords}, "properties": {**{k: v for k, v in props.items() if v is not None}, "bbox": bbox}}] geoms = list(geom.geoms) if geom.geom_type.startswith("Multi") else [geom] out = [] for g in geoms: if g.geom_type != "LineString": continue pts = [list(c[:2]) for c in g.coords] seg = [pts[0]] for i in range(1, len(pts)): dx = abs(pts[i][0] - pts[i-1][0]) dy = abs(pts[i][1] - pts[i-1][1]) if dx > MAX_DEG or dy > MAX_DEG: # salto artificial → emitir segmento acumulado y empezar uno nuevo out.extend(_emit_seg(seg)) seg = [pts[i]] else: seg.append(pts[i]) out.extend(_emit_seg(seg)) return out def _parse_land(enc_path: Path) -> list[dict]: if not _HAS_GPD: return [] out = [] for layer_name, category in _LAND_POLY_LAYERS.items(): try: gdf = gpd.read_file(str(enc_path), layer=layer_name, engine="pyogrio") except Exception: continue if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue props = {"layer": layer_name, "category": category, "name": _safe(row.get("OBJNAM"))} out.extend(_emit_geometry(geom, props)) for layer_name, category in _LAND_LINE_LAYERS.items(): try: gdf = gpd.read_file(str(enc_path), layer=layer_name, engine="pyogrio") except Exception: continue if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue props = {"layer": layer_name, "category": category, "name": _safe(row.get("OBJNAM"))} # COALNE Polygon → ignorar: LANDMASK ya cubre el área de tierra. # COALNE LineString → partir en segmentos eliminando saltos largos (bordes # artificiales de celda que crean líneas diagonales cruzando el océano). if geom.geom_type in ("Polygon", "MultiPolygon"): continue # artefacto: COALNE área → ya cubierto por LANDMASK out.extend(_split_coalne(geom, props)) for layer_name, category in _LAND_POINT_LAYERS.items(): try: gdf = gpd.read_file(str(enc_path), layer=layer_name, engine="pyogrio") except Exception: continue if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue pt = geom.centroid if geom.geom_type != "Point" else geom props = {"layer": layer_name, "category": category, "name": _safe(row.get("OBJNAM")), "height_m": _safe(row.get("HEIGHT"))} if layer_name == "BRIDGE": props["horclr"] = _safe(row.get("HORCLR")) props["verclr"] = _safe(row.get("VERCLR")) out.append({"type": "Feature", "geometry": {"type": "Point", "coordinates": [pt.x, pt.y]}, "properties": {k: v for k, v in props.items() if v is not None}}) return out # ── Navigational hazards ────────────────────────────────────────────────────── _HAZARD_LAYERS = {"WRECKS": "wreck", "OBSTRN": "obstruction", "UWTROC": "rock"} def _parse_hazards(enc_path: Path) -> list[dict]: if not _HAS_GPD: return [] out = [] for layer_name, category in _HAZARD_LAYERS.items(): try: gdf = gpd.read_file(str(enc_path), layer=layer_name, engine="pyogrio") except Exception: continue if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue pt = geom.centroid if geom.geom_type != "Point" else geom props = {"layer": layer_name, "category": category, "name": _safe(row.get("OBJNAM")), "depth": _safe(row.get("VALSOU"))} if layer_name == "WRECKS": props["catwrk"] = _safe(row.get("CATWRK")) if layer_name == "UWTROC": props["catuwr"] = _safe(row.get("CATUWR")) if layer_name == "OBSTRN": props["catobs"] = _safe(row.get("CATOBS")) out.append({"type": "Feature", "geometry": {"type": "Point", "coordinates": [pt.x, pt.y]}, "properties": {k: v for k, v in props.items() if v is not None}}) return out # ── Navigation zones / areas ────────────────────────────────────────────────── _ZONE_LAYERS = { "RESARE": "restricted", "CTNARE": "caution", "ACHARE": "anchorage", "TSSLPT": "traffic_lane","PRCARE": "precautionary","FAIRWY": "fairway", "DMPGRD": "dumpground", "PIPARE": "pipeline_area", } def _parse_zones(enc_path: Path) -> list[dict]: if not _HAS_GPD: return [] out = [] for layer_name, category in _ZONE_LAYERS.items(): try: gdf = gpd.read_file(str(enc_path), layer=layer_name, engine="pyogrio") except Exception: continue if gdf.empty: continue gdf = gdf.to_crs(epsg=4326) for _, row in gdf.iterrows(): geom = row.geometry if geom is None or geom.is_empty: continue props = {"layer": layer_name, "category": category, "name": _safe(row.get("OBJNAM")), "info": _safe(row.get("INFORM"))} if layer_name == "RESARE": props["catrea"] = _safe(row.get("CATREA")) props["restrn"] = _safe_list(row.get("RESTRN")) if layer_name == "ACHARE": props["catach"] = _safe(row.get("CATACH")) out.extend(_emit_geometry(geom, props)) return out def _ensure_meta(cell_id: str): """Write meta.json with auto-detected region if not already set.""" meta = get_meta(cell_id) if "region" not in meta: set_meta(cell_id, region=auto_region(cell_id)) def install_from_zip(zip_path: Path) -> list[str]: if not _HAS_GPD: raise RuntimeError("geopandas not installed — cannot parse .000 ENC files. " "Run: pip install geopandas pyogrio") installed = [] with zipfile.ZipFile(zip_path) as z: enc_files = [n for n in z.namelist() if n.upper().endswith(".000")] if not enc_files: raise ValueError("No .000 ENC files found in ZIP") for enc_name in enc_files: cell_id = Path(enc_name).stem.upper() cell_dir = CHARTS_DIR / cell_id cell_dir.mkdir(exist_ok=True) enc_dest = cell_dir / f"{cell_id}.000" with z.open(enc_name) as src, open(enc_dest, "wb") as dst: shutil.copyfileobj(src, dst) _build_cache(cell_id, enc_dest) _ensure_meta(cell_id) installed.append(cell_id) return installed def _read_cell_id(enc_path: Path) -> str: """Lee el nombre de celda real del header DSID del S-57. Si el archivo tiene nombre temporal (TMP*), usa DSID_DSNM o DSID_DSNM como fallback.""" stem = enc_path.stem.upper() # Si el nombre ya parece un cell ID válido (no empieza con TMP), usarlo directamente if not stem.startswith("TMP"): return stem # Intentar leer DSNM del header S-57 via GDAL try: from osgeo import gdal ds = gdal.OpenEx(str(enc_path)) if ds: md = ds.GetMetadata() or {} dsnm = (md.get("DSID_DSNM") or "").strip().upper() if dsnm and not dsnm.startswith("TMP"): return dsnm except Exception as e: log.debug("Could not read DSNM from %s: %s", enc_path, e) # Fallback: usar el stem aunque sea TMP (se mostrará en lista) return stem def install_from_enc(enc_path: Path, orig_name: str | None = None) -> str: if not _HAS_GPD: raise RuntimeError("geopandas not installed — cannot parse .000 ENC files. " "Run: pip install geopandas pyogrio") # Preferir el nombre original del archivo (sin extensión) como cell_id. # Si no viene o empieza con TMP, intentar leerlo del header S-57. if orig_name and not orig_name.upper().startswith("TMP"): cell_id = orig_name.upper() else: cell_id = _read_cell_id(enc_path) cell_dir = CHARTS_DIR / cell_id cell_dir.mkdir(exist_ok=True) enc_dest = cell_dir / f"{cell_id}.000" shutil.copy2(enc_path, enc_dest) _build_cache(cell_id, enc_dest) _ensure_meta(cell_id) return cell_id def _build_cache(cell_id: str, enc_path: Path): log.info("Parsing ENC %s …", cell_id) cell_dir = CHARTS_DIR / cell_id features = _parse_cell(enc_path) with open(cell_dir / "features.geojson", "w") as f: json.dump({"type": "FeatureCollection", "features": features}, f) log.info("ENC %s → %d aid features", cell_id, len(features)) depths = _parse_depths(enc_path) with open(cell_dir / "depths.geojson", "w") as f: json.dump({"type": "FeatureCollection", "features": depths}, f) log.info("ENC %s → %d depth features", cell_id, len(depths)) land = _parse_land(enc_path) land_path = cell_dir / "land.geojson" # Preserve existing land cache when the .000 has no LNDARE layer # (e.g. custom charts built from CSV that only contain nav aids) if land or not land_path.exists(): with open(land_path, "w") as f: json.dump({"type": "FeatureCollection", "features": land}, f) log.info("ENC %s → %d land features%s", cell_id, len(land), " (preserved existing)" if not land and land_path.exists() else "") hazards = _parse_hazards(enc_path) hazards_path = cell_dir / "hazards.geojson" if hazards or not hazards_path.exists(): with open(hazards_path, "w") as f: json.dump({"type": "FeatureCollection", "features": hazards}, f) log.info("ENC %s → %d hazard features", cell_id, len(hazards)) zones = _parse_zones(enc_path) zones_path = cell_dir / "zones.geojson" if zones or not zones_path.exists(): with open(zones_path, "w") as f: json.dump({"type": "FeatureCollection", "features": zones}, f) log.info("ENC %s → %d zone features", cell_id, len(zones)) # Cache count and bbox in meta.json so list_cells() doesn't need to read features.geojson # Solo usa features de tipo Point para el bbox — los LDLINE (LineString) tienen # coordinates[0] = [lon, lat] (lista), no un float, y rompen min(lons). bbox = None try: lons, lats = [], [] for f in features: geom = f.get("geometry", {}) if geom.get("type") == "Point": coords = geom.get("coordinates", []) if len(coords) >= 2: lons.append(float(coords[0])) lats.append(float(coords[1])) if lons: bbox = [min(lons), min(lats), max(lons), max(lats)] except Exception: pass meta = get_meta(cell_id) meta["feature_count"] = len(features) meta["bbox"] = bbox _meta_path(cell_dir).write_text(json.dumps(meta)) def list_cells() -> list[dict]: cells = [] for cell_dir in sorted(CHARTS_DIR.iterdir()): if not cell_dir.is_dir(): continue # Skip leftover temp directories from failed installs if cell_dir.name.startswith("TMP"): continue enc = list(cell_dir.glob("*.000")) cache = cell_dir / "features.geojson" meta = get_meta(cell_dir.name) count = meta.get("feature_count") bbox = meta.get("bbox") # Recompute bbox si es None (celdas instaladas con código anterior que tenía # el bug LDLINE: coordinates[0] de un LineString es [lon,lat] no float → # min(lons) fallaba → bbox quedaba None → celdas excluidas del viewport). # Solo usa features Point para el bbox y persiste en meta.json para la próxima vez. if bbox is None and cache.exists(): try: fc = json.loads(cache.read_text()) feats = fc.get("features") or [] if count is None: count = len(feats) lons, lats = [], [] for f in feats: geom = f.get("geometry", {}) if geom.get("type") == "Point": coords = geom.get("coordinates", []) if len(coords) >= 2: lons.append(float(coords[0])) lats.append(float(coords[1])) if lons: bbox = [min(lons), min(lats), max(lons), max(lats)] # Persistir para que las siguientes llamadas sean rápidas meta["bbox"] = bbox if count is not None: meta["feature_count"] = count try: _meta_path(cell_dir).write_text(json.dumps(meta)) except Exception: pass except Exception: if count is None: count = 0 if count is None: count = 0 cells.append({ "id": cell_dir.name, "enc_file": enc[0].name if enc else None, "features": count or 0, "cached": cache.exists(), "bbox": bbox, "region": get_region(cell_dir.name), }) return cells def delete_cell(cell_id: str): cell_dir = CHARTS_DIR / cell_id.upper() if cell_dir.exists(): shutil.rmtree(cell_dir) def scan_and_install(directory: str | Path) -> dict: """ Walk *directory* (local path — e.g. an SD card drive letter like E:\\) and install every .000 ENC cell and every .zip archive found. Returns: { "installed": ["CELL1", "CELL2", ...], # successfully installed "skipped": ["already_there.000"], # already installed (skipped) "errors": [{"file": "...", "error": "..."}] # parse failures } """ base = Path(directory) if not base.exists(): raise FileNotFoundError(f"Path not found: {base}") if not base.is_dir(): raise NotADirectoryError(f"Not a directory: {base}") installed: list[str] = [] skipped: list[str] = [] errors: list[dict] = [] # Collect candidates — full recursive walk of the directory tree. # Typical SD card layouts: ENC_ROOT/*.000, ENC_ROOT/US5NY01M/*.000, # NOAA_ENC/US5NY01M/US5NY01M.000, etc. candidates: list[Path] = [] for entry in base.rglob("*"): if entry.is_file(): candidates.append(entry) # Process .zip first (may contain multiple cells), then .000 zips = [c for c in candidates if c.suffix.lower() == ".zip"] encs = [c for c in candidates if c.suffix.upper() == ".000"] existing_ids = {cell_dir.name.upper() for cell_dir in CHARTS_DIR.iterdir() if cell_dir.is_dir()} for zp in zips: try: import zipfile as _zf with _zf.ZipFile(zp) as z: names = z.namelist() has_csv = any(n.lower().endswith(".csv") for n in names) has_enc = any(n.upper().endswith(".000") for n in names) if has_csv and not has_enc: ids = install_from_csv_zip(zp) else: ids = install_from_zip(zp) installed.extend(ids) except Exception as exc: errors.append({"file": zp.name, "error": str(exc)}) for ep in encs: cid = ep.stem.upper() if cid in existing_ids: skipped.append(ep.name) continue try: cid = install_from_enc(ep, ep.stem.upper()) installed.append(cid) except Exception as exc: errors.append({"file": ep.name, "error": str(exc)}) return {"installed": installed, "skipped": skipped, "errors": errors} def get_all_features() -> dict: all_features = [] for cell_dir in CHARTS_DIR.iterdir(): cache = cell_dir / "features.geojson" if not cache.exists(): continue try: fc = json.loads(cache.read_text()) except Exception: continue cell_id = cell_dir.name region = get_region(cell_id) for f in (fc.get("features") or []): p = f.setdefault("properties", {}) p["cell"] = cell_id p["cell_region"] = region # Backfill aid_type for old caches that pre-date the classifier. if "aid_type" not in p: p["aid_type"] = classify(p.get("layer", ""), p) all_features.extend(fc["features"]) return {"type": "FeatureCollection", "features": all_features} def _aggregate_cache(filename: str, bbox=None) -> dict: """Generic aggregator: read from every installed cell.""" all_features = [] if bbox is not None: w, s, e, n = bbox for cell_dir in CHARTS_DIR.iterdir(): cache = cell_dir / filename if not cache.exists(): continue try: fc = json.loads(cache.read_text()) except Exception: continue for f in (fc.get("features") or []): p = f.setdefault("properties", {}) p["cell"] = cell_dir.name if bbox is not None: fb = p.get("bbox") if fb and (fb[2] < w or fb[0] > e or fb[3] < s or fb[1] > n): continue all_features.append(f) return {"type": "FeatureCollection", "features": all_features} def get_all_depths(bbox: tuple[float, float, float, float] | None = None) -> dict: return _aggregate_cache("depths.geojson", bbox) def get_all_land(bbox: tuple[float, float, float, float] | None = None) -> dict: return _aggregate_cache("land.geojson", bbox) def get_all_hazards() -> dict: return _aggregate_cache("hazards.geojson") def get_all_zones(bbox: tuple[float, float, float, float] | None = None) -> dict: return _aggregate_cache("zones.geojson", bbox)