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