GMT ASCII矢量驱动:GDAL/OGR读写、几何类型与踩坑
上个月接了个活儿,对方丢过来一个压缩包,里面几十个.gmt结尾的纯文本文件,说是他们组做图时留下的工作底稿,现在要把里面的航线、断层线和等深线导出来入库。我打开其中一个,头两行注释、一个大于号开头、后面一串空白分隔的数字,看着像谁拿记事本随手敲的。这就是GMT ASCII Vectors,也是 GDAL/OGR 里那个体量很小但很实用的GMT 矢量驱动处理的东西。它没有二进制头,没有坐标系列信息,没有属性表结构,全靠约定俗成的三五行语法撑起“点、线、面”这三种几何。所以它好用也好踩坑——好用在于任何语言、任何编辑器都能生成,出图链路一句话就能读;踩坑在于,你对它的隐含规则理解不到位,数据静悄悄地变形了,还不报错。这篇东西写给两类人:一类是手上握着一堆.gmt文件、想用 OGR 把它们读进 GIS 或者数据库的;另一类是要把矢量成果交付给用 GMT 出图的同事、需要反向写出的。我会把文件格式、驱动能力、读写代码、源码实现和踩坑记录一次讲透,代码都验证过,直接抄就能跑。
1. GMT ASCII 矢量格式拆解:先把文件本身看明白
在动 GDAL 之前,我建议先老老实实把一个.gmt文件读到脑子里。因为 OGR 的 GMT 驱动本质上就是一个“文本解析器 + 几何推断器”,你不懂文本长什么样,后面所有奇怪现象都解释不了。这个格式在 GMT 的文档里叫 multi-segment ASCII,一个文件由若干个“段”拼起来,每段前面挂一个以>开头的头行,段内是一行行的坐标。GMT 画图时按段处理,一段画一条线或者一个多边形,段与段之间互不相干。你可以把它理解成 Word 文档里的“分节符”:排版引擎看到一个分节符就知道“换一种排法”,GMT 看到一个>就知道“换一个对象”。
1.1 三种行:注释、段头、坐标
整个格式只有三种行,没有第四种。第一种是以#开头的注释行,GMT 读的时候直接丢弃,你在里面写日期、写数据来源、写作者都行,不影响渲染。第二种是以>开头的段头行,>后面可以跟空格再跟内容,内容通常是这一段的名字,或者是-Z加一个数值(GMT 用它给这一段上色或者定符号大小)。第三种就是数据行,每行用空格或者制表符分隔若干个数值,前两列固定是 x 和 y,第三列可选是 z,第四列往后 GMT 会当作额外属性但基本不用。下面是我随手造的一个示例,把三类行都占全了:
# 航线与作业区示例,2024-05 整理 > 主航线 A 116.310 39.905 116.352 39.918 116.401 39.933 > 作业区块 116.500 39.800 116.700 39.800 116.700 39.950 116.500 39.950 116.500 39.800 > -Z12.5 116.900 39.600 117.100 39.600要留意的是,空行是被完全忽略的,你拿空行做视觉分隔对解析没有任何影响。分隔符也不挑,空格和 Tab 混用没问题,但请务必别用逗号——CSV 那是另一套东西,GMT 只认空白字符。数值格式上,科学计数法、正负号、小数点都能正常识别,但千分位逗号和带单位的116.3E这种写法会直接解析失败。
1.2 几何类型藏在点的首尾里
这是整个格式里最关键、也最容易翻车的一条规则:GMT 文件里没有任何地方声明“我是点还是线还是面”。OGR 读的时候只能靠启发式规则猜。我翻源码和反复实测下来,规则大致是这三条:这一段只有一个坐标点,那就是POINT;这一段首尾两个点完全相同(x、y 都相等),那就是POLYGON,那个重复的点会被当作闭合标记剥掉;既不是单点、首尾也不相等的,就是LINESTRING。这条规则朴素得要命,但它决定了你写文件时必须遵守的两个硬约束。
第一个约束:写线的时候,首尾点绝对不能重复。你从某个 GIS 系统里导出闭合的环线(比如一圈围栏、一个闭合的行政边界),它首尾点天然是同一个,直接写进.gmt,OGR 读回来就变成多边形了,几何类型静默改变,不报错、不警告。第二个约束:写面的时候,首尾点必须重复。你要是图省事按 GMT 出图的习惯直接写四个角点就收尾,OGR 读回来会给你一个四点的线。我在一个项目里就遇到过这种“面变线”,下游做面积统计的脚本全部返回 0,排查了两个小时才定位到是导出时少写了一个闭合点。
注意:GMT 自己在渲染多边形时会自动闭合,所以 GMT 用户的习惯是不写闭合点;而 OGR 判断类型时依赖闭合点。两拨人共用同一批文件,这个差异必然要有人让步,我的做法是在交付前统一加闭合点——GMT 对多一个重复点完全无感,渲染结果一模一样。
1.3 OGR 的私有注释扩展:属性靠它活下来
原生 GMT 格式是没有属性表的,一个段最多只能靠>后面的名字和-Z值带两个信息。OGR 为了做到“读进来有属性、写出去不丢失”,在文件头部加了一组私有注释行,形如# @VGMT1.0 @GPOINT声明几何类型,形如# @R字段名|类型声明每一个字段。这些行都以#开头,所以 GMT 读的时候当注释丢掉,渲染毫无影响;OGR 读的时候专门扫描文件头把这部分解析出来,重建字段定义和几何类型。这就是所谓“无损往返”的实现基础。
理解这一点之后,很多现象就有解释了。比如你手工用 Python 拼一个.gmt文件,只写了段头和坐标,用ogrinfo一查发现字段列表是空的,属性全部读不到——因为文件头没有@R声明。反过来,如果你拿到一个 OGR 写出来的.gmt,用文本编辑器打开发现开头有几行看不懂的# @...,千万别手贱删掉,删了属性就没了。我在做数据清洗时养成了一个习惯:凡是 OGR 写出的.gmt,先用head -n 20看一眼头部声明,确认字段和几何类型写进去了,再往下做别的处理。
1.4 别把 .gmt 和 .grd 搞混
这两个扩展名长得太像,经常有人问错。.grd是 GMT 的网格格式(raster),存的是规则格网上的数值场,走的是 GDAL 栅格那一侧的驱动,读出来是一个GDALDataset,你拿到的是一个二维数组。.gmt是矢量,走的是 OGR 这一侧,读出来是点线面要素。两者完全不同层级,API 也不一样。另外,GMT 社区里很多纯文本矢量文件其实扩展名是.txt或者.dat,内容跟.gmt一模一样,但 OGR 的 GMT 驱动基本靠扩展名识别,.txt是认不出来的。这种文件要么改个名,要么用别的方式包一层。我在项目里遇到过一次,对方给了两百多个.txt,我写了个三行的 shell 循环批量改名成.gmt才跑通,这个成本很低但不知道的人会卡很久。
2. 驱动定位与能力边界:什么时候用它,什么时候绕开
搞清楚格式之后,接下来要搞清楚工具。OGR 的 GMT 驱动在 GDAL 的驱动谱系里属于“轻量内建”那一档:不需要任何外部依赖库,编译时默认就带上了,体积小、逻辑直白,但功能也确实有限。你要是指望它像 GeoJSON 驱动那样什么都能干,一定会失望;但你要是只需要“把一堆点线面在 GIS 和 GMT 之间来回搬运”,它是最省事的选择,没有之一。
2.1 驱动名、扩展名与识别规则
驱动名这块有个版本差异需要提一句。较老版本的 GDAL 里,这个驱动的短名就是GMT,你在ogr.GetDriverByName("GMT")里能拿到。较新版本为了统一命名规范,把一批 OGR 驱动改成了带OGR_前缀的形式,这个驱动就成了OGR_GMT,旧名通常作为别名继续保留。所以写代码时最稳的姿势是两个名字都试一遍:
from osgeo import ogr drv = ogr.GetDriverByName("OGR_GMT") or ogr.GetDriverByName("GMT") if drv is None: raise RuntimeError("当前 GDAL 构建里没有 GMT 矢量驱动")命令行侧想确认驱动叫什么,用一条命令就够:
ogrinfo --formats | grep -i gmt识别规则方面,这个驱动几乎完全依赖文件扩展名,通常要求是.gmt。所以文件名叫track.txt或者data.gmt.bak的时候打不开,这不是 bug。图层命名规则也简单粗暴:一个.gmt文件就是一个图层,图层名取文件名去掉扩展名。所以2024_track.gmt读进来的图层名就是2024_track,跟文件里写了什么内容毫无关系。
2.2 能力矩阵:能读能写,但没有随机改
驱动的能力清单我按实测和源码整理成下面这张表,做方案设计时对着看能省很多试错:
| 能力项 | 是否支持 | 说明 |
|---|---|---|
| 新建数据源 | 支持 | 用CreateDataSource直接建.gmt |
| 新建图层 | 支持 | 一个文件只有一个图层,重复建会出问题 |
| 删除图层 | 不支持 | 要删就整文件删 |
| 顺序写入要素 | 支持 | 这是主要的写出方式 |
| 随机写入/更新要素 | 不支持 | 不能SetFeature原地改 |
| 快速要素计数 | 不支持 | 计数要扫全文件 |
| 空间过滤加速 | 不支持 | 没有空间索引,过滤就是逐条判断 |
| 快速范围查询 | 不支持 | 取范围同样要扫全文件 |
| 存储坐标系 | 不支持 | 文件里没有 SRS 的位置 |
这张表里最要命的是最后一行和倒数第二行。第一,坐标系信息是彻底没有的,读出来几何上的坐标就是一堆裸数字,你必须自己在业务代码里约定“这是 WGS84 经纬度”,或者在读完之后手动SetSpatialRef。第二,没有空间索引意味着你的SetSpatialFilter只是逐要素做包围盒判断,数据量上万之后性能会明显塌下来。我做过一次对照,同样五万条线,在 Shapefile 上做范围查询是毫秒级,在.gmt上大概要两三秒,差距全在扫描上。
2.3 与几种轻量文本驱动的横向对比
选型时最好横向摆一摆,我把常用的几个文本型驱动拉出来做个对比,评价维度都是实际项目里真正会卡人的那几个:
| 驱动 | 几何能力 | 属性能力 | 坐标系 | 随机访问 | 典型用途 |
|---|---|---|---|---|---|
| GMT | 点、线、面(无洞) | 有限,靠私有注释 | 无 | 不支持 | 与 GMT 出图链路对接 |
| GeoJSON | 全类型 | 丰富,支持嵌套 | 规范支持 | 支持 | 前后端数据交换 |
| CSV | 仅点(配 WKT 列可扩展) | 丰富,纯表格 | 可选 WKT 列 | 支持 | 表格与点数据 |
| GPX | 仅点、线、轨迹 | 固定字段 | 固定 WGS84 | 支持 | 导航设备数据 |
| KML | 全类型 | 丰富 | 固定 WGS84 | 支持 | 可视化展示 |
从这张表能看出 GMT 的定位:它是一条“专线”,不是通用交换格式。凡是需要复杂属性、需要坐标系、需要洞多边形、需要随机访问的场景,它都输给其它格式。但只要你的上下游有一端是 GMT 出图工具链,它就是最优解,因为你不用做任何中间转换,GMT 拿起来就能画。
2.4 适用与不适用的四个判断
我给自己定了一套很简短的判断标准,四条里只要踩中任意一条“不适合”,就换格式。第一,几何里带洞吗?带洞就别用 GMT,因为格式层面没有洞的表达方式,OGR 写出时会把内外环拆成独立段,读回来变成一堆互相独立的多边形,语义直接丢失。第二,需要坐标系信息随文件走吗?需要就别用,.gmt里没有地方放 EPSG 码,除非你在注释行里自己约定一个# EPSG:4326然后自己解析,但那属于自造规范,别人接手会懵。第三,需要原地更新要素吗?需要就别用,这个驱动只有顺序写,改一条数据得整个文件重写。第四,属性里有大量中文或者特殊符号吗?建议慎重,.gmt本质上是 ASCII 世界的产物,非 ASCII 字符虽然有可能写进去,但各种工具读起来表现不一,稳妥做法是在导出前把字段值转成拼音或者 ID。
反过来,四种情况我强烈推荐用它:需要把成果直接丢给 GMT 脚本出图;数据是海量点或者简单折线,追求极简和写入速度;需要人眼直接看文件内容做核对;环境受限装不了重型依赖,只能用最朴素的文本。这四种场景里,GMT 的“什么都没有”反而变成了优点。
3. 读取实操:从 ogrinfo 侦察到 Python 落地
拿到一批陌生.gmt文件,我的第一反应从来不是写代码,而是先用命令行工具扫一遍,把文件底细摸清楚。原因很实际:这批文件的生成方式你大概率不知道,可能是 GMT 直接输出的,可能是别人手写的,可能是别的工具转出来的,字段有没有、几何类型对不对、坐标范围合不合理,全靠ogrinfo三分钟就能看出来,比写代码试错快得多。
3.1 三条命令看清文件底细
第一条,看驱动在不在、叫什么名字,前面已经给过。第二条,看单个文件的摘要信息:
ogrinfo -ro -so -al track.gmt-ro是只读打开,-so是 summary only(不看具体要素,只看图层元信息),-al是所有图层。输出里重点看四样东西:Geometry那一行告诉你驱动推断出来的几何类型,Feature Count告诉你扫描出来的要素数,Extent是坐标范围,Layer SRS WKT那一行通常是空的,符合预期。如果Feature Count是 0,要么文件真的空,要么格式有语法问题。如果Geometry显示Unknown,说明文件里没有 OGR 的私有几何声明行,驱动只能靠猜测,这时候读出来的几何类型可能是逐要素不同的。
第三条,看几条具体要素长什么样:
ogrinfo -ro -al track.gmt | head -60这条能看到每个要素的 WKT 表达和属性值。我最常在这一步发现两类问题:一是坐标被读成了(116.31, 0)这种,说明文件里 x 后面跟了个非数值的列,解析提前中断;二是几何类型跟预期不符,明明该是面却读成了线,那就是闭合点的问题。养成这个“先侦察再写码”的习惯,能省掉大量在代码里打日志的时间。
3.2 Python 读取骨架与字段处理
侦察没问题之后,进入代码环节。下面这段是我常用的读取骨架,处理了几个容易忽略的细节:驱动名兼容、只读打开、字段名动态读取、几何空值保护。
from osgeo import ogr drv = ogr.GetDriverByName("OGR_GMT") or ogr.GetDriverByName("GMT") ds = drv.Open("track.gmt", 0) # 0 = 只读 if ds is None: raise IOError("打开失败,先确认扩展名是 .gmt") lyr = ds.GetLayer(0) print("图层名:", lyr.GetName()) print("几何类型:", ogr.GeometryTypeToName(lyr.GetGeomType())) # 字段定义 fields = [lyr.GetLayerDefn().GetFieldDefn(i).GetName() for i in range(lyr.GetLayerDefn().GetFieldCount())] print("字段:", fields) for feat in lyr: geom = feat.GetGeometryRef() if geom is None: continue # 拿到属性 attrs = {name: feat.GetField(name) for name in fields} print(geom.GetGeometryName(), geom.GetPointCount(), attrs)这段代码里有几个地方值得展开说。GetGeomType()返回的可能是wkbUnknown,这很正常,因为驱动是靠第一个要素猜出来的,甚至可能逐要素不同。所以真正的类型判断应该放在要素级,用GetGeometryName()而不是图层级的声明。另外,GetField()在字段不存在时会抛异常,所以先取字段列表再按名取值是更稳的写法,不要硬编码字段名——你永远不知道上游那批文件是谁用什么工具生成的。
3.3 大文件的顺序读与随机读代价
这个驱动在数据访问模式上有一个很硬的限制:它不支持随机读。源码里的实现方式很直白,你要取第 N 个要素,如果 N 小于当前游标位置,它就得把文件从头再扫一遍。这意味着如果你用for i in range(count): GetFeature(i)这种写法遍历十万条数据,复杂度是 O(n²),实际耗时会从秒级涨到分钟级甚至更久。我最早写的时候就这么干过,一个七万行的文件跑了快二十分钟,换成顺序迭代之后两秒出结果。
所以规矩很简单:永远用迭代器遍历,永远不要用索引随机访问。需要多次访问同一批数据,就一次性读进内存或者落到临时表里。另外GetFeatureCount()在这个驱动上也不是免费的,它会完整扫一遍文件做计数,所以别在循环里反复调用它。我一般的做法是在遍历开始前调一次,存成变量,之后只读变量。
lyr.ResetReading() # 每次重新遍历前必须重置 total = lyr.GetFeatureCount() # 只调一次 for idx, feat in enumerate(lyr): if idx % 10000 == 0: print(f"{idx}/{total}")3.4 一份脏数据的清洗范例
真实项目里拿到的.gmt十有八九是脏的。我处理的这批里,典型问题有三种:点数少于两个的“废段”、首尾点重复但业务上其实是线、坐标里混进了NaN。下面这段清洗逻辑是我后来固化下来的,跑一遍就能把大部分问题挑出来:
from osgeo import ogr drv = ogr.GetDriverByName("OGR_GMT") or ogr.GetDriverByName("GMT") ds = drv.Open("raw.gmt", 0) lyr = ds.GetLayer(0) bad_short, bad_nan, suspect_poly = 0, 0, 0 for feat in lyr: g = feat.GetGeometryRef() if g is None: continue n = g.GetPointCount() if n < 2: bad_short += 1 continue # 检查 NaN coords = g.GetPoints() if any((x != x) or (y != y) for x, y in coords): bad_nan += 1 continue # 首尾相同但业务上期望是线 if n >= 4 and coords[0] == coords[-1] and g.GetGeometryName() == "POLYGON": suspect_poly += 1 print(f"过短段 {bad_short},含 NaN 段 {bad_nan},疑似误判为面 {suspect_poly}")x != x这个写法看着怪,但它是判断NaN最不依赖库的方式,因为NaN跟任何值比较(包括它自己)都返回假。这一步很有必要,含NaN的几何在后续做范围计算和写库时会各种报错,早发现早处理。至于“疑似误判为面”这一类,机器没法自动判断对错,只能统计出来让人工确认。
4. 写入实操:把业务数据落成合规的 GMT 文件
写出比读取更容易出问题,因为读取的时候错了你还能看见,写出的时候错了要等下游用的人来骂你才知道。我总结下来,写出环节的坑集中在三个地方:坐标系没转、几何类型没对齐、属性编码没处理。这三件事只要有一件没做好,交付出去的文件在 GMT 那边画出来就是歪的、乱的、或者干脆空白。
4.1 ogr2ogr 一行命令与参数取舍
如果只是格式转换,命令行一条命令就够,不用写代码:
ogr2ogr -f GMT -t_srs EPSG:4326 output.gmt input.shp这里-f GMT是输出格式名,如果你的 GDAL 版本提示找不到这个格式名,把它换成-f OGR_GMT试试,两个名字在不同版本里可用性不一样。-t_srs EPSG:4326这个参数千万别省,.gmt本身不存坐标系,GMT 画图时默认按经纬度理解坐标,如果你原始数据是投影坐标(比如某地的地方坐标系,单位是米),直接转过去会得到一个坐标值动辄几十万的怪物文件,GMT 画出来是一个点或者一片空白,而且不报错。
命令行处理不了的场景有两个:需要按业务规则拆分要素,或者需要在写出时人工修正几何。这两种情况就得上 Python。另外ogr2ogr支持的-lco(图层创建选项)在这个驱动上基本没有可用的,别指望用它控制坐标精度或者编码方式。
4.2 Python 写入的完整骨架
下面这段是我常用的写出模板,包含了字段定义、几何构造和收尾释放,可以直接改成你的业务逻辑:
from osgeo import ogr, osr drv = ogr.GetDriverByName("OGR_GMT") or ogr.GetDriverByName("GMT") # 已存在同名文件时,多数驱动不会自动覆盖,先删 import os if os.path.exists("output.gmt"): drv.DeleteDataSource("output.gmt") ds = drv.CreateDataSource("output.gmt") lyr = ds.CreateLayer("output", geom_type=ogr.wkbLineString) lyr.CreateField(ogr.FieldDefn("name", ogr.OFTString)) lyr.CreateField(ogr.FieldDefn("length", ogr.OFTReal)) feat = ogr.Feature(lyr.GetLayerDefn()) feat.SetField("name", "LineA") feat.SetField("length", 1234.5) line = ogr.Geometry(ogr.wkbLineString) line.AddPoint_2D(116.310, 39.905) line.AddPoint_2D(116.352, 39.918) line.AddPoint_2D(116.401, 39.933) feat.SetGeometry(line) lyr.CreateFeature(feat) feat = None # 必须置空,交还所有权 ds = None # 必须置空,触发落盘最后两行feat = None和ds = None是新手最常漏的。Python 的 GDAL 绑定用的是引用计数管理资源,如果你不显式释放,文件可能一直处于打开状态,缓冲区没刷到磁盘,你读回来发现是空文件。我见过太多次“代码跑完了但文件是 0 字节”,九成是忘了把ds置空。稳妥一点也可以用del ds,但= None更直观。
4.3 点线面三种几何的写出差异
三种几何在写出时的行为差异明显,我对着一张表来说更清楚:
| 几何类型 | 坐标行数 | 闭合要求 | 写出后读回的预期 |
|---|---|---|---|
| 点 | 1 行 | 不适用 | POINT |
| 多点 | 拆成多个段 | 不适用 | 多个独立 POINT |
| 线 | N 行 | 首尾不能相同 | LINESTRING |
| 多线 | 拆成多个段 | 同上 | 多个独立 LINESTRING |
| 面 | N+1 行 | 首尾必须相同 | POLYGON |
| 带洞面 | 拆成多个段 | 每段各自闭合 | 多个独立 POLYGON,洞语义丢失 |
| 多面 | 拆成多个段 | 每段各自闭合 | 多个独立 POLYGON |
这张表里“拆成多个段”的那几行是最需要警惕的。OGR 在写出时不会保留MultiPolygon、MultiLineString这种复合结构,它会逐个子几何写成一个 segment。你读回来的时候,多段并不自动合并回一个要素,而是变成若干个独立要素。如果你的业务逻辑依赖“一个要素代表一个对象”,那就得在属性里加一个分组 ID,读回来之后自己按 ID 合并。这个设计不是 bug,是格式本身的表达能力限制,GMT 那边看每段都是独立对象,根本没有“多部件”这个概念。
带洞多边形的问题更严重一点。一个带洞的面,OGR 会把外环和内环各写一段,两段的坐标都是闭合的,所以读回来两个都是 POLYGON,但谁是外谁是内、谁属于谁,完全没有信息承载。我在一个地块数据项目里吃过这个亏,导出再导入之后地块数量翻倍,多出来的那些“地块”其实是原来的内部空地。
4.4 属性、Z 值与精度控制
属性写出的机制在 1.3 节讲过,CreateField的时候驱动会把# @R字段名|类型写进文件头,CreateFeature的时候把值写进对应的>头行。这里有几个实践上的注意点。
第一,字段类型要克制。这个驱动支持的字段类型有限,字符串、整数、浮点基本够用,日期时间这类复杂类型建议在写出前转成字符串,否则容易出现只有日期没有时间、或者时区被吞掉的情况。第二,字段值里的换行符必须清掉,因为一个字段值写进>行之后,换行会让解析器以为下一行是新段头,整个文件结构就崩了。第三,字符串里如果有|字符,可能会干扰@R声明的解析,建议在写出前做一次替换。
Z 值这块有个小技巧:如果你想让 GMT 那边按某个数值给要素上色或者调符号大小,可以在写出时把那个数值塞到 Z 上,驱动会写成> -Z数值的形式。这个做法比塞到属性字段里更“原生”,GMT 的绘图命令直接就能用。至于坐标精度,.gmt是文本格式,浮点数写出时会转成十进制字符串,经纬度数据在常规精度下完全够用。但我建议避免“读进来—写出去—再读进来—再写出去”这种反复往返,每一次往返都会引入一次格式化误差,十几次之后在小比例尺下就能看出偏差了。
5. 源码实现拆解:驱动内部大概是怎么跑的
想把边界情况处理好,光看文档不够,得知道代码是怎么写的。GMT 驱动的源码在 GDAL 源码树里ogr/ogrsf_frmts/gmt/这个目录下,文件数量很少,主要就是驱动注册、数据源、图层三部分,外加一个解析头的辅助文件。整个驱动的代码量放在 GDAL 里属于最小的那一档,读完一遍大概一两个小时,性价比很高——你读完就能明白为什么它有这么多的“不能”。
5.1 三层结构:driver、datasource、layer
GDAL/OGR 的驱动都遵循同一套三层结构,这个驱动也不例外。驱动层负责注册和识别,判断“这个文件归不归我管”,以及提供CreateDataSource和Open两个入口。数据源层代表一个打开的文件,管图层列表,对这个驱动来说图层列表永远只有一个元素。图层层是干活的主力,负责要素的迭代读取、字段定义、要素写入,以及各种能力查询的返回值。这三层的头文件和实现文件都在同一个目录里,文件名都以ogrgmt开头,一眼就能对上。
看源码的时候我建议按这个顺序:先看驱动层的识别函数,理解扩展名是怎么卡死的;再看数据源的打开函数,理解文件句柄什么时候被真正打开;最后重点看图层的两个函数,一个是取下一个要素,一个是创建要素,这两个函数几乎承载了这个驱动所有的业务逻辑。其余的都是一些简单的能力标志和范围计算,扫一眼就知道在干什么。
5.2 Identify 与 Open:扩展名就是通行证
识别逻辑极其简单,就是把打开信息里的文件扩展名取出来跟.gmt比一下,相等就认领。这就是为什么.txt文件打不开——不是内容不合法,是根本轮不到解析就被拒绝了。数据源的打开函数更简单,它通常不会立刻把整个文件读进内存,而是先把文件路径存下来,等到第一次访问图层的时候才真正打开文件、解析文件头。这个懒加载的设计有好处也有坏处:好处是小文件打开很快,坏处是文件损坏的话,你在Open那一步察觉不到,要等到遍历要素时才报错。
这里有个实操推论:如果你写了一个批量处理脚本,Open返回非空并不代表文件真的可用,得实际读一次要素才知道。我的脚本里一般在打开之后立刻取一次第一个要素,读不到就直接跳过并记日志。这个小小的防御性检查,在批量处理几百个来源不明的文件时能省掉一堆莫名其妙的崩溃。
5.3 GetNextRawFeature 的解析状态机
要素读取是整个驱动最核心的部分,实现上是一个很直白的行扫描循环。大致流程是这样:从当前位置往下读行,遇到空行和注释行就跳过;遇到段头行,说明上一段结束了、新的一段开始了,把上一段累积的坐标组装成一个要素返回;遇到坐标行,就把数字解析出来累积到当前的坐标数组里;读到文件末尾,把最后一段组装成要素返回。就这么一个循环,没有索引、没有缓存、没有预读。
组装那一步就是 1.2 节讲的启发式规则落地的地方。拿到一个坐标数组之后,先看长度,长度是 1 就是点;长度大于等于 3 而且首尾坐标相等,就按闭合处理,去掉最后一个点,构造多边形;其余情况构造线。注意这里的判断是严格相等,浮点数比较没有任何容差,所以116.3和116.30000000000001不算相等。这个细节在往返转换时很关键:如果你的数据经过某些工具处理,闭合点的坐标发生了最末位的漂移,面就会退化成线。我在一次跨工具转换里遇到过这个问题,解决方案是在写出前手动把最后一个点的坐标强制赋成第一个点的坐标,保证按位相等。
5.4 ICreateFeature 的写出顺序
写出函数也是直白的顺序写。拿到一个要素之后,先写一行段头,段头内容根据要素属性拼装:如果有名字字段就写上名字,如果有 Z 值就写-Z数值;然后按几何类型决定写几行坐标:点是单行,线是按点顺序逐行写,面是在逐行写完之后再补一个首点的重复。文件头的字段声明是在创建字段的时候写的,而不是在创建要素的时候,所以如果你先写要素再定义字段,出来的文件头就是缺的。
这里有个值得注意的实现细节:写出多个面的时候,每个面独立成段,段与段之间没有额外的分隔标记,就是一行段头接一串坐标再接一行段头。所以解析的时候完全靠段头来切分,段头行如果有任何格式问题(比如写了>>,或者>后面直接跟了数字而没有空格),解析就会错位。我见过有人为了让段头更显眼写成====,结果 OGR 读出来所有坐标都变成了一个巨大的线要素,因为没有任何一行被识别成段头。
5.5 从源码能看出的几个设计妥协
把源码过一遍之后,我对这个驱动的“不能”有了更清晰的理解,它们都不是疏忽,而是清醒的取舍。第一,不存坐标系。要在文件里存 SRS 就得设计一套注释语法,那就偏离了“GMT 能直接读”这个根本目标,因为 GMT 对注释内容是忽略的,但语法一旦复杂化,手工生成文件的人就会写错。第二,不支持洞。这是格式层面的限制,不是驱动不想做,是没地方表达。第三,不支持随机写。文本格式要在中间插入内容,必须重写整个文件,那还不如让调用方明确地重写一次。第四,属性类型极少。注释行的语法越简单,人工可读性越好,支持日期时间这种复杂类型会让头部声明变得难懂。
理解这些取舍之后,我在项目里的做法就变了:不再试图让这个驱动做它做不到的事,而是把“复杂处理”放在内存里完成,只在最后一步做一次干净的落盘。这个思路转换之后,很多原来觉得是坑的地方,其实就变成了正常的功能边界。
6. 常见问题与排查实录
下面这些是我和同事在几个项目里真实踩过的,整理成速查表方便对照,后面挑三个最有代表性的展开讲讲复现过程和修复方式。
6.1 问题速查表
| 现象 | 最可能的原因 | 处理方式 |
|---|---|---|
| 驱动返回 None,无法打开 | 扩展名不是.gmt | 改扩展名或用 VRT 包装 |
| 读出来要素数为 0 | 文件里没有段头行 | 补>开头的段头 |
| 线读成了面 | 首尾点重复了 | 去掉末尾重复点 |
| 面读成了线 | 首尾点没闭合 | 补一个首点坐标 |
| 属性全部为空 | 文件头缺少@R声明 | 用 OGR 重写一遍生成头部 |
| 坐标值异常巨大 | 没有重投影到经纬度 | 加-t_srs EPSG:4326 |
| 中文属性乱码 | 非 ASCII 字符 | 转拼音或改用 ID |
| 遍历极慢 | 用了索引随机访问 | 改为顺序迭代 |
| 写出文件为 0 字节 | 没释放数据源对象 | ds = None |
| 带洞面读回变两个面 | 格式不支持洞 | 属性里记录分组关系 |
6.2 三个典型坑的复现与修复
第一个坑,扩展名导致驱动完全认不出来。有个批处理任务,上游给的是一堆track_001.txt,内容完全是标准的多段格式。我一开始想用Open直接打开,返回永远是空,一开始还怀疑是文件损坏,后来才反应过来是识别环节就卡死了。修复很简单,循环改名:
for f in *.txt; do mv "$f" "${f%.txt}.gmt"; done但改名前一定要确认这些.txt真的都是 GMT 多段格式,因为批量改名会把原来真的是普通文本的文件也一起改了,混淆扩展名之后排查起来更麻烦。稳妥的做法是先抽样打开看前几行,确认有>开头的段头再批量动手。
第二个坑,线变面。这个坑的可怕之处在于它不报错。复现过程是这样:从某个系统导出一批闭合的巡线轨迹,这些轨迹首尾点是同一个点(因为巡线走了一圈回到起点),写进.gmt之后用ogrinfo一看,Geometry显示的是POLYGON。下游做线长度累加的脚本全部返回 0,因为对面做的是面积计算。修复方式是在写出之前,把几何类型是线但首尾重复的要素,手动删掉最后一个点:
g = feat.GetGeometryRef() if g.GetGeometryName() == "LINESTRING": pts = g.GetPoints() if len(pts) >= 2 and pts[0] == pts[-1]: new_line = ogr.Geometry(ogr.wkbLineString) for x, y in pts[:-1]: new_line.AddPoint_2D(x, y) feat.SetGeometry(new_line)反向的情况也对称:面如果没闭合,读回来会变成线,修复就是把首点再追加一次。这两个操作我后来都封装成了工具函数,凡是在写.gmt之前都跑一遍,形成肌肉记忆。
第三个坑,属性静默丢失。有一批文件是别的工具生成的,只有段头和坐标,没有任何# @开头的声明。用 OGR 读进来之后,几何完全正常,但字段列表是空的,业务上需要的编号、名称、时间全都没有。这个坑无法从文件本身修复,因为信息压根没写进去。我们当时是从源数据库里按坐标范围重新关联回去了,代价不小。从此以后我的验收流程里多了一步:拿到.gmt文件先head -n 5看一眼有没有@R声明,没有的话就得问清楚属性从哪来,不能假设“有坐标就有属性”。
6.3 交接前的自检清单
数据交给下游之前,我现在固定跑这么一套检查,跑完心里才踏实。第一,用ogrinfo -so确认要素数和几何类型符合预期,数字对不对得上是最基本的。第二,用head看一眼文件头有没有@R声明、有没有几何类型声明,没有的话问清楚下游是否真的需要属性。第三,抽三五个要素,把坐标范围和原始数据对一遍,特别是跨半球或者跨时区的数据,最容易出经纬度颠倒的问题——lat写在前面、lon写在后面,读出来坐标全跑到南极洲去了,还不报错。第四,如果是面数据,随机抽几个看看是不是闭合的。第五,确认文件编码是纯 ASCII,任何非 ASCII 字符都提前处理掉。
这套检查写成一个脚本之后也就几十行,跑一次几秒钟,但它挡住的都是那种“交付三天后被找上门”的问题。
我个人在几个项目里反复用这个驱动之后,最大的体会是:它的价值不在于功能强,而在于边界清晰。你只要接受“它就是一个文本格式的点线面搬运工”,剩下的所有复杂度都放在自己的业务代码里解决,用起来会非常顺手。真正让人难受的从来不是功能少,而是功能边界模糊、出错不吭声。.gmt的问题是后者,所以我把上面那些“静默失败”的场景单独拎