871 lines
34 KiB
Python
871 lines
34 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,
|
||
}
|
||
TABLES = {
|
||
LAYER_20: "navsea_fish_port_grid_cell",
|
||
LAYER_50: "navsea_hazard_grid_cell",
|
||
LAYER_200: "navsea_coast_grid_cell",
|
||
}
|
||
|
||
|
||
def table_name_for_layer(layer_name: str) -> str:
|
||
if layer_name in TABLES:
|
||
return TABLES[layer_name]
|
||
if layer_name.startswith("coast_") and layer_name.endswith("m"):
|
||
return "navsea_coast_grid_cell"
|
||
raise KeyError(layer_name)
|
||
|
||
|
||
def coast_layer_name_for_size(cell_size_m: float) -> str:
|
||
if cell_size_m <= 0:
|
||
raise ValueError("coast cell size must be positive")
|
||
if float(cell_size_m).is_integer():
|
||
return f"coast_{int(cell_size_m)}m"
|
||
return f"coast_{cell_size_m:g}m"
|
||
|
||
|
||
def density_key_for_size(cell_size_m: float) -> str:
|
||
if float(cell_size_m).is_integer():
|
||
return str(int(cell_size_m))
|
||
return f"{cell_size_m:g}"
|
||
|
||
|
||
def unit_tag_for_size(cell_size_m: float) -> str:
|
||
if float(cell_size_m).is_integer():
|
||
return str(int(cell_size_m))
|
||
return str(cell_size_m).replace(".", "_")
|
||
|
||
|
||
def wall_unit_tag(zoom: int, grid_w: int, grid_h: int) -> str:
|
||
if grid_w == grid_h:
|
||
return f"z{zoom}_g{grid_w}"
|
||
return f"z{zoom}_g{grid_w}x{grid_h}"
|
||
|
||
|
||
def build_out_path(path: Path, coast_unit: str, wall_unit: str, *, include_fish_port: bool) -> Path:
|
||
source_unit = f"{coast_unit}_50_20" if include_fish_port else f"{coast_unit}_50"
|
||
stem = path.stem
|
||
if wall_unit != "z12_g40" and stem.endswith("_z12"):
|
||
stem = stem.removesuffix("_z12")
|
||
suffix = f"_{source_unit}" if wall_unit == "z12_g40" else f"_{wall_unit}_{source_unit}"
|
||
return path.with_name(f"{stem}{suffix}{path.suffix}")
|
||
|
||
|
||
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
|
||
Z12_CELL_SIZE_X_M = TILE_SIZE_M / GRID_W
|
||
Z12_CELL_SIZE_Y_M = TILE_SIZE_M / GRID_H
|
||
GLOBAL_CELL_COUNT_X = TILE_COUNT * GRID_W
|
||
GLOBAL_CELL_COUNT_Y = TILE_COUNT * GRID_H
|
||
BLOCK_AXIS_BITS = math.ceil(math.log2(math.ceil(TILE_COUNT / BLOCK_SIZE)))
|
||
|
||
|
||
def configure_wall_grid(*, zoom: int, grid_w: int, grid_h: int, sparse_threshold: int) -> None:
|
||
global ZOOM
|
||
global GRID_W
|
||
global GRID_H
|
||
global CELL_COUNT
|
||
global SPARSE_THRESHOLD
|
||
global TILE_COUNT
|
||
global TILE_SIZE_M
|
||
global Z12_CELL_SIZE_M
|
||
global GLOBAL_CELL_COUNT_X
|
||
global GLOBAL_CELL_COUNT_Y
|
||
global Z12_CELL_SIZE_X_M
|
||
global Z12_CELL_SIZE_Y_M
|
||
global BLOCK_AXIS_BITS
|
||
|
||
if not (0 <= zoom <= 16):
|
||
raise SystemExit("--zoom must be between 0 and 16")
|
||
if not (1 <= grid_w <= 255 and 1 <= grid_h <= 255):
|
||
raise SystemExit("--grid-w/--grid-h must be between 1 and 255")
|
||
if grid_w * grid_h > 65535:
|
||
raise SystemExit("--grid-w * --grid-h must be <= 65535 because cell ids are uint16")
|
||
if sparse_threshold < 0 or sparse_threshold > grid_w * grid_h:
|
||
raise SystemExit("--sparse-threshold must be between 0 and grid_w * grid_h")
|
||
|
||
ZOOM = zoom
|
||
GRID_W = grid_w
|
||
GRID_H = grid_h
|
||
CELL_COUNT = GRID_W * GRID_H
|
||
SPARSE_THRESHOLD = sparse_threshold
|
||
TILE_COUNT = 1 << ZOOM
|
||
TILE_SIZE_M = WORLD_SIZE / TILE_COUNT
|
||
Z12_CELL_SIZE_M = TILE_SIZE_M / GRID_W
|
||
Z12_CELL_SIZE_X_M = TILE_SIZE_M / GRID_W
|
||
Z12_CELL_SIZE_Y_M = TILE_SIZE_M / GRID_H
|
||
GLOBAL_CELL_COUNT_X = TILE_COUNT * GRID_W
|
||
GLOBAL_CELL_COUNT_Y = TILE_COUNT * GRID_H
|
||
BLOCK_AXIS_BITS = math.ceil(math.log2(math.ceil(TILE_COUNT / BLOCK_SIZE)))
|
||
|
||
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 < TILE_COUNT and 0 <= tile_y < TILE_COUNT):
|
||
raise ValueError(f"tile out of z{ZOOM} range: {tile_x}, {tile_y}")
|
||
return (tile_x << ZOOM) | tile_y
|
||
|
||
|
||
def unpack_tile_id(tile_id: int) -> tuple[int, int]:
|
||
mask = (1 << ZOOM) - 1
|
||
return tile_id >> ZOOM, tile_id & mask
|
||
|
||
|
||
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 << BLOCK_AXIS_BITS) | block_y
|
||
|
||
|
||
def global_cell_x_range(minx: float, maxx: float) -> tuple[int, int]:
|
||
start = math.floor((minx - WORLD_MIN) / Z12_CELL_SIZE_X_M)
|
||
end = math.floor(((maxx - 1e-9) - WORLD_MIN) / Z12_CELL_SIZE_X_M)
|
||
return max(0, start), min(GLOBAL_CELL_COUNT_X - 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_Y_M)
|
||
end = math.floor(((WORLD_MAX - miny) - 1e-9) / Z12_CELL_SIZE_Y_M)
|
||
return max(0, start), min(GLOBAL_CELL_COUNT_Y - 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])
|
||
table_name = table_name_for_layer(layer_name)
|
||
sql = f"""
|
||
SELECT cell_id, row_idx, col_idx
|
||
FROM {table_name}
|
||
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,
|
||
*,
|
||
coast_layer_name: str,
|
||
coast_cell_size_m: float,
|
||
include_fish_port: bool,
|
||
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()
|
||
coast_key = density_key_for_size(coast_cell_size_m)
|
||
|
||
summary = {
|
||
"sourceRows": {LAYER_20: 0, LAYER_50: 0, coast_layer_name: 0},
|
||
"sourceDistinctCells": {LAYER_20: 0, LAYER_50: 0, coast_layer_name: 0},
|
||
"keptByDensity": {"20": 0, "50": 0, coast_key: 0},
|
||
"skippedByDensity": {f"50_over_20": 0, f"{coast_key}_over_20": 0, f"{coast_key}_over_50": 0},
|
||
"blockedCellAdds": {"20": 0, "50": 0, coast_key: 0},
|
||
"activeLayers": [LAYER_50, coast_layer_name],
|
||
"disabledLayers": [],
|
||
}
|
||
|
||
if include_fish_port:
|
||
summary["activeLayers"].insert(0, LAYER_20)
|
||
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']}")
|
||
else:
|
||
summary["disabledLayers"].append(
|
||
{
|
||
"layer": LAYER_20,
|
||
"reason": f"coast_cell_size_m={coast_cell_size_m:g} < fish_port_20m precision",
|
||
}
|
||
)
|
||
log(f"跳过 fish_port_20m；原因：coast_cell_size_m={coast_cell_size_m:g} 已细于 20m")
|
||
|
||
hazard_rule = "与 20m 重叠则跳过" if include_fish_port else "不再与 20m 做重叠裁剪"
|
||
log(f"读取 hazard_50m；density 规则：{hazard_rule}")
|
||
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 include_fish_port and 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']}")
|
||
|
||
coast_rule = "与 20m/保留后的 50m 重叠则跳过" if include_fish_port else "与保留后的 50m 重叠则跳过"
|
||
log(f"读取 {coast_layer_name}；density 规则：{coast_rule}")
|
||
for cell_id, row_idx, col_idx in iter_layer_rows(conn, layer_name=coast_layer_name, bbox_mercator=bbox_mercator):
|
||
summary["sourceRows"][coast_layer_name] += 1
|
||
if cell_id in seen_200:
|
||
continue
|
||
seen_200.add(cell_id)
|
||
summary["sourceDistinctCells"][coast_layer_name] += 1
|
||
if include_fish_port and overlaps_finer_cells(row_idx, col_idx, coast_cell_size_m, CELL_SIZES[LAYER_20], keys_20):
|
||
summary["skippedByDensity"][f"{coast_key}_over_20"] += 1
|
||
continue
|
||
if overlaps_finer_cells(row_idx, col_idx, coast_cell_size_m, CELL_SIZES[LAYER_50], keys_50):
|
||
summary["skippedByDensity"][f"{coast_key}_over_50"] += 1
|
||
continue
|
||
summary["keptByDensity"][coast_key] += 1
|
||
summary["blockedCellAdds"][coast_key] += add_mercator_rect_to_blocked(
|
||
blocked,
|
||
*source_cell_mercator_bbox(row_idx, col_idx, coast_cell_size_m),
|
||
)
|
||
if progress_every and summary["sourceRows"][coast_layer_name] % progress_every == 0:
|
||
log(f"{coast_layer_name} rows={summary['sourceRows'][coast_layer_name]} kept={summary['keptByDensity'][coast_key]}")
|
||
|
||
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
|
||
word_count = math.ceil(CELL_COUNT / 32)
|
||
words = [0] * word_count
|
||
for cell_id in sorted_cells:
|
||
word = cell_id >> 5
|
||
bit = cell_id & 31
|
||
words[word] |= 1 << bit
|
||
payload = struct.pack(f"<{word_count}I", *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: "
|
||
f"file version={version} zoom={zoom} grid={grid_w}x{grid_h} block={block_size}; "
|
||
f"reader expects version={VERSION} zoom={ZOOM} grid={GRID_W}x{GRID_H} block={BLOCK_SIZE}"
|
||
)
|
||
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:
|
||
word_count = math.ceil(CELL_COUNT / 32)
|
||
words = struct.unpack_from(f"<{word_count}I", 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_X_M
|
||
maxx = minx + Z12_CELL_SIZE_X_M
|
||
maxy = WORLD_MAX - (tile_y * GRID_H + row) * Z12_CELL_SIZE_Y_M
|
||
miny = maxy - Z12_CELL_SIZE_Y_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_X_M)
|
||
gy = math.floor((WORLD_MAX - y) / Z12_CELL_SIZE_Y_M)
|
||
if not (0 <= gx < GLOBAL_CELL_COUNT_X and 0 <= gy < GLOBAL_CELL_COUNT_Y):
|
||
raise ValueError(f"lon/lat outside z{ZOOM} 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:
|
||
sample_tile_x = min(TILE_COUNT - 1, max(0, TILE_COUNT * 3 // 4))
|
||
sample_tile_y = min(TILE_COUNT - 1, max(0, TILE_COUNT // 5))
|
||
assert unpack_tile_id(pack_tile_id(sample_tile_x, sample_tile_y)) == (sample_tile_x, sample_tile_y)
|
||
assert pack_cell_id(GRID_H - 1, GRID_W - 1) == CELL_COUNT - 1
|
||
assert pack_block_id(32, 64) == ((1 << BLOCK_AXIS_BITS) | 2)
|
||
sparse_cells = {0, min(3, CELL_COUNT - 1), CELL_COUNT - 1}
|
||
encoding, count, payload = encode_tile_cells(sparse_cells)
|
||
assert encoding == 1 and count == len(sparse_cells)
|
||
assert set(struct.unpack(f"<{len(sparse_cells)}H", payload)) == sparse_cells
|
||
if SPARSE_THRESHOLD < CELL_COUNT:
|
||
dense_count = min(CELL_COUNT, max(SPARSE_THRESHOLD + 1, 120))
|
||
dense_cells = set(range(dense_count))
|
||
encoding, count, payload = encode_tile_cells(dense_cells)
|
||
assert encoding == 2 and count == dense_count and len(payload) == math.ceil(CELL_COUNT / 32) * 4
|
||
log("self-test passed")
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description="Build coast wall 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("--zoom", type=int, default=12, help="阻断墙输出 tile zoom，默认 12")
|
||
parser.add_argument(
|
||
"--grid-size",
|
||
type=int,
|
||
default=40,
|
||
help="每个输出 tile 内部切成 N x N 个墙 cell，默认 40；值越大墙 cell 越小",
|
||
)
|
||
parser.add_argument("--grid-w", type=int, default=None, help="覆盖输出 tile 内部横向 cell 数")
|
||
parser.add_argument("--grid-h", type=int, default=None, help="覆盖输出 tile 内部纵向 cell 数")
|
||
parser.add_argument(
|
||
"--sparse-threshold",
|
||
type=int,
|
||
default=96,
|
||
help="单 tile blocked cell 数量小于等于该值时使用 sparse 编码，默认 96",
|
||
)
|
||
parser.add_argument(
|
||
"--coast-cell-size-m",
|
||
type=float,
|
||
default=200.0,
|
||
help="参与 bin 生成的海岸线格网大小，默认 200；可切 100、50、20、10 等已入库层",
|
||
)
|
||
parser.add_argument(
|
||
"--coast-layer-name",
|
||
default=None,
|
||
help="覆盖海岸线层名；不传时按 --coast-cell-size-m 自动生成，例如 coast_200m / coast_50m / coast_10m",
|
||
)
|
||
parser.add_argument("--self-test", action="store_true")
|
||
args = parser.parse_args()
|
||
|
||
grid_w = args.grid_w if args.grid_w is not None else args.grid_size
|
||
grid_h = args.grid_h if args.grid_h is not None else args.grid_size
|
||
configure_wall_grid(
|
||
zoom=args.zoom,
|
||
grid_w=grid_w,
|
||
grid_h=grid_h,
|
||
sparse_threshold=args.sparse_threshold,
|
||
)
|
||
|
||
if args.self_test:
|
||
run_self_test()
|
||
if args.coast_cell_size_m <= 0:
|
||
raise SystemExit("--coast-cell-size-m must be positive")
|
||
|
||
coast_layer_name = args.coast_layer_name or coast_layer_name_for_size(args.coast_cell_size_m)
|
||
coast_cell_size_m = float(args.coast_cell_size_m)
|
||
if coast_layer_name in CELL_SIZES and CELL_SIZES[coast_layer_name] != coast_cell_size_m:
|
||
raise SystemExit(
|
||
f"{coast_layer_name} already maps to {CELL_SIZES[coast_layer_name]}m, "
|
||
f"but --coast-cell-size-m={coast_cell_size_m}"
|
||
)
|
||
CELL_SIZES.setdefault(coast_layer_name, coast_cell_size_m)
|
||
coast_unit_tag = unit_tag_for_size(coast_cell_size_m)
|
||
wall_unit = wall_unit_tag(ZOOM, GRID_W, GRID_H)
|
||
include_fish_port = coast_cell_size_m >= CELL_SIZES[LAYER_20]
|
||
|
||
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)
|
||
out_path = build_out_path(out_path, coast_unit_tag, wall_unit, include_fish_port=include_fish_port)
|
||
meta_path = Path(args.meta_json) if args.meta_json else out_path.with_suffix(out_path.suffix + ".meta.json")
|
||
|
||
log(
|
||
f"开始生成 coast wall bin zoom={ZOOM} grid={GRID_W}x{GRID_H} "
|
||
f"wall_cell_m={Z12_CELL_SIZE_X_M:.3f}x{Z12_CELL_SIZE_Y_M:.3f} "
|
||
f"preset={args.preset} bbox={bbox_lonlat} "
|
||
f"coast_layer={coast_layer_name} include_fish_port={include_fish_port} out={out_path.name}"
|
||
)
|
||
conn = mysql_connect(args.db_name)
|
||
try:
|
||
blocked, density_summary = build_blocked_map(
|
||
conn,
|
||
coast_layer_name=coast_layer_name,
|
||
coast_cell_size_m=coast_cell_size_m,
|
||
include_fish_port=include_fish_port,
|
||
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,
|
||
"coast_layer_name": coast_layer_name,
|
||
"coast_cell_size_m": coast_cell_size_m,
|
||
"include_fish_port_20m": include_fish_port,
|
||
"zoom": ZOOM,
|
||
"grid_w": GRID_W,
|
||
"grid_h": GRID_H,
|
||
"wall_cell_size_m": Z12_CELL_SIZE_M,
|
||
"wall_cell_size_x_m": Z12_CELL_SIZE_X_M,
|
||
"wall_cell_size_y_m": Z12_CELL_SIZE_Y_M,
|
||
"z12_cell_size_m": Z12_CELL_SIZE_M,
|
||
"sparse_threshold": SPARSE_THRESHOLD,
|
||
"density_rule": (
|
||
[
|
||
"keep fish_port_20m first",
|
||
"keep hazard_50m only when it does not overlap kept fish_port_20m",
|
||
f"keep {coast_layer_name} only when it does not overlap kept fish_port_20m or kept hazard_50m",
|
||
]
|
||
if include_fish_port
|
||
else [
|
||
"skip fish_port_20m because coast_cell_size_m is finer than 20m",
|
||
"keep hazard_50m",
|
||
f"keep {coast_layer_name} only when it does not overlap 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()
|