Day 23 把管線串起來之後，偵測結果還停在像素座標，說「框在 (623, 386)」對誰都沒意義。今天要把 Day 15 那支從來沒跑起來過的 inference/geo_utils.py 接上管線，讓每個框都帶著經緯度出來。結果接完之後，前面三天的結論被推翻了一個。

第一撞：pyproj 連 import 都過不了

Day 15 的版本是這樣寫的：

import rasterio
from pyproj import Transformer

rasterio Day 23 已經換成 tifffile 了，pyproj 從來沒裝過，所以先補裝：

Successfully installed pyproj-3.8.0

裝是裝好了，一 import 就炸：

ImportError: DLL load failed while importing _compat: 應用程式控制原則已封鎖此檔案。

Day 23 我猜「pyproj 綁 PROJ，八成會跟 rasterio 一起中彈」，猜對了但猜錯了位置。翻事件記錄：

Code Integrity determined that a process (python.exe) attempted to load
...\.venv\Lib\site-packages\pyproj\_compat.cp312-win_amd64.pyd
that did not meet the Enterprise signing level requirements

被擋的不是 pyproj.libs\proj_9-....dll，是 pyproj 自己編譯的 _compat.pyd。也就是說它還沒走到載入 PROJ 那一步就先死了，Smart App Control 擋的是套件本體的 C extension。

順帶一提，pyproj.libs/ 裡面躺著跟 rasterio 同一套未簽署的 DLL（proj_9、jpeg62、sqlite3、tiff），就算 _compat.pyd 放行，下一關也還是會卡。

解法跟 Day 23 一樣：不關 Smart App Control。關掉是不可逆的，除非重灌 Windows 否則開不回來，為了一個座標轉換把整台機器的安全機制永久關掉，划不來。

所以今天要自己算 UTM 逆投影。

EPSG:32651 這串數字是什麼意思

在動公式之前先把代碼拆開。WGS84 的 UTM 分帶 EPSG 是有規律的：

326xx = WGS84 / UTM zone xx N（北半球）
327xx = WGS84 / UTM zone xx S（南半球）

所以 32651 = 北半球第 51 帶。每帶寬 6 度，第 1 帶的中央經線是 -177°，於是：

central_meridian = zone * 6 - 183      # 51 -> 123°E

南北半球的差別只有一個 false northing：南半球的 y 會加上 10,000,000 m，避免赤道以南出現負數。

def zone_from_epsg(epsg: int) -> UTMZone:
    if 32601 <= epsg <= 32660:
        return UTMZone(epsg - 32600, north=True)
    if 32701 <= epsg <= 32760:
        return UTMZone(epsg - 32700, north=False)
    raise ValueError(f"EPSG:{epsg} 不是 WGS84 UTM 分帶")

逆投影的難點：足點緯度

UTM 是橫麥卡托投影，正投影（經緯度 -> 平面）有封閉解，逆投影沒有。

卡住的地方是子午線弧長。從赤道沿著子午線走到緯度 φ 的距離 M(φ) 是一個橢圓積分，級數展開之後長這樣：

def _meridional_arc(phi):
    return A * (
        (1 - E2/4 - 3*E2**2/64 - 5*E2**3/256) * phi
        - (3*E2/8 + 3*E2**2/32 + 45*E2**3/1024) * np.sin(2*phi)
        + (15*E2**2/256 + 45*E2**3/1024) * np.sin(4*phi)
        - (35*E2**3/3072) * np.sin(6*phi)
    )

正向很好算，但逆投影要做的是已知 M 反求 φ，這個式子解不開。

標準做法是先假裝地球是球體，用 mu = M / (子午線平均半徑) 算出一個近似緯度，再用一組以扁率為參數的級數把它修正回橢球上，這個修正後的緯度叫足點緯度（footpoint latitude）：

e1 = (1 - np.sqrt(1 - E2)) / (1 + np.sqrt(1 - E2))
mu = m / (A * (1 - E2/4 - 3*E2**2/64 - 5*E2**3/256))
phi1 = (mu
        + (3*e1/2 - 27*e1**3/32) * np.sin(2*mu)
        + (21*e1**2/16 - 55*e1**4/32) * np.sin(4*mu)
        + (151*e1**3/96) * np.sin(6*mu)
        + (1097*e1**4/512) * np.sin(8*mu))

足點緯度的物理意義是：如果這個點在中央經線上，它會落在哪個緯度。有了它就能在該點求出兩個曲率半徑，再沿著垂直中央經線的方向展開級數，把真正的經緯度湊出來：

n1 = A / np.sqrt(1 - E2 * sin1**2)              # 卯酉圈曲率半徑
r1 = A * (1 - E2) / (1 - E2 * sin1**2) ** 1.5   # 子午圈曲率半徑
d = x / (n1 * K0)

