我写了一个不依赖 GDAL 的 FileGDB 读写库然后发现 GDAL 有一半的活是 GEOS 干的项目地址https://github.com/lizhiziwang/pyopengdbpyopenfilegdb—— 纯 Python 读写 Esri File Geodatabase零运行时依赖pip install下去不用带一个 DLL。一、为什么想甩掉 GDAL做过交付的人大概都经历过这个循环本地跑得好好的脚本拷到客户机器上就死了。查半天是 GDAL。装 GDAL 要绑 OSGeo4W / conda版本还得跟 Python 对上PyInstaller 打包一个.gdb的处理脚本产物体积从 20 MB 涨到 500 MB最要命的是你只想读三个字段却拖进来一整套投影、栅格、WMS 的机器。而 GDAL 本身没有错。它是地理界的 libc没有它整个行业转不动。问题只是它的粒度是一整个地理引擎而我的需求是读一个目录里的二进制文件。所以我给自己定了个目标写一个库能读写.gdb只用标准库不 importosgeo、不 importfiona、不调任何外部二进制。听起来像不可能——因为 FileGDB 是没有公开格式规范的私有格式。二、关键的转向把 GDAL 当规范读这是整个项目成立的前提。Esri 不公开 FileGDB 的格式规范。全世界范围内对.gdb二进制格式最完整、最权威、最可读的公开实现就是 GDAL 的ogr/ogrsf_frmts/openfilegdb/。所以不依赖 GDAL的正确姿势不是绕过GDAL而是读懂GDAL本库的全部格式知识逆向自 GDAL 的 openfilegdb 驱动 C 源码每一处格式细节在代码注释里都标了对应的 GDAL 出处。这带来一个很舒服的性质不需要猜。.gdbtable头部第几个字节是什么、.gdbtablx的 offset 宽度有 4/5/6 三种怎么判、Esri 的 varint 为什么这么编码、环的绕向怎么存——GDAL 都写清楚了照抄即可。建出来的库可以直接被 ArcGIS 和 QGIS 打开七张系统表与真实 ArcGIS空白库逐段字节一致有测试守着。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())这一档严格按 GDAL 实现是完全成立的——因为 GDAL 的 openfilegdb驱动本身就是纯 C它自己都不依赖 GEOS。但这个成立有个边界而我是在踩上去之后才发现边界在哪。三、天花板GDAL 有一半的活其实是转发给 GEOS 的任务往下走需求变成再给我几个空间计算面积、质心、凸包、相交、包含。按之前的思路我应该去 GDAL 里找对应的 C 抄。结果打开ogrgeometry.cpp看到的是这个形状下面是大意不是逐字引文intOGRGeometry::Intersects(constOGRGeometry*poThis,constOGRGeometry*poOther){#ifdefHAVE_GEOS// 包围盒快速排除剩下的……全部转手给 GEOS#endif}Intersects/Contains/Touches/Within/Covers…… 全部是这个形状一句转发外加一句包围盒预筛。GDAL 自己加的唯一东西是包围盒先快速排除。OGRGeometry::Buffer也一样一行算法都没有。也就是说在空间计算这一层根本没有GDAL 的实现可以照抄。这是这个项目里我最有收获的一课说严格按 GDAL 实现之前先分清你要抄的是格式还是谓词。于是参照系必须换。GDAL 的转发目标是GEOS而 GEOS 是JTSJava Topology Suite的 C 移植。所以这一层的复刻对象变成了 JTS功能复刻对象对拍工具DE-9IM 拓扑谓词OGC 规范 / GEOStools/verify_topology.pybuffer()JTSoperation/buffer/*tools/verify_buffer.pydifference/union/intersection/symmetric_differenceJTSoperation/overlayng/*tools/verify_overlay.py三者共用一套自实现的平面图_planar.py对拍结果verify_buffer.py buffer vs GEOS 1096 对 · 逐位相同 1015 (92.6%) · 超差 0 verify_overlay.py overlay vs GEOS 1264 对 · 集合结构全同 1220 · 超差 0那几十对排除的都是能说清原因的闭合输入上的单侧缓冲GEOS 后面还有一步OverlayNGPolygonizer取最大面以及线结果 JTS 口径是逐条节点边输出、GEOS 会合并成最长链。都写在文档里不是跑不过就跳过。四、写拓扑算法最贵的一课真值表拦不住 bug这一段是整个项目里我觉得最值钱的工程经验值得单独讲。拓扑谓词的实现最自然的验证方式是手推真值表两个矩形相交答案应该是212101212点在线端点答案应该是FF1FF0102……我写了这些全绿。然后它们全都没抓到真正的 bug。原因是手推真值表用的坐标全是整数整数坐标下所有算出来的中间点都逐位精确。而真实数据的坐标不是。真正的 bug 有两个同一个根因一句话重算出来的点不能拿去问浮点。relate()的实现里需要一堆探针点判断拓扑角色——交点、中点、子段采样点。它们是算出来的一般情况下不会精确落在对方线段上。而在不在线上用的是orient(...) 0的精确比较。于是在线上被判成在内部/外部整块 DE-9IM 矩阵跟着塌。最小复现都是极普通的情形不是人造的退化构型# 1) 交点坐标除不尽LINESTRING(-0.00000030.0000004,10.00000029.9999998).relate(LINESTRING(010,100))# 老实现 FF1FF0102 GEOS 0F1FF0102# 2) 共线重叠子段的中点 —— 注意是拿它跟它自己比POLYGON((00,0.10.3,0.40.1,00)).relate(它自己)# 老实现 2F2F11212 正确 2FFF1FFF2bug 2 在真实数据上是能伤人的共线重叠在行政区划数据里到处都是相邻行政区共享界线它把2024年国土行政区划.gdb里两条完全相同的乡/村面判成了部分重叠。而且 bug 2 是修完 bug 1 之后才浮出来的还是把对拍抽样从8 条要素 / 400 对放大到150 条要素 / 900 对才撞出来的。所以还有第二条教训对拍跑过且绿这句话必须带上参数否则它没有意义。更阴的一课判定口径本身也可能错写 overlay 对拍时我拿两个结果的对称差线长当集合是否相同的主判据。920 对里报了 15 对集合不同。逐个查下来全是次 ULP 的伪差同一个环把一个顶点的 y 动了 1 ulp5.68e-14两环周长只差 6.75e-14——而 GEOS 自己给出的symmetric_difference却是一个周长 4.75的退化环。是 GEOS 对近重合输入的数值不稳。于是判定改成三条腿对称差面积两个结果各自的总长之差Hausdorff 距离。换完还逐条注入缺陷实测每条腿各管什么点被挪 0.4 → 面积和线长都是恒 0只有 Hausdorff 抓得住线多一根 0.4 的刺 → 长度差抓面多出一块 0.16 → 面积差抓。这直接引出一个反直觉的结论用 GDAL 当对拍 oracle 时是没有 Hausdorff 的C 层没暴露被挪开的点会静默通过。所以那个工具会明说自己少了这条腿而不是假装验过了。判定口径要有它有活干的证明不然它只是看上去在守。五、性能纯 Python 打 C靠的不是写快是别干活基准语料村行政区划图层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×关键是第一行和第三行的差距为什么这么大。把 4,494 万个顶点的 varint 解码用 C 写一遍能拿到 116× 的单点提速——但整层也只从 59.3 s 降到 4.8 s。真正的正解是惰性几何。read_features()出来的要素只带几何 blob 的未解码引用第一次访问feat.geometry才解析并缓存。所以只读属性的循环一个顶点都不解—— 这正是 GDALSetIgnoredFields省掉的那一段。连几何的字节都不读。几何占记录体的99.1%。先用 32 字节探针探出几何位置再只读几何前后两段中间那段留个文件偏移给惰性解析。判据是跳读读得更少才跳所以读的字节数永不超过优化前。⚠️但这条优化有前提得说清楚下游要是逐点用坐标for x, y in g.coordinates:省下的时间原样回来。前提是下游按数组用坐标。我诚实给出 GDAL 赢的那一栏连几何一起读1.9×。那是 C 语言的胜利不是算法的失败——本库那段解码循环就是 GDAL 算法的逐位复刻GDAL masterfilegdbtable.cpp的ReadVarIntAndAddNoCheck只有一个 1 字节早退ReadXYArray就是个普通模板循环。写路径从 297× 说起写路径原先有个真 bug每条write_feature()都重写整份.gdbtablx10,000 条要33.5 s3.35 ms/条而且N 越大越坏。改法不是批量落盘而是照 GDAL 的逻辑GDAL 是每条都写但每条只写 O(1)只有头部账目是懒的。改完后条数之前之后ms/条1,0001.654 s0.013 s0.0134,0008.879 s0.047 s0.01210,00033.541 s0.113 s0.011耗时与条数成正比了同一条代码路径297×提速。写法上没有花招直接write_feature循环最后layer.sync()一次。一个会让你白测一天的坑⚠️性能必须在终端里测不能按 PyCharm 的 Debug 跑。PyCharm 的调试器pydevd给每个 Python 帧装 line tracer每条字节码回调一次。本库是纯 Python全额上税GDAL 是 C 扩展tracer 进不去税率为 0。在调试器里比较纯 Python 实现 vs C 实现量到的是 tracer 的税率差不是算法差。六、可选加速器装上就快没装照样全功能三个可选依赖性质不一样这个区分很重要没有它影响numpy自动回退纯 Python大环上面积/周长/质心慢×21~25C 扩展_gdbaccel.c自动回退纯 Python整层几何解码慢 ~9×pyproj没有就没有这个功能to_crs()报 ImportError前两个是有就快些dependencies里仍然是零pyproj 是没有这个功能所以单独占一个 extra[crs]。两个加速器都有差分闸门守着拿两个实现互相对拍逐位相同才算过python tools/verify_accel.py# C vs 纯 Python:8 个库 11,885 条 / 6,335,486 顶点python tools/verify_numpy.py# numpy 快路径 vs 顺序路径 Fraction 精确解校准⚠️ 一个具体的坑.pyd绑解释器版本cp311/cp313各一份。只编了一份却用另一个解释器跑会HAS_ACCEL False静默回退——实测因此从 4.4 s 变成 56.5 s。哪个解释器跑代码就用哪个编。七、我诚实列一下什么情况不该用它一个库的价值一半在于它拒绝做什么。以下都是明说的边界不是 bug能替换 GDAL 的部分读写.gdbversion 3 / ArcGIS 10.xversion 4 可读建库、建图层、增删改要素产出能被 ArcGIS/QGIS 直接打开度量、构造、DE-9IM 谓词、buffer、overlay 四算子WKT / GeoJSON 出口__geo_interface__可直接喂 geopandas/folium不能替换 GDAL 的部分请继续用 GDAL任何非 FileGDB格式。这个库只认.gdb没有 Shapefile、GeoPackage、GeoJSON 驱动也没有栅格。需要精确拓扑的场景。本库用浮点方向判定无 snap-rounding、无精确算术。已知的一处1e-9 高的薄片平移到 UTM 量级后只有 2.1 ULP 高双精度网格上根本没有内部点可采样。GEOS 靠组合式拓扑图躲过采样式实现躲不过。is_valid()/is_simple()是部分实现只查环闭合/顶点数下限/单环自交/洞在壳内。is_valid() True不等于OGC 有效。distance()是O(n·m)无线段包围盒预筛。写回不是逐位幂等全语料 21,217 条读→写→读152 条0.7%字节不同环顺序/绕向规范化 ≤1 个量化步长的舍入。面积全部一致是表示差异不是几何差异——但要做哈希去重的场景不能拿本库的输出当原文。buffer/ overlay 结果恒为 2DZ/M 丢弃环是Esri 绕向GEOMETRYCOLLECTION是纯内存类型、写不进 .gdb。八、验证情况383 个用例四配置全绿默认 / 关 numpy / 关 C 扩展 / 换解释器——改一个全局名字最容易漏掉某个引用点所以每条都过。库本体19,143 行18,548 行.py 595 行.c测试6,192 行。差分闸门拿两个实现互相对拍逐位相同才算过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村行政区划import shapely/import osgeo只出现在tools/里不进包、不进tests/、不进pyproject.toml。没装就 skip 并exit 0。九、现在就想试pipinstall.# 纯 Python零运行时依赖pipinstall.[all]# 附带 numpy pyprojfrompyopenfilegdbimportOpenFileGDB,GdbField,Geometry,FGFT_STRING,FGFT_DOUBLE# 读withOpenFileGDB.open(D:/data/x.gdb)asgdb:layergdb.get_layer(村行政区划)featslayer.read_features(where县名 朝阳区,limit100)# 写gdbOpenFileGDB.create(D:/data/new.gdb)layergdb.create_layer(监测点,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()# ← 关闭 结账必须MIT 许可。最后这个项目让我改了一个习惯在说按 X 实现之前先确认 X 真的实现了它。GDAL 在格式上是一座金矿——文档没有代码就是规范逐字节照抄就能得到一个能被 ArcGIS 打开的库。但在空间计算上它是一张贴在 GEOS 上的标签纸。把这两件事分清楚才知道什么时候可以照抄、什么时候必须换参照系、什么时候必须自己对着规范从零写并拿真参照实现对拍。以及那个我最想记住的真值表、不变量、假 oracle 三样都拦不住真正的 bug。只有跑真参照实现才抓得到。仓库https://github.com/lizhiziwang/pyopengdb · 详见README.md· 可选 C 扩展编译见ACCEL.md