676 lines
25 KiB
Python
676 lines
25 KiB
Python
#!/usr/bin/env python3
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import datetime as dt
|
||
import json
|
||
import math
|
||
import struct
|
||
from pathlib import Path
|
||
from typing import Iterable
|
||
|
||
import pymysql
|
||
|
||
|
||
DB_NAME = "navsea_japan_coast_grid"
|
||
DB_USER = "root"
|
||
DB_PASSWORD = "2chi9ks2"
|
||
DB_HOST = "localhost"
|
||
DB_SOCKET = "/tmp/mysql.sock"
|
||
|
||
ZOOM = 12
|
||
GRID_W = 40
|
||
GRID_H = 40
|
||
CELL_COUNT = GRID_W * GRID_H
|
||
BLOCK_SIZE = 32
|
||
SPARSE_THRESHOLD = 96
|
||
MAGIC = b"CWZ12BIN"
|
||
VERSION = 1
|
||
|
||
LAYER_20 = "fish_port_20m"
|
||
LAYER_50 = "hazard_50m"
|
||
LAYER_200 = "coast_200m"
|
||
CELL_SIZES = {
|
||
LAYER_20: 20.0,
|
||
LAYER_50: 50.0,
|
||
LAYER_200: 200.0,
|
||
}
|
||
|
||
RADIUS = 6378137.0
|
||
WORLD_MIN = -math.pi * RADIUS
|
||
WORLD_MAX = math.pi * RADIUS
|
||
WORLD_SIZE = WORLD_MAX - WORLD_MIN
|
||
TILE_COUNT = 1 << ZOOM
|
||
TILE_SIZE_M = WORLD_SIZE / TILE_COUNT
|
||
Z12_CELL_SIZE_M = TILE_SIZE_M / GRID_W
|
||
GLOBAL_CELL_COUNT = TILE_COUNT * GRID_W
|
||
|
||
HEADER_STRUCT = struct.Struct("<8sHBBBBBIII")
|
||
BLOCK_ENTRY_STRUCT = struct.Struct("<IIHH")
|
||
TILE_ENTRY_STRUCT = struct.Struct("<IIHBB")
|
||
|
||
PRESET_BBOXES = {
|
||
# 九州本岛和周边主要岛屿的调试范围，足够覆盖福冈/佐贺/长崎/熊本/鹿儿岛一带。
|
||
"kyushu": (128.0, 30.0, 132.5, 34.9),
|
||
}
|
||
|
||
|
||
def log(message: str) -> None:
|
||
print(f"[{dt.datetime.now().isoformat(timespec='seconds')}] {message}")
|
||
|
||
|
||
def mysql_connect(database: str, *, cursorclass=pymysql.cursors.Cursor):
|
||
kwargs = {
|
||
"host": DB_HOST,
|
||
"user": DB_USER,
|
||
"password": DB_PASSWORD,
|
||
"database": database,
|
||
"charset": "utf8mb4",
|
||
"cursorclass": cursorclass,
|
||
}
|
||
if DB_SOCKET:
|
||
kwargs["unix_socket"] = DB_SOCKET
|
||
return pymysql.connect(**kwargs)
|
||
|
||
|
||
def mercator_x(lon: float) -> float:
|
||
return RADIUS * math.radians(lon)
|
||
|
||
|
||
def mercator_y(lat: float) -> float:
|
||
lat = max(min(lat, 85.0511287798066), -85.0511287798066)
|
||
return RADIUS * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0))
|
||
|
||
|
||
def lon_from_mercator(x: float) -> float:
|
||
return math.degrees(x / RADIUS)
|
||
|
||
|
||
def lat_from_mercator(y: float) -> float:
|
||
return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0)
|
||
|
||
|
||
def parse_bbox(value: str) -> tuple[float, float, float, float]:
|
||
parts = [float(part.strip()) for part in value.split(",") if part.strip()]
|
||
if len(parts) != 4:
|
||
raise argparse.ArgumentTypeError("bbox must be min_lon,min_lat,max_lon,max_lat")
|
||
min_lon, min_lat, max_lon, max_lat = parts
|
||
if min_lon >= max_lon or min_lat >= max_lat:
|
||
raise argparse.ArgumentTypeError("bbox min values must be smaller than max values")
|
||
return min_lon, min_lat, max_lon, max_lat
|
||
|
||
|
||
def parse_lonlat(value: str) -> tuple[float, float]:
|
||
parts = [float(part.strip()) for part in value.split(",") if part.strip()]
|
||
if len(parts) != 2:
|
||
raise argparse.ArgumentTypeError("query point must be lon,lat")
|
||
return parts[0], parts[1]
|
||
|
||
|
||
def mercator_bbox_from_lonlat(bbox: tuple[float, float, float, float]) -> tuple[float, float, float, float]:
|
||
min_lon, min_lat, max_lon, max_lat = bbox
|
||
return mercator_x(min_lon), mercator_y(min_lat), mercator_x(max_lon), mercator_y(max_lat)
|
||
|
||
|
||
def rowcol_range_for_bbox(
|
||
bbox_mercator: tuple[float, float, float, float],
|
||
cell_size_m: float,
|
||
) -> tuple[int, int, int, int]:
|
||
minx, miny, maxx, maxy = bbox_mercator
|
||
col_min = math.floor(minx / cell_size_m)
|
||
col_max = math.floor((maxx - 1e-9) / cell_size_m)
|
||
row_min = math.floor(miny / cell_size_m)
|
||
row_max = math.floor((maxy - 1e-9) / cell_size_m)
|
||
return row_min, row_max, col_min, col_max
|
||
|
||
|
||
def overlap_range(min_value: float, max_value: float, fine_size_m: float) -> tuple[int, int]:
|
||
start = int(math.floor(min_value / fine_size_m))
|
||
end = int(math.floor((max_value - 1e-9) / fine_size_m))
|
||
return start, end
|
||
|
||
|
||
def overlaps_finer_cells(
|
||
row_idx: int,
|
||
col_idx: int,
|
||
coarse_size_m: float,
|
||
fine_size_m: float,
|
||
fine_keys: set[tuple[int, int]],
|
||
) -> bool:
|
||
min_x = col_idx * coarse_size_m
|
||
min_y = row_idx * coarse_size_m
|
||
max_x = min_x + coarse_size_m
|
||
max_y = min_y + coarse_size_m
|
||
row_start, row_end = overlap_range(min_y, max_y, fine_size_m)
|
||
col_start, col_end = overlap_range(min_x, max_x, fine_size_m)
|
||
for fine_row in range(row_start, row_end + 1):
|
||
for fine_col in range(col_start, col_end + 1):
|
||
if (fine_row, fine_col) in fine_keys:
|
||
return True
|
||
return False
|
||
|
||
|
||
def pack_tile_id(tile_x: int, tile_y: int) -> int:
|
||
if not (0 <= tile_x <= 0xFFF and 0 <= tile_y <= 0xFFF):
|
||
raise ValueError(f"tile out of z12 range: {tile_x}, {tile_y}")
|
||
return (tile_x << 12) | tile_y
|
||
|
||
|
||
def unpack_tile_id(tile_id: int) -> tuple[int, int]:
|
||
return tile_id >> 12, tile_id & 0xFFF
|
||
|
||
|
||
def pack_cell_id(row: int, col: int) -> int:
|
||
if not (0 <= row < GRID_H and 0 <= col < GRID_W):
|
||
raise ValueError(f"cell out of tile range: {row}, {col}")
|
||
return row * GRID_W + col
|
||
|
||
|
||
def pack_block_id(tile_x: int, tile_y: int) -> int:
|
||
block_x = tile_x >> 5
|
||
block_y = tile_y >> 5
|
||
return (block_x << 7) | block_y
|
||
|
||
|
||
def global_cell_x_range(minx: float, maxx: float) -> tuple[int, int]:
|
||
start = math.floor((minx - WORLD_MIN) / Z12_CELL_SIZE_M)
|
||
end = math.floor(((maxx - 1e-9) - WORLD_MIN) / Z12_CELL_SIZE_M)
|
||
return max(0, start), min(GLOBAL_CELL_COUNT - 1, end)
|
||
|
||
|
||
def global_cell_y_range(miny: float, maxy: float) -> tuple[int, int]:
|
||
# WebMercator tile rows start at the north edge, while project grid rows use increasing mercator Y.
|
||
start = math.floor((WORLD_MAX - maxy) / Z12_CELL_SIZE_M)
|
||
end = math.floor(((WORLD_MAX - miny) - 1e-9) / Z12_CELL_SIZE_M)
|
||
return max(0, start), min(GLOBAL_CELL_COUNT - 1, end)
|
||
|
||
|
||
def add_mercator_rect_to_blocked(
|
||
blocked: dict[int, set[int]],
|
||
minx: float,
|
||
miny: float,
|
||
maxx: float,
|
||
maxy: float,
|
||
) -> int:
|
||
gx_start, gx_end = global_cell_x_range(minx, maxx)
|
||
gy_start, gy_end = global_cell_y_range(miny, maxy)
|
||
added = 0
|
||
for gy in range(gy_start, gy_end + 1):
|
||
tile_y = gy // GRID_H
|
||
row = gy % GRID_H
|
||
for gx in range(gx_start, gx_end + 1):
|
||
tile_x = gx // GRID_W
|
||
col = gx % GRID_W
|
||
tile_id = pack_tile_id(tile_x, tile_y)
|
||
cell_id = pack_cell_id(row, col)
|
||
cells = blocked.setdefault(tile_id, set())
|
||
before = len(cells)
|
||
cells.add(cell_id)
|
||
if len(cells) != before:
|
||
added += 1
|
||
return added
|
||
|
||
|
||
def iter_layer_rows(
|
||
conn,
|
||
*,
|
||
layer_name: str,
|
||
bbox_mercator: tuple[float, float, float, float] | None,
|
||
):
|
||
where = ["layer_name=%s"]
|
||
params: list[object] = [layer_name]
|
||
if bbox_mercator is not None:
|
||
row_min, row_max, col_min, col_max = rowcol_range_for_bbox(bbox_mercator, CELL_SIZES[layer_name])
|
||
where.append("row_idx BETWEEN %s AND %s")
|
||
where.append("col_idx BETWEEN %s AND %s")
|
||
params.extend([row_min, row_max, col_min, col_max])
|
||
sql = f"""
|
||
SELECT cell_id, row_idx, col_idx
|
||
FROM navsea_grid_cell
|
||
WHERE {' AND '.join(where)}
|
||
ORDER BY row_idx, col_idx, source_name
|
||
"""
|
||
with conn.cursor(pymysql.cursors.SSCursor) as cur:
|
||
cur.execute(sql, tuple(params))
|
||
for cell_id, row_idx, col_idx in cur:
|
||
yield str(cell_id), int(row_idx), int(col_idx)
|
||
|
||
|
||
def source_cell_mercator_bbox(row_idx: int, col_idx: int, cell_size_m: float) -> tuple[float, float, float, float]:
|
||
minx = col_idx * cell_size_m
|
||
miny = row_idx * cell_size_m
|
||
return minx, miny, minx + cell_size_m, miny + cell_size_m
|
||
|
||
|
||
def build_blocked_map(
|
||
conn,
|
||
*,
|
||
bbox_mercator: tuple[float, float, float, float] | None,
|
||
progress_every: int,
|
||
) -> tuple[dict[int, set[int]], dict]:
|
||
blocked: dict[int, set[int]] = {}
|
||
keys_20: set[tuple[int, int]] = set()
|
||
keys_50: set[tuple[int, int]] = set()
|
||
seen_20: set[str] = set()
|
||
seen_50: set[str] = set()
|
||
seen_200: set[str] = set()
|
||
|
||
summary = {
|
||
"sourceRows": {LAYER_20: 0, LAYER_50: 0, LAYER_200: 0},
|
||
"sourceDistinctCells": {LAYER_20: 0, LAYER_50: 0, LAYER_200: 0},
|
||
"keptByDensity": {"20": 0, "50": 0, "200": 0},
|
||
"skippedByDensity": {"50_over_20": 0, "200_over_20": 0, "200_over_50": 0},
|
||
"blockedCellAdds": {"20": 0, "50": 0, "200": 0},
|
||
}
|
||
|
||
log("读取 fish_port_20m；density 规则：20m 永远保留")
|
||
for cell_id, row_idx, col_idx in iter_layer_rows(conn, layer_name=LAYER_20, bbox_mercator=bbox_mercator):
|
||
summary["sourceRows"][LAYER_20] += 1
|
||
if cell_id in seen_20:
|
||
continue
|
||
seen_20.add(cell_id)
|
||
keys_20.add((row_idx, col_idx))
|
||
summary["sourceDistinctCells"][LAYER_20] += 1
|
||
summary["keptByDensity"]["20"] += 1
|
||
summary["blockedCellAdds"]["20"] += add_mercator_rect_to_blocked(
|
||
blocked,
|
||
*source_cell_mercator_bbox(row_idx, col_idx, CELL_SIZES[LAYER_20]),
|
||
)
|
||
if progress_every and summary["sourceRows"][LAYER_20] % progress_every == 0:
|
||
log(f"fish_port_20m rows={summary['sourceRows'][LAYER_20]} kept={summary['keptByDensity']['20']}")
|
||
|
||
log("读取 hazard_50m；density 规则：与 20m 重叠则跳过")
|
||
for cell_id, row_idx, col_idx in iter_layer_rows(conn, layer_name=LAYER_50, bbox_mercator=bbox_mercator):
|
||
summary["sourceRows"][LAYER_50] += 1
|
||
if cell_id in seen_50:
|
||
continue
|
||
seen_50.add(cell_id)
|
||
summary["sourceDistinctCells"][LAYER_50] += 1
|
||
if overlaps_finer_cells(row_idx, col_idx, CELL_SIZES[LAYER_50], CELL_SIZES[LAYER_20], keys_20):
|
||
summary["skippedByDensity"]["50_over_20"] += 1
|
||
continue
|
||
keys_50.add((row_idx, col_idx))
|
||
summary["keptByDensity"]["50"] += 1
|
||
summary["blockedCellAdds"]["50"] += add_mercator_rect_to_blocked(
|
||
blocked,
|
||
*source_cell_mercator_bbox(row_idx, col_idx, CELL_SIZES[LAYER_50]),
|
||
)
|
||
if progress_every and summary["sourceRows"][LAYER_50] % progress_every == 0:
|
||
log(f"hazard_50m rows={summary['sourceRows'][LAYER_50]} kept={summary['keptByDensity']['50']}")
|
||
|
||
log("读取 coast_200m；density 规则：与 20m/保留后的 50m 重叠则跳过")
|
||
for cell_id, row_idx, col_idx in iter_layer_rows(conn, layer_name=LAYER_200, bbox_mercator=bbox_mercator):
|
||
summary["sourceRows"][LAYER_200] += 1
|
||
if cell_id in seen_200:
|
||
continue
|
||
seen_200.add(cell_id)
|
||
summary["sourceDistinctCells"][LAYER_200] += 1
|
||
if overlaps_finer_cells(row_idx, col_idx, CELL_SIZES[LAYER_200], CELL_SIZES[LAYER_20], keys_20):
|
||
summary["skippedByDensity"]["200_over_20"] += 1
|
||
continue
|
||
if overlaps_finer_cells(row_idx, col_idx, CELL_SIZES[LAYER_200], CELL_SIZES[LAYER_50], keys_50):
|
||
summary["skippedByDensity"]["200_over_50"] += 1
|
||
continue
|
||
summary["keptByDensity"]["200"] += 1
|
||
summary["blockedCellAdds"]["200"] += add_mercator_rect_to_blocked(
|
||
blocked,
|
||
*source_cell_mercator_bbox(row_idx, col_idx, CELL_SIZES[LAYER_200]),
|
||
)
|
||
if progress_every and summary["sourceRows"][LAYER_200] % progress_every == 0:
|
||
log(f"coast_200m rows={summary['sourceRows'][LAYER_200]} kept={summary['keptByDensity']['200']}")
|
||
|
||
return blocked, summary
|
||
|
||
|
||
def encode_tile_cells(cells: set[int]) -> tuple[int, int, bytes]:
|
||
sorted_cells = sorted(cells)
|
||
if not sorted_cells:
|
||
raise ValueError("empty tile should not be encoded")
|
||
if sorted_cells[0] < 0 or sorted_cells[-1] >= CELL_COUNT:
|
||
raise ValueError(f"cellId out of range: {sorted_cells[0]}..{sorted_cells[-1]}")
|
||
if len(sorted_cells) <= SPARSE_THRESHOLD:
|
||
payload = struct.pack(f"<{len(sorted_cells)}H", *sorted_cells)
|
||
return 1, len(sorted_cells), payload
|
||
words = [0] * 50
|
||
for cell_id in sorted_cells:
|
||
word = cell_id >> 5
|
||
bit = cell_id & 31
|
||
words[word] |= 1 << bit
|
||
payload = struct.pack("<50I", *words)
|
||
return 2, len(sorted_cells), payload
|
||
|
||
|
||
def write_bin(path: Path, blocked: dict[int, set[int]]) -> dict:
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
grouped: dict[int, list[tuple[int, int, int, bytes]]] = {}
|
||
sparse_count = 0
|
||
bitset_count = 0
|
||
total_blocked_cells = 0
|
||
|
||
for tile_id in sorted(blocked):
|
||
cells = blocked[tile_id]
|
||
if not cells:
|
||
continue
|
||
tile_x, tile_y = unpack_tile_id(tile_id)
|
||
encoding, count, payload = encode_tile_cells(cells)
|
||
if encoding == 1:
|
||
sparse_count += 1
|
||
else:
|
||
bitset_count += 1
|
||
total_blocked_cells += count
|
||
block_id = pack_block_id(tile_x, tile_y)
|
||
grouped.setdefault(block_id, []).append((tile_id, encoding, count, payload))
|
||
|
||
block_ids = sorted(grouped)
|
||
block_count = len(block_ids)
|
||
block_dir_offset = HEADER_STRUCT.size
|
||
tile_dir_offset = block_dir_offset + block_count * BLOCK_ENTRY_STRUCT.size
|
||
tile_entry_count = sum(len(grouped[block_id]) for block_id in block_ids)
|
||
data_offset = tile_dir_offset + tile_entry_count * TILE_ENTRY_STRUCT.size
|
||
|
||
block_entries: list[tuple[int, int, int, int]] = []
|
||
tile_entries_by_block: dict[int, list[tuple[int, int, int, int]]] = {}
|
||
data_parts: list[bytes] = []
|
||
current_tile_dir_offset = tile_dir_offset
|
||
current_data_offset = 0
|
||
|
||
for block_id in block_ids:
|
||
tiles = sorted(grouped[block_id], key=lambda item: item[0])
|
||
block_entries.append((block_id, current_tile_dir_offset, len(tiles), 0))
|
||
current_tile_dir_offset += len(tiles) * TILE_ENTRY_STRUCT.size
|
||
tile_entries: list[tuple[int, int, int, int]] = []
|
||
for tile_id, encoding, count, payload in tiles:
|
||
tile_entries.append((tile_id, current_data_offset, count, encoding))
|
||
data_parts.append(payload)
|
||
current_data_offset += len(payload)
|
||
tile_entries_by_block[block_id] = tile_entries
|
||
|
||
with path.open("wb") as fh:
|
||
fh.write(
|
||
HEADER_STRUCT.pack(
|
||
MAGIC,
|
||
VERSION,
|
||
ZOOM,
|
||
GRID_W,
|
||
GRID_H,
|
||
BLOCK_SIZE,
|
||
0,
|
||
block_count,
|
||
block_dir_offset,
|
||
data_offset,
|
||
)
|
||
)
|
||
for entry in block_entries:
|
||
fh.write(BLOCK_ENTRY_STRUCT.pack(*entry))
|
||
for block_id in block_ids:
|
||
for tile_id, offset, count, encoding in tile_entries_by_block[block_id]:
|
||
fh.write(TILE_ENTRY_STRUCT.pack(tile_id, offset, count, encoding, 0))
|
||
for payload in data_parts:
|
||
fh.write(payload)
|
||
|
||
file_size = path.stat().st_size
|
||
return {
|
||
"file": str(path),
|
||
"fileSizeBytes": file_size,
|
||
"headerSizeBytes": HEADER_STRUCT.size,
|
||
"blockCount": block_count,
|
||
"totalBlockedTiles": sparse_count + bitset_count,
|
||
"totalBlockedCells": total_blocked_cells,
|
||
"sparseTileCount": sparse_count,
|
||
"bitsetTileCount": bitset_count,
|
||
"blockDirOffset": block_dir_offset,
|
||
"dataOffset": data_offset,
|
||
"dataBytes": current_data_offset,
|
||
}
|
||
|
||
|
||
class CoastWallReader:
|
||
def __init__(self, path: Path):
|
||
self.path = path
|
||
self.data = path.read_bytes()
|
||
header = HEADER_STRUCT.unpack_from(self.data, 0)
|
||
magic, version, zoom, grid_w, grid_h, block_size, _flags, block_count, block_dir_offset, data_offset = header
|
||
if magic != MAGIC:
|
||
raise ValueError(f"bad magic: {magic!r}")
|
||
if version != VERSION or zoom != ZOOM or grid_w != GRID_W or grid_h != GRID_H or block_size != BLOCK_SIZE:
|
||
raise ValueError("unsupported coast wall header")
|
||
self.block_count = block_count
|
||
self.block_dir_offset = block_dir_offset
|
||
self.data_offset = data_offset
|
||
self.blocks: dict[int, list[tuple[int, int, int, int]]] = {}
|
||
for idx in range(block_count):
|
||
offset = block_dir_offset + idx * BLOCK_ENTRY_STRUCT.size
|
||
block_id, tile_dir_offset, tile_count, _reserved = BLOCK_ENTRY_STRUCT.unpack_from(self.data, offset)
|
||
entries: list[tuple[int, int, int, int]] = []
|
||
for tile_idx in range(tile_count):
|
||
tile_offset = tile_dir_offset + tile_idx * TILE_ENTRY_STRUCT.size
|
||
tile_id, payload_offset, count, encoding, _tile_reserved = TILE_ENTRY_STRUCT.unpack_from(self.data, tile_offset)
|
||
entries.append((tile_id, payload_offset, count, encoding))
|
||
self.blocks[block_id] = entries
|
||
|
||
def tile_cells(self, tile_id: int) -> set[int]:
|
||
tile_x, tile_y = unpack_tile_id(tile_id)
|
||
block_id = pack_block_id(tile_x, tile_y)
|
||
for entry_tile_id, payload_offset, count, encoding in self.blocks.get(block_id, []):
|
||
if entry_tile_id != tile_id:
|
||
continue
|
||
start = self.data_offset + payload_offset
|
||
if encoding == 1:
|
||
values = struct.unpack_from(f"<{count}H", self.data, start)
|
||
return set(values)
|
||
if encoding == 2:
|
||
words = struct.unpack_from("<50I", self.data, start)
|
||
cells: set[int] = set()
|
||
for word_idx, word in enumerate(words):
|
||
if not word:
|
||
continue
|
||
for bit in range(32):
|
||
if word & (1 << bit):
|
||
cell_id = word_idx * 32 + bit
|
||
if cell_id < CELL_COUNT:
|
||
cells.add(cell_id)
|
||
return cells
|
||
raise ValueError(f"unknown encoding: {encoding}")
|
||
return set()
|
||
|
||
def is_blocked_lonlat(self, lon: float, lat: float) -> bool:
|
||
tile_x, tile_y, _row, _col, cell_id = z12_cell_for_lonlat(lon, lat)
|
||
tile_id = pack_tile_id(tile_x, tile_y)
|
||
return cell_id in self.tile_cells(tile_id)
|
||
|
||
|
||
def validate_bin(path: Path, blocked: dict[int, set[int]], sample_tiles: int = 200) -> None:
|
||
reader = CoastWallReader(path)
|
||
for idx, tile_id in enumerate(sorted(blocked)):
|
||
if idx >= sample_tiles:
|
||
break
|
||
expected = blocked[tile_id]
|
||
actual = reader.tile_cells(tile_id)
|
||
if actual != expected:
|
||
raise RuntimeError(f"reader roundtrip mismatch tile={tile_id}: expected={len(expected)} actual={len(actual)}")
|
||
|
||
|
||
def rect_geometry(min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> dict:
|
||
return {
|
||
"type": "Polygon",
|
||
"coordinates": [
|
||
[
|
||
[min_lon, min_lat],
|
||
[max_lon, min_lat],
|
||
[max_lon, max_lat],
|
||
[min_lon, max_lat],
|
||
[min_lon, min_lat],
|
||
]
|
||
],
|
||
}
|
||
|
||
|
||
def z12_cell_bounds(tile_x: int, tile_y: int, row: int, col: int) -> tuple[float, float, float, float]:
|
||
minx = WORLD_MIN + (tile_x * GRID_W + col) * Z12_CELL_SIZE_M
|
||
maxx = minx + Z12_CELL_SIZE_M
|
||
maxy = WORLD_MAX - (tile_y * GRID_H + row) * Z12_CELL_SIZE_M
|
||
miny = maxy - Z12_CELL_SIZE_M
|
||
return lon_from_mercator(minx), lat_from_mercator(miny), lon_from_mercator(maxx), lat_from_mercator(maxy)
|
||
|
||
|
||
def z12_cell_for_lonlat(lon: float, lat: float) -> tuple[int, int, int, int, int]:
|
||
x = mercator_x(lon)
|
||
y = mercator_y(lat)
|
||
gx = math.floor((x - WORLD_MIN) / Z12_CELL_SIZE_M)
|
||
gy = math.floor((WORLD_MAX - y) / Z12_CELL_SIZE_M)
|
||
if not (0 <= gx < GLOBAL_CELL_COUNT and 0 <= gy < GLOBAL_CELL_COUNT):
|
||
raise ValueError(f"lon/lat outside z12 grid: {lon}, {lat}")
|
||
tile_x = gx // GRID_W
|
||
tile_y = gy // GRID_H
|
||
col = gx % GRID_W
|
||
row = gy % GRID_H
|
||
cell_id = pack_cell_id(row, col)
|
||
return tile_x, tile_y, row, col, cell_id
|
||
|
||
|
||
def write_debug_geojson(path: Path, blocked: dict[int, set[int]], limit: int) -> int:
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
count = 0
|
||
with path.open("w", encoding="utf-8") as fh:
|
||
fh.write('{"type":"FeatureCollection","features":[\n')
|
||
first = True
|
||
for tile_id in sorted(blocked):
|
||
tile_x, tile_y = unpack_tile_id(tile_id)
|
||
for cell_id in sorted(blocked[tile_id]):
|
||
row = cell_id // GRID_W
|
||
col = cell_id % GRID_W
|
||
feature = {
|
||
"type": "Feature",
|
||
"properties": {
|
||
"tile_id": tile_id,
|
||
"tile_x": tile_x,
|
||
"tile_y": tile_y,
|
||
"cell_id": cell_id,
|
||
"row": row,
|
||
"col": col,
|
||
},
|
||
"geometry": rect_geometry(*z12_cell_bounds(tile_x, tile_y, row, col)),
|
||
}
|
||
if not first:
|
||
fh.write(",\n")
|
||
fh.write(json.dumps(feature, ensure_ascii=False))
|
||
first = False
|
||
count += 1
|
||
if limit and count >= limit:
|
||
fh.write("\n]}\n")
|
||
return count
|
||
fh.write("\n]}\n")
|
||
return count
|
||
|
||
|
||
def run_self_test() -> None:
|
||
assert unpack_tile_id(pack_tile_id(3456, 789)) == (3456, 789)
|
||
assert pack_cell_id(39, 39) == 1599
|
||
assert pack_block_id(32, 64) == ((1 << 7) | 2)
|
||
sparse_cells = {0, 3, 1599}
|
||
encoding, count, payload = encode_tile_cells(sparse_cells)
|
||
assert encoding == 1 and count == 3
|
||
assert set(struct.unpack("<3H", payload)) == sparse_cells
|
||
dense_cells = set(range(120))
|
||
encoding, count, payload = encode_tile_cells(dense_cells)
|
||
assert encoding == 2 and count == 120 and len(payload) == 200
|
||
log("self-test passed")
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description="Build coast_wall_z12.bin from MySQL density-overview grid layers")
|
||
parser.add_argument("--db-name", default=DB_NAME)
|
||
parser.add_argument("--out", default="out/navmask/coast_wall_z12.bin")
|
||
parser.add_argument("--meta-json", default=None, help="默认写到 <out>.meta.json")
|
||
parser.add_argument("--debug-geojson", default=None)
|
||
parser.add_argument("--debug-limit", type=int, default=5000)
|
||
parser.add_argument("--preset", choices=["all", *sorted(PRESET_BBOXES)], default="kyushu")
|
||
parser.add_argument("--bbox", type=parse_bbox, default=None, help="min_lon,min_lat,max_lon,max_lat；优先于 --preset")
|
||
parser.add_argument("--query", type=parse_lonlat, default=None, help="生成后用 reader 查询一个 lon,lat 点")
|
||
parser.add_argument("--progress-every", type=int, default=200000)
|
||
parser.add_argument("--self-test", action="store_true")
|
||
args = parser.parse_args()
|
||
|
||
if args.self_test:
|
||
run_self_test()
|
||
|
||
bbox_lonlat = args.bbox
|
||
if bbox_lonlat is None and args.preset != "all":
|
||
bbox_lonlat = PRESET_BBOXES[args.preset]
|
||
bbox_mercator = mercator_bbox_from_lonlat(bbox_lonlat) if bbox_lonlat is not None else None
|
||
|
||
out_path = Path(args.out)
|
||
meta_path = Path(args.meta_json) if args.meta_json else out_path.with_suffix(out_path.suffix + ".meta.json")
|
||
|
||
log(f"开始生成 coast wall z12 bin preset={args.preset} bbox={bbox_lonlat}")
|
||
conn = mysql_connect(args.db_name)
|
||
try:
|
||
blocked, density_summary = build_blocked_map(
|
||
conn,
|
||
bbox_mercator=bbox_mercator,
|
||
progress_every=args.progress_every,
|
||
)
|
||
finally:
|
||
conn.close()
|
||
|
||
if not blocked:
|
||
raise SystemExit("no blocked cells produced")
|
||
|
||
write_summary = write_bin(out_path, blocked)
|
||
validate_bin(out_path, blocked)
|
||
debug_count = 0
|
||
if args.debug_geojson:
|
||
debug_count = write_debug_geojson(Path(args.debug_geojson), blocked, args.debug_limit)
|
||
query_result = None
|
||
if args.query is not None:
|
||
lon, lat = args.query
|
||
reader = CoastWallReader(out_path)
|
||
tile_x, tile_y, row, col, cell_id = z12_cell_for_lonlat(lon, lat)
|
||
query_result = {
|
||
"lon": lon,
|
||
"lat": lat,
|
||
"tile_x": tile_x,
|
||
"tile_y": tile_y,
|
||
"row": row,
|
||
"col": col,
|
||
"cell_id": cell_id,
|
||
"blocked": reader.is_blocked_lonlat(lon, lat),
|
||
}
|
||
|
||
meta = {
|
||
"generated_at": dt.datetime.now().isoformat(timespec="seconds"),
|
||
"db_name": args.db_name,
|
||
"preset": args.preset,
|
||
"bbox_lonlat": bbox_lonlat,
|
||
"zoom": ZOOM,
|
||
"grid_w": GRID_W,
|
||
"grid_h": GRID_H,
|
||
"z12_cell_size_m": Z12_CELL_SIZE_M,
|
||
"density_rule": [
|
||
"keep fish_port_20m first",
|
||
"keep hazard_50m only when it does not overlap kept fish_port_20m",
|
||
"keep coast_200m only when it does not overlap kept fish_port_20m or kept hazard_50m",
|
||
],
|
||
"density_summary": density_summary,
|
||
"binary_summary": write_summary,
|
||
"debug_geojson": args.debug_geojson,
|
||
"debug_feature_count": debug_count,
|
||
"query_result": query_result,
|
||
}
|
||
meta_path.parent.mkdir(parents=True, exist_ok=True)
|
||
meta_path.write_text(json.dumps(meta, ensure_ascii=False, indent=2), encoding="utf-8")
|
||
|
||
log(f"完成 bin={out_path} meta={meta_path}")
|
||
log(
|
||
"统计 "
|
||
f"tiles={write_summary['totalBlockedTiles']} cells={write_summary['totalBlockedCells']} "
|
||
f"sparse={write_summary['sparseTileCount']} bitset={write_summary['bitsetTileCount']} "
|
||
f"blocks={write_summary['blockCount']} size={write_summary['fileSizeBytes']}"
|
||
)
|
||
if query_result is not None:
|
||
log(f"query lon={query_result['lon']} lat={query_result['lat']} blocked={query_result['blocked']}")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|