lat = phi1 - (n1 * tan1 / r1) * (
    d**2/2
    - (5 + 3*t1 + 10*c1 - 4*c1**2 - 9*EP2) * d**4/24
    + (61 + 90*t1 + 298*c1 + 45*t1**2 - 252*EP2 - 3*c1**2) * d**6/720
)

公式抄自 Snyder 的《Map Projections: A Working Manual》，在分帶範圍（中央經線 ±3°）內誤差小於 1 mm。

這種扣要怎麼驗

這是今天最花時間的部分。沒有 pyproj 可以對答案，而這種一大串係數的級數展開，打錯一個 45 跟 55 根本看不出來，跑起來也不會報錯，只會默默偏掉幾十公尺。

所以我疊了三層驗證，一層比一層獨立。

第一層：往返一致

正投影再逆回來，看回不回得到原點。順手把正投影也寫了（反正之後畫地圖要用）：

2000 個隨機點，最大往返誤差 0.8275 mm

但這一層只能抓打錯字。如果我把同一個係數在正逆兩邊都寫錯，往返還是會完美閉合。這是自洽，不是正確。

第二層：子午線弧長用數值積分獨立驗

這一層才有意義。子午線弧長除了級數展開，還可以沿著緯度對子午圈曲率半徑做數值積分得到，兩條路的數學完全不同：

def radius(phi):
    return A * (1 - E2) / (1 - E2 * np.sin(phi)**2) ** 1.5

numeric = np.trapezoid(radius(np.linspace(0, phi, 200001)), ...)

最後一列還可以跟一個公開常數對答案，WGS84 從赤道到極點的子午線象限弧長是 10,001,965.729 m：

象限弧長 10001965.730 m  vs 公開值 10001965.729 m  差 0.5 mm

這個數字不是我算的，是大地測量學的標準值，對得上代表整組係數沒抄錯。

第三層：中央經線上的退化解

在中央經線上，橫麥卡托應該退化成一個很單純的形式，easting 恰好等於 500000，northing 恰好等於 k0 x M(φ)。這是一個不需要任何外部資料就能檢查的硬約束：

lat=  0.0  E= 500000.0000  N=       0.0000  誤差 0.0000 mm
lat= 25.0  E= 500000.0000  N= 2764947.7475  誤差 0.0000 mm
lat= 60.0  E= 500000.0000  N= 6651411.1911  誤差 0.0000 mm

三層都過了，這才敢拿去用。

地面真相

最後把 GeoTIFF 的四個角丟進去：

EPSG:32651 -> UTM 51N  中央經線 123°E
左上  (   0.0,   0.0) -> 25.123231N, 121.165839E
右上  (1349.0,   0.0) -> 25.124826N, 121.299593E
左下  (   0.0, 923.0) -> 25.039918N, 121.167080E
右下  (1349.0, 923.0) -> 25.041507N, 121.300744E
中心  ( 674.5, 461.5) -> 25.082386N, 121.233313E

檔名叫 som_san_air.tif，我一直以為是索馬利亞。25.08N 121.23E 是桃園，這張圖從頭到尾都是桃園國際機場一帶。

接上管線：GeoJSON

有了經緯度就可以吐標準格式。RFC 7946 有三個坑：

一、座標順序是 (經度, 緯度)，跟中文講「經緯度」的順序一樣，但跟大家報座標習慣講的「25.08N, 121.23E」相反。寫錯的話點會跑到索馬利亞外海，而且不會有任何錯誤訊息。

二、多邊形的環必須閉合，第一個點要重複放在最後一個。

三、外環必須是逆時針（右手定則）。這一條在我們的情境下特別容易錯，因為：

# 像素座標的 row 往下增加、經緯度的 lat 往上增加，
# 換算之後繞行方向會`翻面`，不能沿用像素空間的順序
area = sum(ring[i][0]*ring[i+1][1] - ring[i+1][0]*ring[i][1]
           for i in range(len(ring)-1))
return ring if area > 0 else ring[::-1]

另外畫地圖要用到反向轉換（經緯度 -> 像素），這裡也有個小坑，GeoRef 存的仿射矩陣是 2x3 不是方陣，沒辦法直接求逆，要補一列 [0, 0, 1] 湊成 3x3 齊次矩陣：

m = np.vstack([georef.transform, [0.0, 0.0, 1.0]])
inv = np.linalg.inv(m)

輸出長這樣：

