diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts b/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts index 914792aef..b49bdf12a 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts +++ b/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts @@ -236,10 +236,21 @@ export function drawRaster( map: ScreenMap, raster: { raster: PlanRaster; image: HTMLImageElement }, ): void { + // 세 귀퉁이로 놓음 — TM 격자는 경위도 바탕에 대해 수렴각만큼 돌아 있음(586-0) const [x0, y0, x1, y1] = raster.raster.box_m; const [left, top] = map.toScreen([x0, y1]); - const [right, bottom] = map.toScreen([x1, y0]); - context.drawImage(raster.image, left, top, right - left, bottom - top); + const [right, rightY] = map.toScreen([x1, y1]); + const [bottomX, bottom] = map.toScreen([x0, y0]); + const { naturalWidth: w, naturalHeight: h } = raster.image; + context.transform( + (right - left) / w, + (rightY - top) / w, + (bottomX - left) / h, + (bottom - top) / h, + left, + top, + ); + context.drawImage(raster.image, 0, 0); } /** 켬끔 취향(`key`) — 없거나 깨지면 처음 값. */ diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_UI_Map.ts b/D03_PlanAnalysis/D03_PlanAnalysis_UI_Map.ts index ac8828bb4..e539e1e2e 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_UI_Map.ts +++ b/D03_PlanAnalysis/D03_PlanAnalysis_UI_Map.ts @@ -10,6 +10,7 @@ * 겹은 돌리기 전 그림 좌표로 그리고 받음 — 누른 자리는 `unrotate` 로 되돌려 넘김. * 캔버스 두 장(582) — 아래 = 위성 · 바탕 겹(`base`) · 위 = 그은 선 · 재기 · 누름 받음. * 마우스 움직임 = 위만 다시 그림 · 화면이 바뀌면(옮김 · 확대 · 돌림 · 겹) 둘 다 · 한 화면 한 번(rAF). + * TM 겹(그은 선 · 재기 · 그림 겹)은 경위도로 바꿔 바탕과 같은 자리(586-0 · `localTmFrame`) — 좌표계를 모르면 옛 꼴. * ========================================================================== */ import { createMapTools, mapToolButton } from "@ui/ui_template_map_tools"; @@ -19,6 +20,8 @@ import { computeMapRect, computeRouteView, createMetricProjector, + createNormalizer, + lonLatToScreen, normalizedToScreen, routeLineColor, type ViewState, @@ -27,6 +30,7 @@ import { bindRouteEditNavigation } from "../B05_Profile/B05_Profile_UI_RouteEdit import { createMapRotation } from "../B05_Profile/B05_Profile_UI_RouteEdit_Rotate"; import { fitCanvas, screenToMetric } from "../B05_Profile/B05_Profile_UI_RouteEdit_View"; import type { Vertex } from "./logic/D03_PlanAnalysis_Plan"; +import { localTmFrame, tmParams } from "./logic/D03_PlanAnalysis_Tm"; /** 측점 간격(m) — B 규칙 측점과 같음. */ export const PLAN_STATION_INTERVAL_M = 20; @@ -48,6 +52,8 @@ export interface PlanMapScene { image: HTMLImageElement | null; /** 처음 맞출 점(시점 · 종점 · 그은 선). */ fit: Vertex[]; + /** 작업 좌표계 — TM 겹을 경위도 바탕에 맞춤(586-0). */ + epsg?: number | null; } /** 겹이 쓰는 좌표 — 그림 px(돌리기 전) ↔ 작업 좌표계 m · 글자 세우기. */ @@ -105,6 +111,8 @@ export interface PlanMap { fit: () => void; /** 검증 창구 — 작업 좌표계 점 → 화면에 보이는 px(돌린 뒤 · 그리기 전이면 null). */ screenOf: (point: Vertex) => [number, number] | null; + /** 검증 창구 — 경위도 → 화면 px(바탕 겹 꼴 · 586-0 맞춤 대조). */ + screenOfLonLat: (lon: number, lat: number) => [number, number] | null; /** 검증 창구 — 지금 돌린 각(도). */ rotationDeg: () => number; dispose: () => void; @@ -204,12 +212,40 @@ export function createPlanMap(events: PlanMapEvents): PlanMap { onChange: () => draw(), }); - const screenMap = (): ScreenMap => { - const projector = createMetricProjector(meta!, view); + /** TM ↔ 그림 px — 경위도 바탕에 맞춤(586-0) · 가짜 도엽 · 모르는 좌표계 = 옛 「같은 직사각형」 꼴. */ + const metricFrame = (): Pick => { + const params = tmParams(scene.epsg); + if (!params || !scene.meta || meta !== scene.meta) { + const projector = createMetricProjector(meta!, view); + return { + toScreen: (point) => projector.toScreen(point[0], point[1]), + toMetric: (px, py) => screenToMetric(meta!, view, px, py), + pxPerMeter: projector.pxPerMeter, + }; + } + const normalizer = createNormalizer(scene.meta); + const [x0, y0] = normalizedToScreen(view, 0, 0); + const [x1, y1] = normalizedToScreen(view, 1, 1); + const frame = localTmFrame( + params, + (lon, lat) => lonLatToScreen(normalizer, view, lon, lat), + (px, py) => [ + normalizer.lonMin + ((px - x0) / (x1 - x0)) * normalizer.lonRange, + normalizer.latMin + (1 - (py - y0) / (y1 - y0)) * normalizer.latRange, + ], + view.width / 2, + view.height / 2, + ); return { - toScreen: (point) => projector.toScreen(point[0], point[1]), - toMetric: (px, py) => screenToMetric(meta!, view, px, py), - pxPerMeter: projector.pxPerMeter, + toScreen: (point) => frame.toScreen(point[0], point[1]), + toMetric: (px, py) => frame.toMetric(px, py), + pxPerMeter: frame.pxPerMeter, + }; + }; + + const screenMap = (): ScreenMap => { + return { + ...metricFrame(), view, ...marks, uprightRad: rotation.uprightRad(), @@ -380,6 +416,10 @@ export function createPlanMap(events: PlanMapEvents): PlanMap { const [x, y] = screenMap().toScreen(point); return rotation.rerotate(x, y); }, + screenOfLonLat: (lon, lat) => { + if (!scene.meta) return null; + return rotation.rerotate(...lonLatToScreen(createNormalizer(scene.meta), view, lon, lat)); + }, rotationDeg: () => (rotation.radians() * 180) / Math.PI, dispose: () => { observer.disconnect(); diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_UI_Page.ts b/D03_PlanAnalysis/D03_PlanAnalysis_UI_Page.ts index 00e4952f6..9a9d25c5a 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_UI_Page.ts +++ b/D03_PlanAnalysis/D03_PlanAnalysis_UI_Page.ts @@ -577,6 +577,7 @@ export async function renderD03PlanAnalysis(root: HTMLElement): Promise { base: () => points, lineNames: () => lineNames(network.lines), screenOf: map.screenOf, + screenOfLonLat: map.screenOfLonLat, preview: draw.preview, stroke: draw.stroke, selected: edit.selected, @@ -632,7 +633,7 @@ export async function renderD03PlanAnalysis(root: HTMLElement): Promise { drawYearPick.setOptions(yearOptions(), String(drawYear)); const draft = readState("d03-plan", projectId); if (draft?.version === DRAFT_VERSION && draft.network?.lines?.length) network = draft.network; - map.update({ fit: [...points, ...allVertices()] }); + map.update({ fit: [...points, ...allVertices()], epsg }); refresh(); await loadSurface(); refresh(); diff --git a/D03_PlanAnalysis/logic/D03_PlanAnalysis_Tm.ts b/D03_PlanAnalysis/logic/D03_PlanAnalysis_Tm.ts new file mode 100644 index 000000000..1adbd6da8 --- /dev/null +++ b/D03_PlanAnalysis/logic/D03_PlanAnalysis_Tm.ts @@ -0,0 +1,138 @@ +/* ============================================================================= + * D03_PlanAnalysis_Tm.ts + * 작업 좌표계(TM m) ↔ 경위도 — 한국 GRS80 TM(EPSG:5185~5188 · 5179) · 화면 한 장 TM 틀(PLAN 586-0). + * 짝 = pyproj(`common_util_crs.tm_to_lonlat`) · 거울 시험 `test_586_tm_frame.py`. + * 바탕(위성 · 도엽 겹)은 경위도로 펴고 그은 선 · 급경사 그림은 TM 이라, 두 상자를 같은 직사각형으로 보면 + * 자오선 수렴각만큼 어긋남(D 프로젝트 귀퉁이 63 m) → TM 점을 경위도로 바꿔 바탕과 같은 자리에 놓음. + * DOM · fetch · sessionStorage 없음. + * ========================================================================== */ + +export interface TmParams { + lon0: number; + lat0: number; + k0: number; + fe: number; + fn: number; +} + +const KOREA_TM: Record = { + 5185: { lon0: 125, lat0: 38, k0: 1, fe: 200000, fn: 600000 }, + 5186: { lon0: 127, lat0: 38, k0: 1, fe: 200000, fn: 600000 }, + 5187: { lon0: 129, lat0: 38, k0: 1, fe: 200000, fn: 600000 }, + 5188: { lon0: 131, lat0: 38, k0: 1, fe: 200000, fn: 600000 }, + 5179: { lon0: 127.5, lat0: 38, k0: 0.9996, fe: 1000000, fn: 2000000 }, +}; + +/** 아는 좌표계만(GRS80 TM) — 그 밖(옛 Bessel 등) = null. */ +export const tmParams = (epsg: number | null | undefined): TmParams | null => + (epsg && KOREA_TM[epsg]) || null; + +// GRS80 · Krüger 6차(Karney 2011) — 수 km 안 mm 아래 +const F = 1 / 298.257222101; +const N = F / (2 - F); +const E = Math.sqrt(F * (2 - F)); +const A = (6378137 / (1 + N)) * (1 + N ** 2 / 4 + N ** 4 / 64 + N ** 6 / 256); +const ALPHA = [ + N / 2 - (2 * N ** 2) / 3 + (5 * N ** 3) / 16 + (41 * N ** 4) / 180 - (127 * N ** 5) / 288, + (13 * N ** 2) / 48 - (3 * N ** 3) / 5 + (557 * N ** 4) / 1440 + (281 * N ** 5) / 630, + (61 * N ** 3) / 240 - (103 * N ** 4) / 140 + (15061 * N ** 5) / 26880, + (49561 * N ** 4) / 161280 - (179 * N ** 5) / 168, + (34729 * N ** 5) / 80640, +]; +const BETA = [ + N / 2 - (2 * N ** 2) / 3 + (37 * N ** 3) / 96 - N ** 4 / 360 - (81 * N ** 5) / 512, + N ** 2 / 48 + N ** 3 / 15 - (437 * N ** 4) / 1440 + (46 * N ** 5) / 105, + (17 * N ** 3) / 480 - (37 * N ** 4) / 840 - (209 * N ** 5) / 4480, + (4397 * N ** 4) / 161280 - (11 * N ** 5) / 504, + (4583 * N ** 5) / 161280, +]; +const DELTA = [ + 2 * N - (2 * N ** 2) / 3 - 2 * N ** 3 + (116 * N ** 4) / 45, + (7 * N ** 2) / 3 - (8 * N ** 3) / 5 - (227 * N ** 4) / 45, + (56 * N ** 3) / 15 - (136 * N ** 4) / 35, + (4279 * N ** 4) / 630, +]; +const RAD = Math.PI / 180; + +/** 경위도 → 등각 (ξ, η). */ +function conformal(lonRad: number, latRad: number): [number, number] { + const s = Math.sin(latRad); + const t = Math.sinh(Math.atanh(s) - E * Math.atanh(E * s)); + const xi1 = Math.atan2(t, Math.cos(lonRad)); + const eta1 = Math.atanh(Math.sin(lonRad) / Math.sqrt(1 + t * t)); + let xi = xi1; + let eta = eta1; + ALPHA.forEach((a, i) => { + const j = 2 * (i + 1); + xi += a * Math.sin(j * xi1) * Math.cosh(j * eta1); + eta += a * Math.cos(j * xi1) * Math.sinh(j * eta1); + }); + return [xi, eta]; +} + +const xi0Of = (p: TmParams): number => conformal(0, p.lat0 * RAD)[0]; + +export function lonLatToTm(p: TmParams, lon: number, lat: number): [number, number] { + const [xi, eta] = conformal((lon - p.lon0) * RAD, lat * RAD); + return [p.fe + p.k0 * A * eta, p.fn + p.k0 * A * (xi - xi0Of(p))]; +} + +export function tmToLonLat(p: TmParams, x: number, y: number): [number, number] { + const xi = (y - p.fn) / (p.k0 * A) + xi0Of(p); + const eta = (x - p.fe) / (p.k0 * A); + let xi1 = xi; + let eta1 = eta; + BETA.forEach((b, i) => { + const j = 2 * (i + 1); + xi1 -= b * Math.sin(j * xi) * Math.cosh(j * eta); + eta1 -= b * Math.cos(j * xi) * Math.sinh(j * eta); + }); + const chi = Math.asin(Math.sin(xi1) / Math.cosh(eta1)); + let lat = chi; + DELTA.forEach((d, i) => (lat += d * Math.sin(2 * (i + 1) * chi))); + return [p.lon0 + Math.atan2(Math.sinh(eta1), Math.cos(xi1)) / RAD, lat / RAD]; +} + +/** 화면 한 장 TM 틀 — 화면 가운데 둘레 1차 맞춤(그 화면 안 1 px 아래 · 그림 겹은 세 귀퉁이로 놓음). */ +export interface TmFrame { + toScreen: (x: number, y: number) => [number, number]; + toMetric: (px: number, py: number) => [number, number]; + pxPerMeter: number; +} + +/** 미분 걸음(m) — 가운데 차분. */ +const STEP_M = 25; + +/** + * 바탕이 경위도를 화면에 펴는 길(`project` · `unproject`)을 받아, 화면 가운데(`cx`, `cy`) 둘레에서 + * TM → 화면을 1차식으로 맞춤. 바탕이 어떤 꼴(캔버스 직사각형 · MapLibre 메르카토르)이든 같은 식. + */ +export function localTmFrame( + p: TmParams, + project: (lon: number, lat: number) => [number, number], + unproject: (px: number, py: number) => [number, number], + cx: number, + cy: number, +): TmFrame { + const [mx, my] = lonLatToTm(p, ...unproject(cx, cy)); + const at = (x: number, y: number): [number, number] => project(...tmToLonLat(p, x, y)); + const [ox, oy] = at(mx, my); + const [x1, y1] = at(mx + STEP_M, my); + const [x0, y0] = at(mx - STEP_M, my); + const [x3, y3] = at(mx, my + STEP_M); + const [x2, y2] = at(mx, my - STEP_M); + const a = (x1 - x0) / (2 * STEP_M); + const b = (y1 - y0) / (2 * STEP_M); + const c = (x3 - x2) / (2 * STEP_M); + const d = (y3 - y2) / (2 * STEP_M); + const det = a * d - b * c || 1; + return { + toScreen: (x, y) => [ox + a * (x - mx) + c * (y - my), oy + b * (x - mx) + d * (y - my)], + toMetric: (px, py) => { + const dx = px - ox; + const dy = py - oy; + return [mx + (d * dx - c * dy) / det, my + (a * dy - b * dx) / det]; + }, + pxPerMeter: Math.sqrt(Math.abs(det)), + }; +} diff --git a/resources/tester/helper_586_tm_frame.cjs b/resources/tester/helper_586_tm_frame.cjs new file mode 100644 index 000000000..df8ed978c --- /dev/null +++ b/resources/tester/helper_586_tm_frame.cjs @@ -0,0 +1,41 @@ +/* 586 TM ↔ 경위도 · 화면 TM 틀 헬퍼 — `D03_PlanAnalysis_Tm.ts` 를 stdin 점 목록으로 돌려 JSON 으로 냄. + * 입력 { points: [[epsg, x, y]], lonlat: [[epsg, lon, lat]], frames: [{ epsg, box, view, probes }] } + */ +const path = require("path"); +const { makeLoader } = require("./helper_ts_stub_loader.cjs"); + +const ROOT = path.join(__dirname, "..", ".."); +const { load } = makeLoader({ root: ROOT }); +const tm = load("D03_PlanAnalysis", "logic", "D03_PlanAnalysis_Tm.ts"); + +const input = JSON.parse(require("fs").readFileSync(0, "utf-8")); +const out = { + unknown: tm.tmParams(5174), + toLonLat: input.points.map(([epsg, x, y]) => tm.tmToLonLat(tm.tmParams(epsg), x, y)), + toTm: input.lonlat.map(([epsg, lon, lat]) => tm.lonLatToTm(tm.tmParams(epsg), lon, lat)), + frames: input.frames.map(({ epsg, box, view, probes }) => { + // 경위도 상자를 화면 직사각형에 펴는 바탕(캔버스 지금 꼴) · 가운데 = (cx, cy) + const [lonMin, latMin, lonMax, latMax] = box; + const [w, h, x0, y0, sx, sy] = view; + const project = (lon, lat) => [ + x0 + ((lon - lonMin) / (lonMax - lonMin)) * sx, + y0 + ((latMax - lat) / (latMax - latMin)) * sy, + ]; + const unproject = (px, py) => [ + lonMin + ((px - x0) / sx) * (lonMax - lonMin), + latMax - ((py - y0) / sy) * (latMax - latMin), + ]; + const p = tm.tmParams(epsg); + const frame = tm.localTmFrame(p, project, unproject, w / 2, h / 2); + return { + pxPerMeter: frame.pxPerMeter, + probes: probes.map(([x, y]) => { + const exact = project(...tm.tmToLonLat(p, x, y)); + const got = frame.toScreen(x, y); + const back = frame.toMetric(...got); + return { exact, got, back }; + }), + }; + }), +}; +process.stdout.write(JSON.stringify(out)); diff --git a/resources/tester/test_582_contour_lod.py b/resources/tester/test_582_contour_lod.py index a6f76ac10..132a3dfca 100644 --- a/resources/tester/test_582_contour_lod.py +++ b/resources/tester/test_582_contour_lod.py @@ -75,7 +75,7 @@ def test_D03_만_켬(): layers = _read("D03_PlanAnalysis/D03_PlanAnalysis_UI_Layers.ts") assert "indexLayerLod(prepare(collection" in layers # 581-1 합침 — 겹 모두 바탕 캔버스 - assert layers.count("base: true") == layers.count("draw: (context, map)") == 12 + assert layers.count("base: true") == layers.count("draw: (context, map)") >= 12 # 겹 모두 바탕(581 · 584 로 늘어남) plan_map = _read("D03_PlanAnalysis/D03_PlanAnalysis_UI_Map.ts") assert "requestAnimationFrame(flush)" in plan_map assert "drawTop();" in plan_map diff --git a/resources/tester/test_586_tm_frame.py b/resources/tester/test_586_tm_frame.py new file mode 100644 index 000000000..6f9f14d88 --- /dev/null +++ b/resources/tester/test_586_tm_frame.py @@ -0,0 +1,120 @@ +# -*- coding: utf-8 -*- +"""586-0 TM ↔ 경위도 거울 · 화면 TM 틀 — 그은 선 · 급경사 그림을 위성 · 도엽 겹과 같은 자리에. + +- `D03_PlanAnalysis_Tm.ts` 의 TM ↔ 경위도 = pyproj(`common_util_crs`) 와 1 mm 안(5185~5188 · 5179) +- 모르는 좌표계(5174 Bessel) = null(지금 꼴로 그림) +- 화면 TM 틀(가운데 1차 맞춤) = 정식 변환 자리와 화면 안 0.5 px 안 · 되돌림 = 1 mm 안 +- 옛 「같은 직사각형」 꼴은 D 프로젝트 귀퉁이에서 수십 m 어긋남(고친 까닭) + +Node 헬퍼(`helper_586_tm_frame.cjs`)가 진짜 TS 를 돌린다. +""" + +import json +import math +import os +import subprocess + +import pytest + +from common_util.common_util_crs import lonlat_to_tm, tm_to_lonlat + +HERE = os.path.dirname(os.path.abspath(__file__)) + +# 검증_sub2_D 위성 상자(EPSG:5186 · vworld-meta) +META = { + "x_min": 282572.566652281, + "x_max": 288034.8885857916, + "y_min": 408543.914277056, + "y_max": 411930.1592677671, + "lon_min": 127.9193115234375, + "lon_max": 127.979736328125, + "lat_min": 36.270850207239036, + "lat_max": 36.30184530368433, +} +POINTS = [ + (5186, 285303.7, 410237.0), + (5186, 282572.6, 411930.2), + (5186, 288034.9, 408543.9), + (5185, 180000.0, 450000.0), + (5187, 210000.0, 350000.0), + (5188, 190000.0, 520000.0), + (5179, 960000.0, 1800000.0), +] +LONLAT = [(5186, 127.95, 36.28), (5187, 129.1, 35.5), (5179, 127.0, 37.5), (5185, 126.3, 34.8)] + + +def _view(w, h, scale, ox, oy): + """지금 캔버스 꼴 — 상자 비율 그대로 화면에 맞추고 배율 · 옮김.""" + ratio = (META["x_max"] - META["x_min"]) / (META["y_max"] - META["y_min"]) + mw, mh = (w, w / ratio) if ratio > w / h else (h * ratio, h) + x0 = (w - mw) / 2 * scale + w / 2 * (1 - scale) + ox + y0 = (h - mh) / 2 * scale + h / 2 * (1 - scale) + oy + return [w, h, x0, y0, mw * scale, mh * scale] + + +@pytest.fixture(scope="module") +def r(): + box = [META["lon_min"], META["lat_min"], META["lon_max"], META["lat_max"]] + corners = [ + (META["x_min"], META["y_max"]), + (META["x_max"], META["y_max"]), + (META["x_min"], META["y_min"]), + (META["x_max"], META["y_min"]), + (285303.7, 410237.0), + ] + frames = [ + {"epsg": 5186, "box": box, "view": _view(1400, 900, 1, 0, 0), "probes": corners}, + { + "epsg": 5186, + "box": box, + "view": _view(1400, 900, 40, -3000, 2000), + "probes": [(285303.7 + dx, 410237.0 + dy) for dx in (-40, 0, 40) for dy in (-25, 25)], + }, + ] + payload = {"points": POINTS, "lonlat": LONLAT, "frames": frames} + proc = subprocess.run( + ["node", os.path.join(HERE, "helper_586_tm_frame.cjs")], + input=json.dumps(payload), + capture_output=True, + text=True, + encoding="utf-8", + timeout=60, + ) + assert proc.returncode == 0, proc.stderr + return json.loads(proc.stdout) + + +def test_TM_경위도_pyproj_와_같음(r): + for (epsg, x, y), (lon, lat) in zip(POINTS, r["toLonLat"]): + want_lon, want_lat = tm_to_lonlat(x, y, epsg) + assert abs(lon - want_lon) < 1e-8 and abs(lat - want_lat) < 1e-8, (epsg, x, y) + for (epsg, lon, lat), (x, y) in zip(LONLAT, r["toTm"]): + want_x, want_y = lonlat_to_tm(lon, lat, epsg) + assert math.hypot(x - want_x, y - want_y) < 1e-3, (epsg, lon, lat) + + +def test_모르는_좌표계는_null(r): + assert r["unknown"] is None + + +def test_화면_TM_틀이_정식_자리와_같음(r): + for frame in r["frames"]: + for probe in frame["probes"]: + (ex, ey), (gx, gy) = probe["exact"], probe["got"] + assert math.hypot(ex - gx, ey - gy) < 0.5, probe + # 되돌림(화면 → TM) + probes = [(285303.7 + dx, 410237.0 + dy) for dx in (-40, 0, 40) for dy in (-25, 25)] + for (x, y), probe in zip(probes, r["frames"][1]["probes"]): + bx, by = probe["back"] + assert math.hypot(bx - x, by - y) < 1e-3 + + +def test_옛_같은_직사각형_꼴은_귀퉁이에서_어긋남(): + """고친 까닭 — TM 상자 · 경위도 상자를 같은 직사각형으로 보면 위 오른쪽 귀퉁이 수십 m.""" + x, y = META["x_max"], META["y_max"] + lon, lat = tm_to_lonlat(x, y, 5186) + w = META["x_max"] - META["x_min"] + h = META["y_max"] - META["y_min"] + dx = ((lon - META["lon_min"]) / (META["lon_max"] - META["lon_min"]) - 1) * w + dy = ((META["lat_max"] - lat) / (META["lat_max"] - META["lat_min"])) * h + assert math.hypot(dx, dy) > 50