Files
pbf/coastline/export_prc20_test_assets.py

219 lines
6.6 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
#!/usr/bin/env python3
from __future__ import annotations
import argparse
import datetime as dt
import json
import math
from pathlib import Path
import pymysql
DB_NAME = "navsea_japan_coast_grid"
DB_USER = "root"
DB_PASSWORD = "2chi9ks2"
DB_HOST = "localhost"
DB_SOCKET = "/tmp/mysql.sock"
DEFAULT_OUT_DIR = "src/pbf/coastline-mysql/prc50_test"
RADIUS = 6378137.0
TARGET_LAYER = "fish_port_20m"
DEFAULT_PRC = "50"
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 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 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 mercator_bbox_to_lonlat(minx: float, miny: float, maxx: float, maxy: float) -> tuple[float, float, float, float]:
return (
lon_from_mercator(minx),
lat_from_mercator(miny),
lon_from_mercator(maxx),
lat_from_mercator(maxy),
)
def normalize_bbox(layer_name: str, min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> tuple[float, float, float, float]:
if layer_name == "coast_200m" and (
abs(min_lon) > 180.0 or abs(max_lon) > 180.0 or abs(min_lat) > 90.0 or abs(max_lat) > 90.0
):
return mercator_bbox_to_lonlat(min_lon, min_lat, max_lon, max_lat)
return min_lon, min_lat, max_lon, max_lat
def export_layer(
conn,
*,
path: Path,
where_sql: str,
params: tuple,
) -> int:
path.parent.mkdir(parents=True, exist_ok=True)
sql = f"""
SELECT layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name,
source_name, min_lon, min_lat, max_lon, max_lat
FROM navsea_grid_cell
WHERE {where_sql}
ORDER BY row_idx, col_idx
"""
count = 0
with conn.cursor(pymysql.cursors.SSCursor) as cur, path.open("w", encoding="utf-8") as fh:
cur.execute(sql, params)
fh.write('{"type":"FeatureCollection","features":[\n')
first = True
for row in cur:
(
layer_name,
cell_id,
row_idx,
col_idx,
cell_size_m,
state_name,
class_name,
source_name,
min_lon,
min_lat,
max_lon,
max_lat,
) = row
min_lon, min_lat, max_lon, max_lat = normalize_bbox(
row[0],
float(min_lon),
float(min_lat),
float(max_lon),
float(max_lat),
)
feature = {
"type": "Feature",
"properties": {
"layer_name": layer_name,
"cell_id": cell_id,
"row": int(row_idx),
"col": int(col_idx),
"cell_size_m": float(cell_size_m),
"state_name": state_name,
"class_name": class_name,
"source_name": source_name,
},
"geometry": rect_geometry(min_lon, min_lat, max_lon, max_lat),
}
if not first:
fh.write(",\n")
fh.write(json.dumps(feature, ensure_ascii=False))
first = False
count += 1
fh.write("\n]}\n")
return count
def fetch_stats(conn, where_sql: str, params: tuple) -> tuple[int, list[float] | None]:
sql = f"""
SELECT COUNT(*), MIN(min_lon), MIN(min_lat), MAX(max_lon), MAX(max_lat)
FROM navsea_grid_cell
WHERE {where_sql}
"""
with conn.cursor() as cur:
cur.execute(sql, params)
row = cur.fetchone()
count = int(row[0] or 0)
if count == 0:
return 0, None
min_lon, min_lat, max_lon, max_lat = normalize_bbox(
params[0] if params else "",
float(row[1]),
float(row[2]),
float(row[3]),
float(row[4]),
)
return count, [min_lon, min_lat, max_lon, max_lat]
def main() -> None:
parser = argparse.ArgumentParser(description="Export a single PRC fish_port_20m test subset")
parser.add_argument("--db-name", default=DB_NAME)
parser.add_argument("--out-dir", default=DEFAULT_OUT_DIR)
parser.add_argument("--prc", default=DEFAULT_PRC, help="只导出这个 PRC 的 fish_port_20m,默认 50")
args = parser.parse_args()
project_root = Path(__file__).resolve().parent.parent
out_dir = Path(args.out_dir)
if not out_dir.is_absolute():
out_dir = project_root / out_dir
out_dir.mkdir(parents=True, exist_ok=True)
prc = str(args.prc).strip()
if prc.isdigit():
prc = f"{int(prc):02d}"
where_sql = "layer_name=%s AND state_name=%s AND source_name LIKE %s"
params = (TARGET_LAYER, "COASTLINE", f"C09-06+coastline PRC={prc}%")
out_path = out_dir / f"fish_port_20m_prc{prc}_grid.geojson"
conn = mysql_connect(args.db_name)
try:
count = export_layer(conn, path=out_path, where_sql=where_sql, params=params)
stat_count, bbox = fetch_stats(conn, where_sql, params)
if count != stat_count:
raise RuntimeError(f"export count mismatch: {count} != {stat_count}")
finally:
conn.close()
manifest = {
"database": args.db_name,
"generated_at": dt.datetime.now().isoformat(timespec="seconds"),
"prc": prc,
"layers": {
TARGET_LAYER: {
"count": count,
"bbox_lonlat": bbox,
"cell_size_m": 20.0,
"export_file": str(out_path.relative_to(project_root)),
"filter": f"state_name='COASTLINE' AND source_name LIKE 'C09-06+coastline PRC={prc}%'",
}
},
}
manifest_path = out_dir / "manifest.json"
manifest_path.write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8")
print(f"{TARGET_LAYER}: {count} -> {out_path}")
print(f"manifest: {manifest_path}")
if __name__ == "__main__":
main()