{
  "type": "Feature",
  "geometry": {"type": "Polygon", "coordinates": [[
    [121.249044, 25.099070], [121.243307, 25.105579],
    [121.207183, 25.079172], [121.212920, 25.072665],
    [121.249044, 25.099070]]]},
  "properties": {
    "score": 0.5067, "n_fragments": 2, "aspect": 5.05,
    "center": [121.228112, 25.089122],
    "gmaps": "https://www.google.com/maps/search/?api=1&query=25.089122,121.228112"
  }
}

換成公尺之後，事情變得很難看

這才是今天真正的重點。以前框的尺寸只能講像素，現在可以講公尺了：

[0] conf=0.507 碎片x2
    中心 25.089122N, 121.228112E
    實地尺寸 4,675 m x 926 m （像素 467 x 92）
    框住面積 4.32 km²

一條跑道寬 926 公尺。

桃園機場 05L/23R 的實際規格是 3,660 m x 60 m。我們的框寬了15 倍。

Day 22 我寫「長寬比從 2.52 拉到 5.05，形狀終於像一條跑道」，那是在像素空間裡看圖說故事。換成公尺，5.05 的長寬比對應的是一個4.7 公里長、將近 1 公里寬的東西，那不是跑道，那是一整片工業區。

地面控制點：框根本不在機場上

既然有經緯度了，就可以拿公開座標來驗。桃園國際機場（RCTP）的機場基準點 ARP 是 ICAO 公告的 25°04'36"N 121°13'58"E，把它反投影回像素：

TPE ARP 25.0767N 121.2328E -> 像素 (668, 524)
偵測中心 25.0891N 121.2281E -> 像素 (623, 386)

灰色是 0.02 度一格的經緯網格，青色十字是 ARP，綠色是模型的框，黃色是框的中心。放大看：

青色十字精準地落在跑道上，這代表投影算對了。但綠色的框：

ARP 在框內?            False
ARP 到框邊界的距離      907 m
框中心到 ARP 的距離   1,455 m

機場基準點根本不在框裡面，差了快一公里。

框的方向是對的，跟跑道平行；框的長度也還行；但它整個往西北偏了 1.4 公里，壓在跑道北側的儲槽區與工業區上，只有右下那條長邊擦到機場的西北邊界。

Day 22 我寫「一個框、一條跑道、零誤判」。加上地面控制點之後，正確的說法是一個框、零命中。

回頭看標註：我們一直在往錯的方向優化

為什麼模型會框出一個 4.7 km x 0.9 km 的東西？回去看資料集的類別名稱：

names:
- airport
nc: 1

是 airport 不是 runway。標註框的本來就是整座機場，不是一條跑道。

把訓練集的標註全部丟進 minAreaRect 量一次形狀：

標註的長寬比中位數是 1.86，是一個接近正方形的區塊。

這下 Day 21 到 Day 22 那條「長寬比越高越好」的路線整個站不住：

• Day 21 幾何合併把長寬比從 1.54 拉到 2.52，我當時嫌還不夠像跑道
• Day 22 混合資料集拉到 5.05，我寫形狀終於像一條跑道

但標註教模型的目標是 1.86。Day 21 的 2.52 其實比 Day 22 的 5.05 更接近標註。我連續兩天把「偏離訓練目標」當成進步在慶祝，因為我腦中想的是跑道，而資料集講的是機場。

這跟 Day 23 的資料洩漏是同一類錯誤：沒有先把量測對象搞清楚就開始優化。

小結

今天本來只是要把像素換成經緯度，結果換出三件事。

pyproj 跟 rasterio 一起被 Smart App Control 擋掉，而且擋在更前面的 _compat.pyd，所以 UTM 逆投影只能自己寫。寫這種級數展開最大的風險不是跑不動而是默默算錯，所以驗證疊了三層：往返只能抓打錯字，子午線弧長用數值積分獨立算一次才是真的在驗，最後對上公開的象限弧長常數 10,001,965.729 m，差 0.5 mm。

接上管線吐 GeoJSON，RFC 7946 的三個坑（座標順序、環要閉合、外環逆時針）裡最陰的是第三個，因為像素的 row 往下、緯度往上，繞行方向會翻面。

但真正的收穫是換算成公尺之後看到的東西。我們引以為傲的那個框是 4,675 m x 926 m，用公開的機場基準點一驗，ARP 根本不在框內，差了 907 公尺。再回頭看標註才發現類別叫 airport，長寬比中位數只有 1.86，而我花了兩天把它推到 5.05 還很高興。

一個座標系統就把三天的自我感覺良好戳破了，這大概是今天最值得記下來的一件事。

明天要把 GeoJSON 丟上 folium 做成可以互動的網頁地圖，順便把今天這張經緯網格換成真的底圖，看看那個框到底壓在什麼東西上面，那我們明天見。