我寫了一個不依賴 GDAL 的 FileGDB 讀寫庫然后發(fā)現(xiàn) GDAL 有一半的活是 GEOS 干的項目地址https://github.com/lizhiziwang/pyopengdbpyopenfilegdb—— 純 Python 讀寫 Esri File Geodatabase零運行時依賴pip install下去不用帶一個 DLL。一、為什么想甩掉 GDAL做過交付的人大概都經(jīng)歷過這個循環(huán)本地跑得好好的腳本拷到客戶機器上就死了。查半天是 GDAL。裝 GDAL 要綁 OSGeo4W / conda版本還得跟 Python 對上PyInstaller 打包一個.gdb的處理腳本產(chǎn)物體積從 20 MB 漲到 500 MB最要命的是你只想讀三個字段卻拖進來一整套投影、柵格、WMS 的機器。而 GDAL 本身沒有錯。它是地理界的 libc沒有它整個行業(yè)轉(zhuǎn)不動。問題只是它的粒度是一整個地理引擎而我的需求是讀一個目錄里的二進制文件。所以我給自己定了個目標(biāo)寫一個庫能讀寫.gdb只用標(biāo)準(zhǔn)庫不 importosgeo、不 importfiona、不調(diào)任何外部二進制。聽起來像不可能——因為 FileGDB 是沒有公開格式規(guī)范的私有格式。二、關(guān)鍵的轉(zhuǎn)向把 GDAL 當(dāng)規(guī)范讀這是整個項目成立的前提。Esri 不公開 FileGDB 的格式規(guī)范。全世界范圍內(nèi)對.gdb二進制格式最完整、最權(quán)威、最可讀的公開實現(xiàn)就是 GDAL 的ogr/ogrsf_frmts/openfilegdb/。所以不依賴 GDAL的正確姿勢不是繞過GDAL而是讀懂GDAL本庫的全部格式知識逆向自 GDAL 的 openfilegdb 驅(qū)動 C 源碼每一處格式細節(jié)在代碼注釋里都標(biāo)了對應(yīng)的 GDAL 出處。這帶來一個很舒服的性質(zhì)不需要猜。.gdbtable頭部第幾個字節(jié)是什么、.gdbtablx的 offset 寬度有 4/5/6 三種怎么判、Esri 的 varint 為什么這么編碼、環(huán)的繞向怎么存——GDAL 都寫清楚了照抄即可。建出來的庫可以直接被 ArcGIS 和 QGIS 打開七張系統(tǒng)表與真實 ArcGIS空白庫逐段字節(jié)一致有測試守著。frompyopenfilegdbimportOpenFileGDBwithOpenFileGDB.open(xxxxx)asgdb:layergdb.get_layer(xxxxx)print(layer.record_count)# 21217forfeatinlayer.read_features(limit3):print(feat.oid,feat.attributes[xxxxx],feat.geometry.area())這一檔嚴(yán)格按 GDAL 實現(xiàn)是完全成立的——因為 GDAL 的 openfilegdb驅(qū)動本身就是純 C它自己都不依賴 GEOS。但這個成立有個邊界而我是在踩上去之后才發(fā)現(xiàn)邊界在哪。三、天花板GDAL 有一半的活其實是轉(zhuǎn)發(fā)給 GEOS 的任務(wù)往下走需求變成再給我?guī)讉€空間計算面積、質(zhì)心、凸包、相交、包含。按之前的思路我應(yīng)該去 GDAL 里找對應(yīng)的 C 抄。結(jié)果打開ogrgeometry.cpp看到的是這個形狀下面是大意不是逐字引文intOGRGeometry::Intersects(constOGRGeometry*poThis,constOGRGeometry*poOther){#ifdefHAVE_GEOS// 包圍盒快速排除剩下的……全部轉(zhuǎn)手給 GEOS#endif}Intersects/Contains/Touches/Within/Covers…… 全部是這個形狀一句轉(zhuǎn)發(fā)外加一句包圍盒預(yù)篩。GDAL 自己加的唯一東西是包圍盒先快速排除。OGRGeometry::Buffer也一樣一行算法都沒有。也就是說在空間計算這一層根本沒有GDAL 的實現(xiàn)可以照抄。這是這個項目里我最有收獲的一課說嚴(yán)格按 GDAL 實現(xiàn)之前先分清你要抄的是格式還是謂詞。于是參照系必須換。GDAL 的轉(zhuǎn)發(fā)目標(biāo)是GEOS而 GEOS 是JTSJava Topology Suite的 C 移植。所以這一層的復(fù)刻對象變成了 JTS功能復(fù)刻對象對拍工具DE-9IM 拓撲謂詞OGC 規(guī)范 / GEOStools/verify_topology.pybuffer()JTSoperation/buffer/*tools/verify_buffer.pydifference/union/intersection/symmetric_differenceJTSoperation/overlayng/*tools/verify_overlay.py三者共用一套自實現(xiàn)的平面圖_planar.py對拍結(jié)果verify_buffer.py buffer vs GEOS 1096 對 · 逐位相同 1015 (92.6%) · 超差 0 verify_overlay.py overlay vs GEOS 1264 對 · 集合結(jié)構(gòu)全同 1220 · 超差 0那幾十對排除的都是能說清原因的閉合輸入上的單側(cè)緩沖GEOS 后面還有一步OverlayNGPolygonizer取最大面以及線結(jié)果 JTS 口徑是逐條節(jié)點邊輸出、GEOS 會合并成最長鏈。都寫在文檔里不是跑不過就跳過。四、寫拓撲算法最貴的一課真值表攔不住 bug這一段是整個項目里我覺得最值錢的工程經(jīng)驗值得單獨講。拓撲謂詞的實現(xiàn)最自然的驗證方式是手推真值表兩個矩形相交答案應(yīng)該是212101212點在線端點答案應(yīng)該是FF1FF0102……我寫了這些全綠。然后它們?nèi)紱]抓到真正的 bug。原因是手推真值表用的坐標(biāo)全是整數(shù)整數(shù)坐標(biāo)下所有算出來的中間點都逐位精確。而真實數(shù)據(jù)的坐標(biāo)不是。真正的 bug 有兩個同一個根因一句話重算出來的點不能拿去問浮點。relate()的實現(xiàn)里需要一堆探針點判斷拓撲角色——交點、中點、子段采樣點。它們是算出來的一般情況下不會精確落在對方線段上。而在不在線上用的是orient(...) 0的精確比較。于是在線上被判成在內(nèi)部/外部整塊 DE-9IM 矩陣跟著塌。最小復(fù)現(xiàn)都是極普通的情形不是人造的退化構(gòu)型# 1) 交點坐標(biāo)除不盡LINESTRING(-0.00000030.0000004,10.00000029.9999998).relate(LINESTRING(010,100))# 老實現(xiàn) FF1FF0102 GEOS 0F1FF0102# 2) 共線重疊子段的中點 —— 注意是拿它跟它自己比POLYGON((00,0.10.3,0.40.1,00)).relate(它自己)# 老實現(xiàn) 2F2F11212 正確 2FFF1FFF2bug 2 在真實數(shù)據(jù)上是能傷人的共線重疊在行政區(qū)劃數(shù)據(jù)里到處都是相鄰行政區(qū)共享界線它把2024年國土行政區(qū)劃.gdb里兩條完全相同的鄉(xiāng)/村面判成了部分重疊。而且 bug 2 是修完 bug 1 之后才浮出來的還是把對拍抽樣從8 條要素 / 400 對放大到150 條要素 / 900 對才撞出來的。所以還有第二條教訓(xùn)對拍跑過且綠這句話必須帶上參數(shù)否則它沒有意義。更陰的一課判定口徑本身也可能錯寫 overlay 對拍時我拿兩個結(jié)果的對稱差線長當(dāng)集合是否相同的主判據(jù)。920 對里報了 15 對集合不同。逐個查下來全是次 ULP 的偽差同一個環(huán)把一個頂點的 y 動了 1 ulp5.68e-14兩環(huán)周長只差 6.75e-14——而 GEOS 自己給出的symmetric_difference卻是一個周長 4.75的退化環(huán)。是 GEOS 對近重合輸入的數(shù)值不穩(wěn)。于是判定改成三條腿對稱差面積兩個結(jié)果各自的總長之差Hausdorff 距離。換完還逐條注入缺陷實測每條腿各管什么點被挪 0.4 → 面積和線長都是恒 0只有 Hausdorff 抓得住線多一根 0.4 的刺 → 長度差抓面多出一塊 0.16 → 面積差抓。這直接引出一個反直覺的結(jié)論用 GDAL 當(dāng)對拍 oracle 時是沒有 Hausdorff 的C 層沒暴露被挪開的點會靜默通過。所以那個工具會明說自己少了這條腿而不是假裝驗過了。判定口徑要有它有活干的證明不然它只是看上去在守。五、性能純 Python 打 C靠的不是寫快是別干活基準(zhǔn)語料村行政區(qū)劃圖層21,217 條面要素、4,494 萬個頂點、記錄體 261.7 MB。場景本庫對照只讀屬性不解幾何0.42 sGDAL/OGR 同口徑0.3 ~ 0.742 s只要包圍盒0.892 s——連幾何一起讀1.50 sGDAL/OGR 同口徑0.80 s→1.9×連幾何一起讀強制純 Python不開 C 擴展65.4 s比開 C 擴展慢44×關(guān)鍵是第一行和第三行的差距為什么這么大。把 4,494 萬個頂點的 varint 解碼用 C 寫一遍能拿到 116× 的單點提速——但整層也只從 59.3 s 降到 4.8 s。真正的正解是惰性幾何。read_features()出來的要素只帶幾何 blob 的未解碼引用第一次訪問feat.geometry才解析并緩存。所以只讀屬性的循環(huán)一個頂點都不解—— 這正是 GDALSetIgnoredFields省掉的那一段。連幾何的字節(jié)都不讀。幾何占記錄體的99.1%。先用 32 字節(jié)探針探出幾何位置再只讀幾何前后兩段中間那段留個文件偏移給惰性解析。判據(jù)是跳讀讀得更少才跳所以讀的字節(jié)數(shù)永不超過優(yōu)化前。??但這條優(yōu)化有前提得說清楚下游要是逐點用坐標(biāo)for x, y in g.coordinates:省下的時間原樣回來。前提是下游按數(shù)組用坐標(biāo)。我誠實給出 GDAL 贏的那一欄連幾何一起讀1.9×。那是 C 語言的勝利不是算法的失敗——本庫那段解碼循環(huán)就是 GDAL 算法的逐位復(fù)刻GDAL masterfilegdbtable.cpp的ReadVarIntAndAddNoCheck只有一個 1 字節(jié)早退ReadXYArray就是個普通模板循環(huán)。寫路徑從 297× 說起寫路徑原先有個真 bug每條write_feature()都重寫整份.gdbtablx10,000 條要33.5 s3.35 ms/條而且N 越大越壞。改法不是批量落盤而是照 GDAL 的邏輯GDAL 是每條都寫但每條只寫 O(1)只有頭部賬目是懶的。改完后條數(shù)之前之后ms/條1,0001.654 s0.013 s0.0134,0008.879 s0.047 s0.01210,00033.541 s0.113 s0.011耗時與條數(shù)成正比了同一條代碼路徑297×提速。寫法上沒有花招直接write_feature循環(huán)最后layer.sync()一次。一個會讓你白測一天的坑??性能必須在終端里測不能按 PyCharm 的 Debug 跑。PyCharm 的調(diào)試器pydevd給每個 Python 幀裝 line tracer每條字節(jié)碼回調(diào)一次。本庫是純 Python全額上稅GDAL 是 C 擴展tracer 進不去稅率為 0。在調(diào)試器里比較純 Python 實現(xiàn) vs C 實現(xiàn)量到的是 tracer 的稅率差不是算法差。六、可選加速器裝上就快沒裝照樣全功能三個可選依賴性質(zhì)不一樣這個區(qū)分很重要沒有它影響numpy自動回退純 Python大環(huán)上面積/周長/質(zhì)心慢×21~25C 擴展_gdbaccel.c自動回退純 Python整層幾何解碼慢 ~9×pyproj沒有就沒有這個功能to_crs()報 ImportError前兩個是有就快些dependencies里仍然是零pyproj 是沒有這個功能所以單獨占一個 extra[crs]。兩個加速器都有差分閘門守著拿兩個實現(xiàn)互相對拍逐位相同才算過python tools/verify_accel.py# C vs 純 Python:8 個庫 11,885 條 / 6,335,486 頂點python tools/verify_numpy.py# numpy 快路徑 vs 順序路徑 Fraction 精確解校準(zhǔn)?? 一個具體的坑.pyd綁解釋器版本cp311/cp313各一份。只編了一份卻用另一個解釋器跑會HAS_ACCEL False靜默回退——實測因此從 4.4 s 變成 56.5 s。哪個解釋器跑代碼就用哪個編。七、我誠實列一下什么情況不該用它一個庫的價值一半在于它拒絕做什么。以下都是明說的邊界不是 bug能替換 GDAL 的部分讀寫.gdbversion 3 / ArcGIS 10.xversion 4 可讀建庫、建圖層、增刪改要素產(chǎn)出能被 ArcGIS/QGIS 直接打開度量、構(gòu)造、DE-9IM 謂詞、buffer、overlay 四算子WKT / GeoJSON 出口__geo_interface__可直接喂 geopandas/folium不能替換 GDAL 的部分請繼續(xù)用 GDAL任何非 FileGDB格式。這個庫只認.gdb沒有 Shapefile、GeoPackage、GeoJSON 驅(qū)動也沒有柵格。需要精確拓撲的場景。本庫用浮點方向判定無 snap-rounding、無精確算術(shù)。已知的一處1e-9 高的薄片平移到 UTM 量級后只有 2.1 ULP 高雙精度網(wǎng)格上根本沒有內(nèi)部點可采樣。GEOS 靠組合式拓撲圖躲過采樣式實現(xiàn)躲不過。is_valid()/is_simple()是部分實現(xiàn)只查環(huán)閉合/頂點數(shù)下限/單環(huán)自交/洞在殼內(nèi)。is_valid() True不等于OGC 有效。distance()是O(n·m)無線段包圍盒預(yù)篩。寫回不是逐位冪等全語料 21,217 條讀→寫→讀152 條0.7%字節(jié)不同環(huán)順序/繞向規(guī)范化 ≤1 個量化步長的舍入。面積全部一致是表示差異不是幾何差異——但要做哈希去重的場景不能拿本庫的輸出當(dāng)原文。buffer/ overlay 結(jié)果恒為 2DZ/M 丟棄環(huán)是Esri 繞向GEOMETRYCOLLECTION是純內(nèi)存類型、寫不進 .gdb。八、驗證情況383 個用例四配置全綠默認 / 關(guān) numpy / 關(guān) C 擴展 / 換解釋器——改一個全局名字最容易漏掉某個引用點所以每條都過。庫本體19,143 行18,548 行.py 595 行.c測試6,192 行。差分閘門拿兩個實現(xiàn)互相對拍逐位相同才算過python tools/verify_accel.py# C vs 純 Pythonpython tools/verify_numpy.py# numpy vs 順序路徑python tools/verify_topology.py# relate() 10 個謂詞 vs GEOSpython tools/verify_buffer.py# buffer vs GEOSpython tools/verify_overlay.py# 4 個 overlay 算子 vs GEOSpython tools/verify_wkt_roundtrip.py# 全語料 WKT 出口往返體檢python tools/bench_read.pyD:/work/x.gdb村行政區(qū)劃import shapely/import osgeo只出現(xiàn)在tools/里不進包、不進tests/、不進pyproject.toml。沒裝就 skip 并exit 0。九、現(xiàn)在就想試pipinstall.# 純 Python零運行時依賴pipinstall.[all]# 附帶 numpy pyprojfrompyopenfilegdbimportOpenFileGDB,GdbField,Geometry,FGFT_STRING,FGFT_DOUBLE# 讀withOpenFileGDB.open(D:/data/x.gdb)asgdb:layergdb.get_layer(村行政區(qū)劃)featslayer.read_features(where縣名 朝陽區(qū),limit100)# 寫gdbOpenFileGDB.create(D:/data/new.gdb)layergdb.create_layer(監(jiān)測點,geometry_typepoint,fields[GdbField(NAME,FGFT_STRING,length64),GdbField(HEIGHT,FGFT_DOUBLE),])layer.write_feature({NAME:A1,HEIGHT:12.5,Shape:Geometry.from_wkt(POINT(116.4 39.9))})gdb.close()# ← 關(guān)閉 結(jié)賬必須MIT 許可。最后這個項目讓我改了一個習(xí)慣在說按 X 實現(xiàn)之前先確認 X 真的實現(xiàn)了它。GDAL 在格式上是一座金礦——文檔沒有代碼就是規(guī)范逐字節(jié)照抄就能得到一個能被 ArcGIS 打開的庫。但在空間計算上它是一張貼在 GEOS 上的標(biāo)簽紙。把這兩件事分清楚才知道什么時候可以照抄、什么時候必須換參照系、什么時候必須自己對著規(guī)范從零寫并拿真參照實現(xiàn)對拍。以及那個我最想記住的真值表、不變量、假 oracle 三樣都攔不住真正的 bug。只有跑真參照實現(xiàn)才抓得到。倉庫https://github.com/lizhiziwang/pyopengdb · 詳見README.md· 可選 C 擴展編譯見ACCEL.md