""" S-57 ENC writer — produces structurally valid `.000` files that the GDAL S-57 driver can read. Reuses a NOAA DDR template verbatim (the schema is fixed by IHO S-57 ed. 3.1) and only encodes data records. ISO 8211 record format (data records, leader entry-map "3404"): Leader (24 bytes) + Directory (11 bytes per entry + FT) + Field area. Directory entry layout: tag(4) + length(3) + position(4). """ from __future__ import annotations import struct from pathlib import Path # ── ISO 8211 control characters ──────────────────────────────────────────────── FT = b"\x1e" # Field terminator UT = b"\x1f" # Unit terminator (subfield delimiter) DDR_TEMPLATE = Path(__file__).parent / "noaa_ddr_template.bin" # S-57 record names (RCNM) RCNM_VI = 110 # Isolated Node (point primitive) RCNM_VC = 120 # Connected Node RCNM_VE = 130 # Edge (line primitive) RCNM_VF = 140 # Face RCNM_FE = 100 # Feature # Object class codes — from s57objectclasses.csv (IHO S-57 Ed.3.1, Google Earth GDAL copy) # fmt: off OBJL_ACHBRT = 3 # Anchor berth OBJL_ACHARE = 4 # Anchorage area OBJL_BCNCAR = 5 # Beacon, cardinal OBJL_BCNISD = 6 # Beacon, isolated danger OBJL_BCNLAT = 7 # Beacon, lateral OBJL_BCNSAW = 8 # Beacon, safe water OBJL_BCNSPP = 9 # Beacon, special purpose/general OBJL_BERTHS = 10 # Berth OBJL_BRIDGE = 11 # Bridge OBJL_BUISGL = 12 # Building, single OBJL_BUAARE = 13 # Built-up area OBJL_BOYCAR = 14 # Buoy, cardinal OBJL_BOYISD = 16 # Buoy, isolated danger OBJL_BOYLAT = 17 # Buoy, lateral OBJL_BOYSAW = 18 # Buoy, safe water ← was 19 (BUG fixed) OBJL_BOYSPP = 19 # Buoy, special purpose/general (S-57 std acronym) OBJL_CBLOHD = 21 # Cable, overhead OBJL_CBLSUB = 22 # Cable, submarine OBJL_CANALS = 23 # Canal OBJL_CTSARE = 25 # Cargo transshipment area OBJL_COALNE = 30 # Coastline OBJL_CRANES = 35 # Crane OBJL_DAYMAR = 39 # Daymark OBJL_DWRTCL = 40 # Deep water route centre line OBJL_DWRTPT = 41 # Deep water route part OBJL_DEPARE = 42 # Depth area OBJL_DEPCNT = 43 # Depth contour OBJL_DOCARE = 45 # Dock area OBJL_DRGARE = 46 # Dredged area OBJL_DRYDOC = 47 # Dry dock OBJL_DMPGRD = 48 # Dumping ground OBJL_EXEZNE = 50 # Exclusive Economic Zone OBJL_FAIRWY = 51 # Fairway OBJL_FERYRT = 53 # Ferry route OBJL_FSHZNE = 54 # Fishery zone OBJL_FOGSIG = 58 # Fog signal OBJL_FRPARE = 60 # Free port area OBJL_GATCON = 61 # Gate / sluice OBJL_HRBARE = 63 # Harbour area (administrative) OBJL_HRBFAC = 64 # Harbour facility OBJL_HULKES = 65 # Hulk OBJL_ISTZNE = 68 # Inshore traffic zone OBJL_LAKARE = 69 # Lake OBJL_LNDARE = 71 # Land area OBJL_LNDRGN = 73 # Land region OBJL_LNDMRK = 74 # Landmark OBJL_LIGHTS = 75 # Light OBJL_LITFLT = 76 # Light float OBJL_LITVES = 77 # Light vessel OBJL_LOKBSN = 79 # Lock basin OBJL_MAGVAR = 81 # Magnetic variation OBJL_MARCUL = 82 # Marine farm/culture OBJL_MIPARE = 83 # Military practice area OBJL_MORFAC = 84 # Mooring/warping facility OBJL_NAVLNE = 85 # Navigation line OBJL_OBSTRN = 86 # Obstruction OBJL_OFSPLF = 87 # Offshore platform OBJL_OSPARE = 88 # Offshore production area OBJL_PILPNT = 90 # Pile OBJL_PILBOP = 91 # Pilot boarding place OBJL_PIPARE = 92 # Pipeline area OBJL_PIPOHD = 93 # Pipeline, overhead OBJL_PIPSOL = 94 # Pipeline, submarine/on land OBJL_PONTON = 95 # Pontoon OBJL_PRCARE = 96 # Precautionary area OBJL_PYLONS = 98 # Pylon/bridge support OBJL_RADRFL = 101 # Radar reflector OBJL_RADSTA = 102 # Radar station OBJL_RDOCAL = 104 # Radio calling-in point OBJL_RDOSTA = 105 # Radio station OBJL_RECTRC = 109 # Recommended track OBJL_RESARE = 112 # Restricted area OBJL_RETRFL = 113 # Retro-reflector OBJL_RIVERS = 114 # River OBJL_RIVBNK = 115 # River bank OBJL_ROADWY = 116 # Road OBJL_RUNWAY = 117 # Runway OBJL_SEAARE = 119 # Sea area / named water area OBJL_SPLARE = 120 # Sea-plane landing area (SPLARE) OBJL_SBDARE = 121 # Seabed area OBJL_SLCONS = 122 # Shoreline construction OBJL_SLOTOP = 126 # Slope topline OBJL_SLOGRD = 127 # Sloping ground OBJL_SOUNDG = 129 # Sounding OBJL_STSLNE = 132 # Straight territorial sea baseline OBJL_SUBTLN = 133 # Submarine transit lane OBJL_SWPARE = 134 # Swept area OBJL_TESARE = 135 # Territorial sea area OBJL_TIDEWY = 143 # Tideway OBJL_TOPMAR = 144 # Top mark OBJL_TSELNE = 145 # Traffic Separation Line OBJL_TSSBND = 146 # Traffic Separation Scheme Boundary OBJL_TSSCRS = 147 # Traffic Separation Scheme Crossing OBJL_TSSLPT = 148 # Traffic Separation Scheme Lane part OBJL_TSSRON = 149 # Traffic Separation Scheme Roundabout OBJL_TUNNEL = 151 # Tunnel OBJL_TWRTPT = 152 # Two-way route part OBJL_UWTROC = 153 # Underwater rock / awash rock OBJL_VEGATN = 155 # Vegetation OBJL_WATTUR = 156 # Water turbulence OBJL_WEDKLP = 158 # Weed/Kelp OBJL_WRECKS = 159 # Wreck OBJL_M_ACCY = 300 # Accuracy of data (meta) OBJL_M_CSCL = 301 # Compilation scale (meta) OBJL_M_COVR = 302 # Coverage (meta) OBJL_M_HDAT = 303 # Horizontal datum (meta) OBJL_M_NSYS = 306 # Navigational system of marks (meta) OBJL_M_QUAL = 308 # Quality of data (meta) OBJL_M_SDAT = 309 # Sounding datum (meta) OBJL_M_UNIT = 311 # Units of measurement (meta) # fmt: on # Lookup: S-57 acronym → OBJL code (complete IHO S-57 Ed 3.1 set) OBJL_BY_ACRONYM: dict[str, int] = { "ACHBRT": OBJL_ACHBRT, "ACHARE": OBJL_ACHARE, "BCNCAR": OBJL_BCNCAR, "BCNISD": OBJL_BCNISD, "BCNLAT": OBJL_BCNLAT, "BCNSAW": OBJL_BCNSAW, "BCNSPP": OBJL_BCNSPP, "BERTHS": OBJL_BERTHS, "BRIDGE": OBJL_BRIDGE, "BUISGL": OBJL_BUISGL, "BUAARE": OBJL_BUAARE, "BOYCAR": OBJL_BOYCAR, "BOYISD": OBJL_BOYISD, "BOYLAT": OBJL_BOYLAT, "BOYSAW": OBJL_BOYSAW, "BOYSPP": OBJL_BOYSPP, "BOYSPEC": OBJL_BOYSPP, # user alias → BOYSPP "CBLOHD": OBJL_CBLOHD, "CBLSUB": OBJL_CBLSUB, "CANALS": OBJL_CANALS, "CTSARE": OBJL_CTSARE, "COALNE": OBJL_COALNE, "CRANES": OBJL_CRANES, "DAYMAR": OBJL_DAYMAR, "DWRTCL": OBJL_DWRTCL, "DWRTPT": OBJL_DWRTPT, "DEPARE": OBJL_DEPARE, "DEPCNT": OBJL_DEPCNT, "DOCARE": OBJL_DOCARE, "DRGARE": OBJL_DRGARE, "DRYDOC": OBJL_DRYDOC, "DMPGRD": OBJL_DMPGRD, "EXEZNE": OBJL_EXEZNE, "FAIRWY": OBJL_FAIRWY, "FERYRT": OBJL_FERYRT, "FSHZNE": OBJL_FSHZNE, "FRPARE": OBJL_FRPARE, "FOGSIG": OBJL_FOGSIG, "GATCON": OBJL_GATCON, "HRBARE": OBJL_HRBARE, "HRBFAC": OBJL_HRBFAC, "HULKES": OBJL_HULKES, "ISTZNE": OBJL_ISTZNE, "LAKARE": OBJL_LAKARE, "LNDARE": OBJL_LNDARE, "LNDRGN": OBJL_LNDRGN, "LNDMRK": OBJL_LNDMRK, "LIGHTS": OBJL_LIGHTS, "LITFLT": OBJL_LITFLT, "LITVES": OBJL_LITVES, "LOKBSN": OBJL_LOKBSN, "MAGVAR": OBJL_MAGVAR, "MARCUL": OBJL_MARCUL, "MIPARE": OBJL_MIPARE, "MORFAC": OBJL_MORFAC, "NAVLNE": OBJL_NAVLNE, "OBSTRN": OBJL_OBSTRN, "OFSPLF": OBJL_OFSPLF, "OSPARE": OBJL_OSPARE, "PILPNT": OBJL_PILPNT, "PILBOP": OBJL_PILBOP, "PIPARE": OBJL_PIPARE, "PIPOHD": OBJL_PIPOHD, "PIPSOL": OBJL_PIPSOL, "PONTON": OBJL_PONTON, "PRCARE": OBJL_PRCARE, "PYLONS": OBJL_PYLONS, "RADRFL": OBJL_RADRFL, "RADSTA": OBJL_RADSTA, "RDOCAL": OBJL_RDOCAL, "RDOSTA": OBJL_RDOSTA, "RECTRC": OBJL_RECTRC, "RESARE": OBJL_RESARE, "RETRFL": OBJL_RETRFL, "RIVERS": OBJL_RIVERS, "RIVBNK": OBJL_RIVBNK, "ROADWY": OBJL_ROADWY, "RUNWAY": OBJL_RUNWAY, "SEAARE": OBJL_SEAARE, "SPLARE": OBJL_SPLARE, "SBDARE": OBJL_SBDARE, "SLCONS": OBJL_SLCONS, "SLOTOP": OBJL_SLOTOP, "SLOGRD": OBJL_SLOGRD, "SOUNDG": OBJL_SOUNDG, "STSLNE": OBJL_STSLNE, "SUBTLN": OBJL_SUBTLN, "SWPARE": OBJL_SWPARE, "TESARE": OBJL_TESARE, "TIDEWY": OBJL_TIDEWY, "TOPMAR": OBJL_TOPMAR, "TSELNE": OBJL_TSELNE, "TSSBND": OBJL_TSSBND, "TSSCRS": OBJL_TSSCRS, "TSSLPT": OBJL_TSSLPT, "TSSRON": OBJL_TSSRON, "TUNNEL": OBJL_TUNNEL, "TWRTPT": OBJL_TWRTPT, "UWTROC": OBJL_UWTROC, "VEGATN": OBJL_VEGATN, "WATTUR": OBJL_WATTUR, "WEDKLP": OBJL_WEDKLP, "WRECKS": OBJL_WRECKS, "M_ACCY": OBJL_M_ACCY, "M_CSCL": OBJL_M_CSCL, "M_COVR": OBJL_M_COVR, "M_HDAT": OBJL_M_HDAT, "M_NSYS": OBJL_M_NSYS, "M_QUAL": OBJL_M_QUAL, "M_SDAT": OBJL_M_SDAT, "M_UNIT": OBJL_M_UNIT, } # PRIM codes for FRID PRIM_POINT = 1 PRIM_LINE = 2 PRIM_AREA = 3 # ── S-57 attribute codes (ATTL) — from s57attributes.csv (IHO S-57 Ed.3.1) ─── # fmt: off ATTL_AGENCY = 1 # AGENCY — agency responsible for production ATTL_BCNSHP = 2 # BCNSHP — beacon shape ATTL_BOYSHP = 4 # BOYSHP — buoy shape ATTL_BURDEP = 5 # BURDEP — buried depth ATTL_CATCAM = 13 # CATCAM — category of cardinal mark ATTL_CATCOV = 18 # CATCOV — category of coverage (1=coverage, 2=no coverage) ATTL_CATFOG = 27 # CATFOG — category of fog signal ATTL_CATLAM = 36 # CATLAM — category of lateral mark (1=port, 2=starboard) ATTL_CATLMK = 35 # CATLMK — category of landmark ATTL_CATLIT = 37 # CATLIT — category of light ATTL_CATOBS = 42 # CATOBS — category of obstruction ATTL_CATREA = 54 # CATREA — category of restricted area ATTL_CATSIW = 63 # CATSIW — category of signal station warning ATTL_COLOUR = 75 # COLOUR — colour (list) ATTL_COLPAT = 76 # COLPAT — colour pattern ATTL_CONDTN = 81 # CONDTN — condition ATTL_DRVAL1 = 90 # DRVAL1 — draft value 1 ATTL_DRVAL2 = 91 # DRVAL2 — draft value 2 ATTL_HEIGHT = 95 # HEIGHT — height above sea surface ATTL_LITCHR = 107 # LITCHR — light characteristic ATTL_MLTYLT = 110 # MLTYLT — multiplicity of lights ATTL_NATCON = 112 # NATCON — nature of construction ATTL_OBJNAM = 116 # OBJNAM — object name (free text label) ATTL_ORIENT = 117 # ORIENT — orientation (degrees) ATTL_PEREND = 119 # PEREND — period end (season) ATTL_PERSTA = 120 # PERSTA — period start (season) ATTL_QUASOU = 126 # QUASOU — quality of sounding measurement ATTL_SIGGRP = 141 # SIGGRP — signal group ("(2)", etc.) ATTL_SIGPER = 142 # SIGPER — signal period (seconds) ATTL_STATUS = 149 # STATUS — status (permanent, occasional…) ATTL_VALSOU = 179 # VALSOU — value of sounding (metres) ATTL_VERACC = 180 # VERACC — vertical accuracy ATTL_VERCLR = 181 # VERCLR — vertical clearance ATTL_VALNMR = 178 # VALNMR — value of nominal range (nautical miles) ATTL_VERDAT = 187 # VERDAT — vertical datum # fmt: on # Convenience dict: acronym → ATTL code (superset of the named constants above) ATTR_CODE: dict[str, int] = { "AGENCY": 1, "BCNSHP": 2, "BUISHP": 3, "BOYSHP": 4, "BURDEP": 5, "CALSGN": 6, "CATAIR": 7, "CATACH": 8, "CATBRG": 9, "CATBUA": 10, "CATCBL": 11, "CATCAN": 12, "CATCAM": 13, "CATCHP": 14, "CATCOA": 15, "CATCTR": 16, "CATCON": 17, "CATCOV": 18, "CATCRN": 19, "CATDAM": 20, "CATDIS": 21, "CATDOC": 22, "CATDPG": 23, "CATFNC": 24, "CATFRY": 25, "CATFIF": 26, "CATFOG": 27, "CATFOR": 28, "CATGAT": 29, "CATHAF": 30, "CATHLK": 31, "CATICE": 32, "CATINB": 33, "CATLND": 34, "CATLMK": 35, "CATLAM": 36, "CATLIT": 37, "CATMFA": 38, "CATMPA": 39, "CATMOR": 40, "CATNAV": 41, "CATOBS": 42, "CATOFP": 43, "CATOLB": 44, "CATPLE": 45, "CATPIL": 46, "CATPIP": 47, "CATPRA": 48, "CATPYL": 49, "CATQUA": 50, "CATRAS": 51, "CATRTB": 52, "CATROS": 53, "CATREA": 54, "CATSEA": 55, "CATSIT": 56, "CATSLC": 57, "CATSPM": 58, "CATSCF": 59, "CATSUB": 60, "CAATTS": 61, "CATTSS": 62, "CATSIW": 63, "CATTRK": 64, "CATVEG": 65, "CATWED": 66, "CATWRK": 67, "CATZOC": 68, "CATWAT": 69, "COLOUR": 75, "COLPAT": 76, "CONDTN": 81, "CONRAD": 82, "CONVIS": 83, "CURVEL": 84, "DATEND": 85, "DATSTA": 86, "DRVAL1": 90, "DRVAL2": 91, "ELEVAT": 93, "ESTRNG": 94, "HEIGHT": 95, "HORACC": 96, "HORCLR": 97, "HORLEN": 98, "HORWID": 99, "ICEFAC": 100,"INFORM": 101,"JRSDTN": 102,"LIFCAP": 103, "LITCHR": 107, "LITVIS": 108,"MARSYS": 109,"MLTYLT": 110,"NATION": 111, "NATCON": 112, "NATQUA": 113,"NATSUR": 114,"NOBJNM": 115,"OBJNAM": 116, "ORIENT": 117, "PEREND": 119,"PERSTA": 120,"PICREP": 121,"POSACC": 122, "PRCTRY": 124, "PRODCT": 125,"QUASOU": 126,"RADWAL": 127,"RADIUS": 128, "RYRMGV": 133, "SCAMAX": 134,"SCAMIN": 135,"SCVAL1": 136,"SCVAL2": 137, "SECTR1": 138, "SECTR2": 139,"SHIPAM": 140,"SIGGRP": 141,"SIGPER": 142, "SIGSEQ": 143, "SOUACC": 144,"SDISMX": 145,"SDISMN": 146,"SORDAT": 147, "SORIND": 148, "STATUS": 149,"SURATH": 150,"SUREND": 151,"SURSTA": 152, "SURTYP": 153, "TECSOU": 156,"TXTDSC": 157,"TS_TSP": 158,"TS_TSV": 159, "T_ACWL": 160, "T_HWLW": 161,"T_MTOD": 162,"T_THDF": 163,"T_TIMS": 164, "T_TRNP": 165, "T_VAHF": 166,"T_VAVL": 167,"TIMEND": 168,"TIMSTA": 169, "TOPSHP": 170, "TRAFIC": 171,"VALDCO": 172,"VERACC": 180,"VERCLR": 181, "VERCCL": 182, "VERCOP": 183,"VERCSA": 184,"VERDAT": 187,"VERLEN": 188, "WATLEV": 187, "VALSOU": 179, "VALNMR": 178, } # ── Binary primitives (little-endian, S-57 convention) ───────────────────────── def b11(n): return int(n).to_bytes(1, "little", signed=False) def b12(n): return int(n).to_bytes(2, "little", signed=False) def b14(n): return int(n).to_bytes(4, "little", signed=False) def b21(n): return int(n).to_bytes(1, "little", signed=True) def b22(n): return int(n).to_bytes(2, "little", signed=True) def b24(n): return int(n).to_bytes(4, "little", signed=True) def A_var(s: str) -> bytes: """Variable-length ASCII subfield, delimited by UT.""" return s.encode("latin-1") + UT def A_fixed(s: str, n: int) -> bytes: """Fixed-width ASCII (no UT).""" return s.encode("latin-1").ljust(n, b" ")[:n] def R_fixed(value: float, n: int = 4) -> bytes: """Fixed-width ASCII real, e.g. R(4) = '03.1'.""" s = f"{value:.{max(0, n-2)}f}" return s.encode("latin-1").ljust(n, b" ")[:n] def name_5(rcnm: int, rcid: int) -> bytes: """S-57 NAME field: RCNM (1 byte) + RCID (4 bytes LE) = 5 bytes (= 40 bits).""" return b11(rcnm) + b14(rcid) # ── Field builders (one per S-57 field tag) ──────────────────────────────────── def field_0001(rcid: int) -> bytes: """Record header (b12 RCID) + FT.""" return b12(rcid) + FT def field_DSID(*, rcnm=10, rcid=1, expp=1, intu=3, dsnm="", edtn="1", updn="0", uadt="", isdt="", sted=3.1, prsp=1, psdn="", pred="2.0", prof=1, agen=999, comt="") -> bytes: return (b11(rcnm) + b14(rcid) + b11(expp) + b11(intu) + A_var(dsnm) + A_var(edtn) + A_var(updn) + A_fixed(uadt, 8) + A_fixed(isdt, 8) + R_fixed(sted, 4) + b11(prsp) + A_var(psdn) + A_var(pred) + b11(prof) + b12(agen) + A_var(comt) + FT) def field_DSSI(*, dstr=2, aall=1, nall=1, nomr=0, nocr=0, nogr=0, nolr=0, noin=0, nocn=0, noed=0, nofa=0) -> bytes: """Format: (3b11, 8b14)""" return (b11(dstr) + b11(aall) + b11(nall) + b14(nomr) + b14(nocr) + b14(nogr) + b14(nolr) + b14(noin) + b14(nocn) + b14(noed) + b14(nofa) + FT) def field_DSPM(*, rcnm=20, rcid=1, hdat=2, vdat=17, sdat=23, cscl=50000, duni=1, huni=1, puni=1, coun=1, comf=10000000, somf=10, comt="") -> bytes: """Format: (b11, b14, 3b11, b14, 4b11, 2b14, A) HDAT=2 (WGS84), VDAT=17 (MLLW), SDAT=23 (MLLW), DUNI=1 (m), HUNI=1 (m), PUNI=1 (m), COUN=1 (lat/long) COMF=10^7 means coords stored as int = (degrees * 10^7). SOMF=10 means depths stored as int = (metres * 10).""" return (b11(rcnm) + b14(rcid) + b11(hdat) + b11(vdat) + b11(sdat) + b14(cscl) + b11(duni) + b11(huni) + b11(puni) + b11(coun) + b14(comf) + b14(somf) + A_var(comt) + FT) def field_VRID(*, rcnm: int, rcid: int, rver: int = 1, ruin: int = 1) -> bytes: """Format: (b11, b14, b12, b11). RUIN=1 means insert.""" return b11(rcnm) + b14(rcid) + b12(rver) + b11(ruin) + FT def field_SG2D(coords_deg: list[tuple[float, float]], comf: int) -> bytes: """Format: *(b24 YCOO, b24 XCOO). Coords are (lon, lat) in degrees → ints multiplied by COMF (typically 10^7). Note S-57 stores YCOO (lat) before XCOO (lon) for each pair.""" out = b"" for lon, lat in coords_deg: y_int = int(round(lat * comf)) x_int = int(round(lon * comf)) out += b24(y_int) + b24(x_int) return out + FT def field_VRPT(pointers: list[tuple[int, int, int, int, int, int]]) -> bytes: """Format: *(B(40) NAME, b11 ORNT, b11 USAG, b11 TOPI, b11 MASK). Each item: (rcnm, rcid, ornt, usag, topi, mask).""" out = b"" for rcnm, rcid, ornt, usag, topi, mask in pointers: out += name_5(rcnm, rcid) + b11(ornt) + b11(usag) + b11(topi) + b11(mask) return out + FT def field_FRID(*, rcnm: int = RCNM_FE, rcid: int, prim: int, grup: int = 1, objl: int, rver: int = 1, ruin: int = 1) -> bytes: """Format: (b11, b14, 2b11, 2b12, b11).""" return (b11(rcnm) + b14(rcid) + b11(prim) + b11(grup) + b12(objl) + b12(rver) + b11(ruin) + FT) def field_FOID(*, agen: int = 999, fidn: int, fids: int = 1) -> bytes: """Format: (b12, b14, b12).""" return b12(agen) + b14(fidn) + b12(fids) + FT def field_ATTF(attrs: list[tuple[int, str]]) -> bytes: """Format: *(b12 ATTL, A ATVL). Empty list → no field at all (caller skips).""" out = b"" for attl, atvl in attrs: out += b12(attl) + A_var(str(atvl)) return out + FT def field_FSPT(pointers: list[tuple[int, int, int, int, int]]) -> bytes: """Format: *(B(40) NAME, b11 ORNT, b11 USAG, b11 MASK). Each: (rcnm, rcid, ornt, usag, mask). ORNT=1=forward,2=reverse,255=null; USAG=1=exterior,2=interior,3=both,255=null; MASK=1=mask,2=show,255=null.""" out = b"" for rcnm, rcid, ornt, usag, mask in pointers: out += name_5(rcnm, rcid) + b11(ornt) + b11(usag) + b11(mask) return out + FT # ── Record builder (data record with entry-map 5504) ─────────────────────────── def build_record(fields: list[tuple[str, bytes]]) -> bytes: """fields: list of (4-char tag, field bytes). Field bytes include trailing FT. Uses entry-map "5504": field-length=5 digits (max 99999 bytes per field), position=5 digits (max 99999). Each directory entry is 14 bytes: tag(4)+len(5)+pos(5). History of bugs: "3404" (original): 3-digit length = max 999 bytes → LNDARE polygons produce SG2D fields of 2000-53000 bytes → length truncated → GDAL: "Not enough byte to initialize field 'SG2D'". "4404": 4-digit = max 9999 bytes → still not enough for large coastal outlines. "5504": 5-digit = max 99999 bytes → handles any realistic polygon after RDP. The DDR template uses "3404" for its own directory, but that only applies to how the DDR is parsed. GDAL reads each data record's leader independently, so data records may use a different entry-map than the DDR — this is valid ISO 8211. """ field_area = b"" dir_str = "" pos = 0 for tag, data in fields: assert len(tag) == 4, f"tag must be 4 chars: {tag!r}" flen = len(data) assert flen <= 99999, f"field {tag} too large ({flen} bytes); apply more RDP" assert pos <= 99999, f"field area position overflow at {tag} (pos={pos})" dir_str += f"{tag:<4s}{flen:05d}{pos:05d}" field_area += data pos += flen directory = dir_str.encode("latin-1") + FT base_addr = 24 + len(directory) record_len = base_addr + len(field_area) assert record_len <= 99999, f"record too large ({record_len} bytes)" leader = ( f"{record_len:05d}" # 0-4: record length " " # 5: interchange level (space for DR) "D" # 6: leader id " " # 7: in-line code extension (space for DR) " " # 8: version " " # 9: app indicator " " # 10-11: field control length (00 for DR) f"{base_addr:05d}" # 12-16: base address of field area " " # 17-19: extended character set indicator "5504" # 20-23: entry map (len=5, pos=5, reserved=0, tag=4) ).encode("latin-1") assert len(leader) == 24, f"leader len {len(leader)}" return leader + directory + field_area # ── High-level builder for a complete `.000` ────────────────────────────────── class S57Cell: """Builds a complete S-57 cell. Call add_* methods then write().""" def __init__(self, *, dsnm: str, edition: int = 1, intu: int = 5, scale: int = 50000, agen: int = 999, comt: str = "AR ECDIS custom", issue_date: str = "20260428"): self.dsnm = dsnm self.edition = edition self.intu = intu self.scale = scale self.agen = agen self.comt = comt self.issue_date = issue_date self.comf = 10_000_000 # 10^7 — coords stored as int(deg * COMF) self.somf = 10 # depths × 10 # Spatial primitives & features accumulated by callers self._next_vi = 1 self._next_vc = 1 self._next_ve = 1 self._next_fid = 1 self._records: list[tuple[str, list[tuple[str, bytes]]]] = [] # Counters for DSSI self._n_meta = 0 self._n_pt = 0 # NOMR isolated nodes self._n_cn = 0 # NOCN connected nodes self._n_ed = 0 # NOED edges self._n_geo = 0 self._n_col = 0 # NOCR collection self._n_in = 0 # NOIN isolated nodes alt count self._n_lr = 0 # NOLR self._n_gr = 0 # NOGR self._n_fa = 0 # NOFA faces # --- Add an isolated node (for points: buoys, lights, landmarks) --- def add_isolated_node(self, lon: float, lat: float) -> tuple[int, int]: rcid = self._next_vi self._next_vi += 1 fields = [ ("0001", field_0001(rcid)), ("VRID", field_VRID(rcnm=RCNM_VI, rcid=rcid)), ("SG2D", field_SG2D([(lon, lat)], self.comf)), ] self._records.append(("VI", fields)) self._n_pt += 1 self._n_in += 1 return (RCNM_VI, rcid) # --- Add an edge connecting (start_lon, start_lat) → (end_lon, end_lat) --- # First creates two connected nodes, then the edge that points to them with # intermediate coordinates. def add_edge(self, coords: list[tuple[float, float]]) -> tuple[int, int]: if len(coords) < 2: raise ValueError("edge needs ≥2 coordinates") # Create connected nodes for the endpoints cn_start = self._add_connected_node(*coords[0]) cn_end = self._add_connected_node(*coords[-1]) # The edge: VRID + VRPT (pointers to both CNs) + SG2D (intermediate coords) rcid = self._next_ve self._next_ve += 1 # Intermediate points only (S-57 stores edges as: start_CN, intermediates..., end_CN) intermediate = coords[1:-1] fields = [ ("0001", field_0001(rcid)), ("VRID", field_VRID(rcnm=RCNM_VE, rcid=rcid)), # VRPT: two pointers — to start node, to end node ("VRPT", field_VRPT([ (RCNM_VC, cn_start[1], 1, 1, 1, 255), # ORNT=1, USAG=1=ext, TOPI=1=begin, MASK=255 (RCNM_VC, cn_end[1], 2, 1, 2, 255), # ORNT=2=reverse, TOPI=2=end ])), ] if intermediate: fields.append(("SG2D", field_SG2D(intermediate, self.comf))) self._records.append(("VE", fields)) self._n_ed += 1 return (RCNM_VE, rcid) def _add_connected_node(self, lon: float, lat: float) -> tuple[int, int]: rcid = self._next_vc self._next_vc += 1 fields = [ ("0001", field_0001(rcid)), ("VRID", field_VRID(rcnm=RCNM_VC, rcid=rcid)), ("SG2D", field_SG2D([(lon, lat)], self.comf)), ] self._records.append(("VC", fields)) self._n_cn += 1 return (RCNM_VC, rcid) # --- Feature: point --- def add_point_feature(self, *, objl: int, lon: float, lat: float, attrs: list[tuple[int, str]] | None = None) -> int: rcid = self._next_fid self._next_fid += 1 fid = rcid node_name = self.add_isolated_node(lon, lat) fields = [ ("0001", field_0001(rcid)), ("FRID", field_FRID(rcid=rcid, prim=PRIM_POINT, objl=objl)), ("FOID", field_FOID(agen=self.agen, fidn=fid)), ] if attrs: fields.append(("ATTF", field_ATTF(attrs))) fields.append(("FSPT", field_FSPT([ (node_name[0], node_name[1], 255, 255, 255), # ORNT/USAG/MASK = NULL for point ]))) self._records.append(("FE", fields)) self._n_geo += 1 return rcid # --- Feature: line (e.g. COALNE, DEPCNT) --- def add_line_feature(self, *, objl: int, coords: list[tuple[float, float]], attrs: list[tuple[int, str]] | None = None) -> int: """coords: ordered list of (lon, lat) vertices (at least 2). The whole polyline is encoded as a single VE edge so that GDAL reassembles it into a LineString geometry. For very long coastlines you can split into multiple calls, each producing a separate feature.""" if len(coords) < 2: raise ValueError("line feature needs ≥2 coordinates") edge_name = self.add_edge(coords) rcid = self._next_fid self._next_fid += 1 fields = [ ("0001", field_0001(rcid)), ("FRID", field_FRID(rcid=rcid, prim=PRIM_LINE, objl=objl)), ("FOID", field_FOID(agen=self.agen, fidn=rcid)), ] if attrs: fields.append(("ATTF", field_ATTF(attrs))) fields.append(("FSPT", field_FSPT([ (edge_name[0], edge_name[1], 1, 255, 255), # ORNT=1=fwd, USAG=null, MASK=null ]))) self._records.append(("FE", fields)) self._n_geo += 1 return rcid # --- Feature: area defined by closed boundary polygon --- def add_area_feature(self, *, objl: int, ring: list[tuple[float, float]], attrs: list[tuple[int, str]] | None = None) -> int: """ring: list of (lon, lat) — must be closed (first == last) or will be.""" if ring[0] != ring[-1]: ring = ring + [ring[0]] edge_name = self.add_edge(ring) rcid = self._next_fid self._next_fid += 1 fields = [ ("0001", field_0001(rcid)), ("FRID", field_FRID(rcid=rcid, prim=PRIM_AREA, objl=objl)), ("FOID", field_FOID(agen=self.agen, fidn=rcid)), ] if attrs: fields.append(("ATTF", field_ATTF(attrs))) fields.append(("FSPT", field_FSPT([ (edge_name[0], edge_name[1], 1, 1, 255), # ORNT=1=fwd, USAG=1=exterior, MASK=null ]))) self._records.append(("FE", fields)) self._n_geo += 1 return rcid # --- Build the full file --- def write(self, output_path: str | Path) -> None: if not DDR_TEMPLATE.exists(): raise FileNotFoundError( f"DDR template not found: {DDR_TEMPLATE}. " "Extract noaa_ddr_template.bin from a real ENC first." ) with open(DDR_TEMPLATE, "rb") as f: ddr = f.read() # Build DSID record (record 2) dsid_record = build_record([ ("0001", field_0001(2)), ("DSID", field_DSID( rcnm=10, rcid=1, expp=1, intu=self.intu, dsnm=self.dsnm, edtn=str(self.edition), updn="0", uadt=self.issue_date, isdt=self.issue_date, sted=3.1, prsp=1, psdn="", pred="2.0", prof=1, agen=self.agen, comt=self.comt, )), ("DSSI", field_DSSI( dstr=2, aall=1, nall=1, nomr=self._n_meta, nocr=self._n_col, nogr=self._n_gr, nolr=self._n_lr, noin=self._n_pt, nocn=self._n_cn, noed=self._n_ed, nofa=self._n_fa, )), ]) # Build DSPM record (record 3) dspm_record = build_record([ ("0001", field_0001(3)), ("DSPM", field_DSPM( rcnm=20, rcid=1, hdat=2, vdat=17, sdat=23, cscl=self.scale, duni=1, huni=1, puni=1, coun=1, comf=self.comf, somf=self.somf, comt="", )), ]) # Build all data records collected by the caller rcid_seq = 4 # records 1=DDR, 2=DSID, 3=DSPM, then sequential body_records = b"" for kind, fields in self._records: # Replace the 0001 RCID with global record sequence (per-record-type # numbering would also be valid; simplest is global sequence) new_fields = [ ("0001", field_0001(rcid_seq)) if t == "0001" else (t, d) for t, d in fields ] body_records += build_record(new_fields) rcid_seq += 1 with open(output_path, "wb") as f: f.write(ddr + dsid_record + dspm_record + body_records)