Java读取SHP文件并提取点线坐标:GeoTools与手写解析实践
1. shp 不是单个文件:先把三件套的关系捋清楚
很多人第一次接到"把 shp 文件里的点和线坐标提出来"这个需求时,直觉反应是打开文件读文本,结果用记事本一开,满屏乱码。这不是编码问题,而是 shp 本身就是二进制格式,它的坐标是按 IEEE 754 双精度浮点数直接写进字节流里的,用文本编辑器看当然是一片混乱。我在实际项目里见过不止一次有人试图用正则去匹配坐标,最后耗掉两天才回头老老实实写解析逻辑。所以动手写 Java 代码之前,先把 shp 这组文件的构成搞清楚,后面 80% 的坑都能提前避开。
shp 文件本质上是 ESRI 定义的一种矢量几何存储格式,它只存几何,不存属性,也不存坐标系定义。一个能被 GIS 软件正常打开的 shp 数据集,通常是一组同名不同后缀的文件凑在一起,缺一个就可能出问题。你从别人那里拿到的压缩包里如果只有 .shp,那这份数据大概率是不完整的,读取时要么报错,要么拿到的坐标没有参照系,用起来会非常别扭。
1.1 .shp/.shx/.dbf/.prj 各自管什么
这四个文件的职责划分很清晰,理解之后你就能判断"我只想要坐标,到底需要读几个文件"。
| 后缀 | 作用 | 只取坐标时是否必需 |
|---|---|---|
| .shp | 存几何图形本体,点线面的坐标全在这里 | 必需 |
| .shx | 几何记录索引,记录每条记录在 .shp 中的字节偏移 | 可选,顺序遍历时用不上 |
| .dbf | 存属性表,dBASE 格式,一条记录对应一个几何 | 只在需要属性时读 |
| .prj | 存坐标参考系的 WKT 文本描述 | 涉及坐标转换时必需 |
| .cpg | 声明 .dbf 的字符编码 | 读中文属性时必需 |
如果你的需求纯粹是"把线的拐点坐标捞出来",那么只读 .shp 就完全够了,.shx 的存在意义是让你能按记录号随机定位,做顺序全量遍历时它反而是多余的。这一点很关键,因为很多人用 GeoTools 打开数据源之后,会不自觉地把属性也一起读了,几百万条记录的属性表加载进来,内存直接顶到几个 G,其实完全没必要。
.shx 的结构其实很简单,100 字节的文件头(和 .shp 的头几乎一样)加上每 8 字节一条的索引项:前 4 字节是偏移量(以 16 位字为单位,大端),后 4 字节是内容长度。真要做随机读取的时候,自己解析 .shx 也就二十行代码的事。
1.2 ShapeType 决定了字节怎么排
shp 文件头的第 33 到 36 字节(从 0 开始数是索引 32)存了一个整型的 ShapeType,这个值决定了后面每一条记录体的字节布局。常见的有这么几个:
- 0:Null Shape,空几何,占位用
- 1:Point,单个点,两个 double
- 3:PolyLine,折线,可以多段
- 5:Polygon,面
- 8:MultiPoint,多点
- 11 / 13 / 15:带 Z 值(高程)的 Point / PolyLine / Polygon
- 21 / 23 / 25:带 M 值(度量值)的对应类型
这里有一个非常容易踩的坑:同一个 shp 文件里所有记录的 ShapeType 必须一致,但每条记录体内部还会再写一遍 ShapeType。有些工具生成的文件,记录体的 ShapeType 和文件头不一致,解析时以记录体为准更安全。我就遇到过一次,文件头写 3(PolyLine),实际记录体里混了几条 0(Null Shape),如果不做判断直接按 PolyLine 解析,会直接抛数组越界。
另外,Z 类型和 M 类型的字段是追加在标准坐标之后的,如果你只关心平面 XY,读的时候可以把多余部分跳过,但偏移量必须算准,否则后面所有记录全错位。
1.3 文件头里两个长度字段的双字节序陷阱
shp 格式最反人类的设计之一,是它在同一个文件里混用了两种字节序。文件头里的文件长度字段(第 24 到 27 字节)是大端存储,单位是 16 位字;而版本号、ShapeType、包围盒的四个 double,全都是小端存储。记录头里的记录号和内容长度是大端,记录体里的 ShapeType 和坐标又是小端。
这意味着你在 Java 里用DataInputStream读的时候,readInt()默认是大端,读文件头那个长度字段正好合适,但读版本号就必须Integer.reverseBytes()转一下。我见过有人为了省事全程用ByteBuffer然后忘了order()切换,结果包围盒读出来是1.0E-322这种离谱的值,还以为是文件损坏。
还有一个隐蔽的问题:文件头里的"文件长度"字段,单位是 16 位字(2 字节),所以真实字节数要乘以 2。有些工具在写出时把这个字段写错,或者文件被截断过,导致这个值和实际文件大小对不上。稳妥的做法是以文件实际长度为准做循环终止判断,而不是依赖头部声明。
提示:判断一个文件是不是合法 shp,最快的办法是读前 4 字节,转成大端整数,应该等于 9994。这个魔数不对,后面怎么调都是白费力气。
2. Java 读 shp 的三条路线与选型判断
搞清楚格式之后,就要决定用什么工具去读。这个决策没那么简单,因为不同项目的约束条件差别很大:有的项目只要几十个点、跑一次就完事;有的要每天定时处理上千个文件、单文件几百万条记录;还有的项目运行在依赖管控极严的环境里,加一个新 jar 都要走审批。我在不同场景下用过三种方案,各有各的活法,下面把判断依据摊开讲。
先说结论性的判断逻辑:如果你需要坐标转换、需要读属性、需要处理面和多部件几何,闭眼选 GeoTools;如果你只要 XY 坐标、文件量大、想控制内存和启动时间,手写解析更省心;如果你的数据源混杂着 shp、GeoJSON、PostGIS 甚至 CAD 格式,那应该考虑更上层的抽象而不是纠结单个格式。
2.1 GeoTools:一把梭,但依赖树不小
GeoTools 是 Java 生态里做 GIS 绕不开的库,gt-shapefile模块把 shp 的读写、坐标系识别、几何模型(基于 JTS)全都封装好了。用它读一个 shp,核心代码不超过十五行,坐标转换一行CRS.findMathTransform加一行JTS.transform就完事。对于需要处理 .prj 定义的投影、需要把投影坐标转成经纬度的场景,这是最省事的路。
代价是依赖体积。gt-shapefile会带进来gt-main、gt-referencing、gt-metadata、gt-opengis一大串,还有 JTS 和一堆坐标系定义资源文件。整个依赖树拉下来,几十兆是常态。如果你的应用是个轻量的命令行工具或者要打进 Android 包,这个体积就有点难受了。
第二个代价是版本对齐。GeoTools 的版本和 JDK 版本绑定得比较紧,新版本普遍要求 JDK 11 以上,一些老项目还卡在 JDK 8,就只能用 2.x 系列的旧版本。而旧版本和新版本的 API 有差异,网上搜到的示例代码经常对不上,复制过来编译不过是很常见的事。我建议先把版本号在pom.xml里锁死,然后所有示例都按这个版本来,别东拼西凑。
第三个要注意的是几何对象的生命周期。GeoTools 的SimpleFeature在迭代器下一次next()之后,内部缓冲区可能被复用,如果你把 Feature 直接塞进 List 里晚点再处理,拿到的可能是错乱的数据。正确做法是在循环体内部就把需要的坐标和属性提取成自己的 POJO,不要持有 Feature 引用。
2.2 纯手写二进制解析:只取坐标时最轻
如果需求就是"给一个 shp,返回所有点和线的坐标数组",手写解析其实非常划算。整个逻辑就是读 100 字节文件头 + 循环读记录,核心代码一百多行,零第三方依赖,启动时间几乎为零,内存占用可控,还能顺便做成流式处理——读一条吐一条,几百万条记录也不会 OOM。
手写方案的优势在大文件上体现得特别明显。GeoTools 的getFeatures()默认行为会做一些缓存和索引准备,虽然有SimpleFeatureIterator做流式,但底层还是会预读一部分。而自己写的话,用BufferedInputStream包一层,每次只读一条记录的内容长度,处理完就丢,内存曲线是平的。
风险也很清楚:一是几何类型要自己覆盖全,Polygon 的多环结构、MultiPoint、Z/M 值,都得单独处理;二是容错要自己做,遇到截断文件、非法 ShapeType 得有降级策略;三是坐标系只能靠自己读 .prj 再找库转换,这部分反而又绕回 GeoTools 或者 Proj4j 了。所以我的经验是:只要涉及坐标转换,手写方案就不划算了,因为投影转换这块的坑远比读二进制深。
2.3 其他方案与适用边界
还有一个思路是通过 GDAL 的 Java 绑定。它的格式支持面最广,shp 只是其中一种,遇到编码混乱、文件损坏的疑难文件时,GDAL 的容错能力比纯 Java 实现强不少。但它的部署麻烦,需要本地库文件跟着走,跨平台打包要分别准备 Windows、Linux、macOS 的二进制,CI 环境里配置起来挺折腾。如果服务要跑在容器里,还得注意基础镜像里有没有对应的依赖。
另外有些人会考虑先把 shp 转成 GeoJSON 再用 JSON 库读。这个路子在小数据量、一次性处理时确实省事,但转换本身需要工具,而且 GeoJSON 体积比 shp 大好几倍,几百万条记录转出来可能就是几个 G 的文本,解析时内存直接爆掉。我只在"临时看一眼数据长什么样"的时候用这招,正式流程里不推荐。
选型这块还有个隐性成本要算:团队里其他人能不能接手。GeoTools 的方案虽然重,但网上资料多,出了问题好查;手写解析器一旦原作者离职,后面的人读那一堆reverseBytes会非常痛苦。所以如果代码要长期维护,注释一定要写足,尤其是每个字段的字节偏移和字节序。
3. 用 GeoTools 把点和线的坐标完整取出来
假设你已经决定走 GeoTools 这条路,下面把从依赖配置到坐标提取的完整链路走一遍。这段代码我在几个项目里反复用过,改动不大,可以直接抄。需要注意的是版本号,我这里用的是较新的稳定版,如果你的 JDK 版本低,把版本号往下降,但 API 主体逻辑基本一致。
整个流程分四步:配置依赖 → 打开 DataStore → 拿到 FeatureSource → 遍历并提取几何。每一步都有容易出问题的地方,我按顺序讲。
3.1 依赖与仓库配置
Maven 依赖只需要一个gt-shapefile,它会传递带入其他模块。但要注意,GeoTools 的正式发布版本放在官方的 OSGeo 仓库里,中央仓库不一定有全量版本,所以pom.xml里要额外声明仓库地址。
<repositories> <repository> <id>osgeo</id> <name>OSGeo Release Repository</name> <url>https://repo.osgeo.org/repository/release/</url> <snapshots><enabled>false</enabled></snapshots> <releases><enabled>true</enabled></releases> </repository> </repositories> <dependencies> <dependency> <groupId>org.geotools</groupId> <artifactId>gt-shapefile</artifactId> <version>31.0</version> </dependency> <!-- 需要坐标转换时加这个,通常会被传递依赖带入,显式声明更稳 --> <dependency> <groupId>org.geotools</groupId> <artifactId>gt-referencing</artifactId> <version>31.0</version> </dependency> </dependencies>版本号一定要所有gt-*模块保持一致,混用版本是依赖冲突的头号来源。如果项目里已经引入了其他 GIS 相关库,比如某个版本的 JTS,很容易出现NoSuchMethodError。排查方法是mvn dependency:tree看 JTS 是不是被降级了,GeoTools 对 JTS 版本有明确要求,必要时用dependencyManagement锁死。
如果最终产物是个 fat jar,还要注意META-INF/services下的 SPI 文件在打包时会被覆盖,导致 GeoTools 找不到数据源实现,报No factory found之类的错。用 shade 插件的话要加ServicesResourceTransformer,用 assembly 插件则要合并而不是覆盖。
3.2 打开数据源到遍历 Feature 的主干代码
打开 shp 有两种写法,一种是传 URL,一种是直接传 File 对象。传 URL 更通用,因为 GeoTools 支持从 jar 包内、网络地址读取,但本地文件用file.toURI().toURL()转换一下就行。参数 map 里的 key 在不同版本有细微差别,url是通用的那个。
public static void readShapefile(String shpPath) throws Exception { File file = new File(shpPath); if (!file.exists()) { throw new IllegalArgumentException("文件不存在: " + shpPath); } Map<String, Object> params = new HashMap<>(); params.put("url", file.toURI().toURL()); DataStore dataStore = DataStoreFinder.getDataStore(params); if (dataStore == null) { throw new IllegalStateException("无法识别该数据源,检查是否为合法 shp"); } String typeName = dataStore.getTypeNames()[0]; SimpleFeatureSource source = dataStore.getFeatureSource(typeName); SimpleFeatureCollection collection = source.getFeatures(); try (SimpleFeatureIterator it = collection.features()) { while (it.hasNext()) { SimpleFeature feature = it.next(); Object geomObj = feature.getDefaultGeometry(); if (!(geomObj instanceof Geometry)) { continue; // Null Shape 直接跳过 } Geometry geometry = (Geometry) geomObj; handleGeometry(geometry, feature); } } finally { dataStore.dispose(); // 释放文件句柄,Windows 下不调用会删不掉文件 } }这段代码有两个细节值得强调。第一是dataStore.dispose(),在 Windows 上如果不调用,文件句柄会一直被占用,后续想覆盖或删除 shp 文件会失败,这个坑我踩过好几次,日志里只会显示"文件被另一个进程占用",很难联想到是这里没释放。第二是feature.getDefaultGeometry()可能返回 null,对应 Null Shape 记录,必须判空,否则后面强转直接 NPE。
另外collection.features()返回的迭代器一定要用 try-with-resources 关闭,它内部持有文件游标。有些人习惯先collection.toArray()再遍历,这在十万条以内还行,上百万条就是灾难,内存直接翻几倍。
3.3 Point / LineString / MultiLineString 的取值写法
拿到Geometry之后,不同类型的取值方式不一样。最通用的方式是getCoordinates(),它返回一个Coordinate[]数组,不管什么几何类型都能拿到所有坐标点,但会丢失几何结构信息——你不知道哪些点属于同一条线。
private static void handleGeometry(Geometry g, SimpleFeature feature) { String geomType = g.getGeometryType(); if (g instanceof Point) { Point p = (Point) g; double x = p.getX(); // 经度或投影 X double y = p.getY(); // 纬度或投影 Y System.out.printf("点: x=%.6f, y=%.6f%n", x, y); } else if (g instanceof LineString) { LineString line = (LineString) g; int n = line.getNumPoints(); for (int i = 0; i < n; i++) { Coordinate c = line.getCoordinateN(i); System.out.printf("顶点[%d]: %.6f, %.6f%n", i, c.x, c.y); } } else if (g instanceof MultiLineString) { MultiLineString multi = (MultiLineString) g; for (int i = 0; i < multi.getNumGeometries(); i++) { LineString part = (LineString) multi.getGeometryN(i); // 每一段单独处理,不要直接拍平 for (int j = 0; j < part.getNumPoints(); j++) { Coordinate c = part.getCoordinateN(j); System.out.printf("段[%d] 点[%d]: %.6f, %.6f%n", i, j, c.x, c.y); } } } else if (g instanceof MultiPoint) { MultiPoint mp = (MultiPoint) g; for (int i = 0; i < mp.getNumGeometries(); i++) { Point p = (Point) mp.getGeometryN(i); System.out.printf("多点[%d]: %.6f, %.6f%n", i, p.getX(), p.getY()); } } }这里要提醒一句:不要图省事直接用getCoordinates()然后丢掉结构。线的分段信息在业务里经常有用,比如一条道路的 shp 里,同一条记录的 MultiLineString 可能包含了好几段互不相连的路段,拍平之后你没法再还原。我见过有人把多段线拍平后做路径计算,结果算出"从北京瞬移到广州"的鬼畜结果,就是因为把不连续的两段当成了同一条线。
还有个容易忽略的点:Coordinate有时会带 Z 值,c.z默认是NaN,如果业务需要高程,得判断一下是不是二维数据。另外getCoordinateN返回的是内部引用还是副本,不同版本行为不同,需要保存储存的话用new Coordinate(c.x, c.y)复制一份更保险。
3.4 属性字段与几何一起带出
实际业务里很少只要坐标,通常还要把属性表的字段一起拿出来,比如道路名称、编号、等级。用 GeoTools 的话这部分很自然,feature.getAttribute("字段名")就能取,但字段名大小写敏感,且不同工具写出的 dbf 里字段名可能被截断成 10 个字符。
SimpleFeatureType schema = feature.getFeatureType(); StringBuilder sb = new StringBuilder(); for (AttributeDescriptor ad : schema.getAttributeDescriptors()) { String name = ad.getLocalName(); if ("the_geom".equals(name)) continue; Object val = feature.getAttribute(name); sb.append(name).append("=").append(val).append("; "); }把字段清单先打出来是最省事的排查手段,我一般会在正式处理前先跑一遍source.getSchema(),把字段名和类型全部打印,确认没问题再写业务逻辑。类型上要注意,dbf 里的数值字段读出来可能是Integer、Long或BigDecimal,直接强转会炸,用String.valueOf()兜底最稳。日期字段更麻烦,dbf 的日期格式是yyyyMMdd字符串,GeoTools 有时会转成java.util.Date,有时保持字符串,需要看具体版本。
4. 手写解析器:100 字节文件头加记录循环
有些场景下 GeoTools 那套确实带不动,比如要在一个启动时间要求 200ms 以内的函数计算里读 shp,或者要处理几百 G 的批量文件。这时候手写解析器就体现出价值了。下面这段代码我在几个数据清洗任务里用过,处理过单文件 400 多万条记录,内存占用稳定在几十兆。
核心思路很简单:DataInputStream包BufferedInputStream,先读 100 字节头,然后循环读 8 字节记录头 + N 字节记录体。难点全在字节序和字段偏移上,我把每一段拆开讲。
4.1 文件头逐字段拆解
文件头固定 100 字节,字段排布如下表。注意"字节序"这一列,混着来是常态。
| 偏移(字节) | 长度 | 内容 | 字节序 |
|---|---|---|---|
| 0 | 4 | 文件码,固定 9994 | 大端 |
| 4 | 20 | 保留,全 0 | - |
| 24 | 4 | 文件长度(单位 16 位字) | 大端 |
| 28 | 4 | 版本号,固定 1000 | 小端 |
| 32 | 4 | ShapeType | 小端 |
| 36 | 32 | Xmin, Ymin, Xmax, Ymax | 小端 double |
| 68 | 32 | Zmin, Zmax, Mmin, Mmax | 小端 double |
读的时候要按顺序来,因为DataInputStream是流式的,跳过一个字段就会错位。读到 100 字节正好结束,第一条记录从第 101 字节开始。
我个人习惯是先把这 100 字节读进byte[],然后用ByteBuffer按需切换 order 解析,代码可读性比一路reverseBytes好得多。但循环体里为了性能,还是用DataInputStream直接读更省事。
4.2 记录体里 Point 与 PolyLine 的排布
每条记录的结构是:8 字节记录头(记录号 + 内容长度,都是大端 int,长度单位是 16 位字),然后是记录体。记录体开头 4 字节是这条记录的 ShapeType(小端 int),后面按类型不同排布。
Point(ShapeType = 1)最简单,就两个小端 double:
[4字节 ShapeType][8字节 X][8字节 Y]PolyLine(ShapeType = 3)复杂一些:
[4字节 ShapeType] [8字节 Xmin][8字节 Ymin][8字节 Xmax][8字节 Ymax] // 包围盒 [4字节 NumParts] // 分段数 [4字节 NumPoints] // 总点数 [4字节 * NumParts] // Parts 数组,每段起始点在 Points 中的索引 [8字节 * NumPoints * 2] // 所有点的 X,Y 依次排列关键理解Parts数组的含义:它不是每段的点数,而是每段第一个点在全局点数组里的下标。比如 Parts = [0, 5, 12],NumPoints = 20,意思就是第 0 段是点 0~4,第 1 段是点 5~11,第 2 段是点 12~19。这个设计很省空间,但第一次看容易理解错。
Polygon(ShapeType = 5)的字节布局和 PolyLine 完全一样,区别只在语义:每个环必须闭合(首尾点相同),且外环顺时针、内环逆时针。如果只取坐标,处理方式和 PolyLine 一致。
4.3 完整可运行代码
下面这段可以直接复制运行,覆盖了 Point、PolyLine、MultiPoint 和多段线的分段还原。
import java.io.*; import java.nio.ByteBuffer; import java.nio.ByteOrder; import java.util.ArrayList; import java.util.List; public class ShpParser { private static int readIntBE(DataInputStream in) throws IOException { return Integer.reverseBytes(in.readInt()); } private static int readIntLE(DataInputStream in) throws IOException { return in.readInt(); } private static double readDoubleLE(DataInputStream in) throws IOException { return Double.longBitsToDouble(Long.reverseBytes(in.readLong())); } public static void parse(String path) throws IOException { try (DataInputStream in = new DataInputStream( new BufferedInputStream(new FileInputStream(path), 1 << 16))) { // ---- 文件头 100 字节 ---- int fileCode = readIntBE(in); if (fileCode != 9994) { throw new IOException("不是合法的 shp 文件, fileCode=" + fileCode); } in.skipBytes(20); int fileLengthWords = readIntBE(in); int version = readIntLE(in); int shapeType = readIntLE(in); double xmin = readDoubleLE(in); double ymin = readDoubleLE(in); double xmax = readDoubleLE(in); double ymax = readDoubleLE(in); // 剩下 32 字节 Z/M 范围,直接跳过 in.skipBytes(32); System.out.printf("版本=%d 类型=%d 范围=[%.3f,%.3f,%.3f,%.3f]%n", version, shapeType, xmin, ymin, xmax, ymax); // ---- 记录循环 ---- long total = new File(path).length(); long read = 100; while (read + 8 <= total) { int recordNumber = readIntBE(in); int contentLengthWords = readIntBE(in); int contentBytes = contentLengthWords * 2; byte[] buf = new byte[contentBytes]; in.readFully(buf); read += 8 + contentBytes; ByteBuffer bb = ByteBuffer.wrap(buf).order(ByteOrder.LITTLE_ENDIAN); int recType = bb.getInt(); if (recType == 0) { continue; // Null Shape } else if (recType == 1) { double x = bb.getDouble(); double y = bb.getDouble(); System.out.printf("记录%d 点: %.6f, %.6f%n", recordNumber, x, y); } else if (recType == 3 || recType == 5) { bb.getDouble(); bb.getDouble(); bb.getDouble(); bb.getDouble(); // 包围盒 int numParts = bb.getInt(); int numPoints = bb.getInt(); int[] parts = new int[numParts]; for (int i = 0; i < numParts; i++) { parts[i] = bb.getInt(); } double[] xs = new double[numPoints]; double[] ys = new double[numPoints]; for (int i = 0; i < numPoints; i++) { xs[i] = bb.getDouble(); ys[i] = bb.getDouble(); } // 按 parts 还原分段 for (int p = 0; p < numParts; p++) { int start = parts[p]; int end = (p == numParts - 1) ? numPoints : parts[p + 1]; System.out.printf("记录%d 段%d 点数=%d%n", recordNumber, p, end - start); for (int i = start; i < end; i++) { System.out.printf(" %.6f, %.6f%n", xs[i], ys[i]); } } } else if (recType == 8) { bb.getDouble(); bb.getDouble(); bb.getDouble(); bb.getDouble(); int numPoints = bb.getInt(); for (int i = 0; i < numPoints; i++) { double x = bb.getDouble(); double y = bb.getDouble(); System.out.printf("记录%d 多点: %.6f, %.6f%n", recordNumber, x, y); } } else { System.out.println("未处理的类型: " + recType); } } } } }这段代码有几个设计取舍值得说明。第一,记录体先整体读进byte[]再解,而不是直接从流里连续读 double,好处是万一某条记录解析出问题,不会影响后续记录的位置——流的位置已经按contentLength走完了。第二,循环终止条件用文件实际长度而不是头部声明,避开了头部字段写错的文件。第三,用readFully而不是read,保证读满,避免短读导致的错位。
如果数据量特别大,new byte[contentBytes]稍微有点浪费,可以复用一个足够大的缓冲区,但单条记录长度可能变化,复用时要注意清零或按需截取。我一般不去抠这点性能,几百万条记录的 GC 压力也还好。
4.4 自测与交叉验证方法
手写解析最怕的是"看起来能跑但结果少了一位小数"。验证方法我总结了三层:
第一层,对照文件头里的包围盒。解析完所有坐标后,自己算一遍 X、Y 的最小最大值,和文件头读出来的四个值比对,误差应该在 1e-9 以内。如果对不上,说明有记录漏读或错位。
第二层,记录数与 .shx 对照。.shx 的文件长度减 100 再除以 8,就是记录总数,和你解析出的条数比一下,数量对得上说明循环没提前退出。
第三层,抽样比对。拿几个坐标去任意一个桌面 GIS 工具里点一下,看数值是否一致。这一步能发现坐标系层面的问题(比如数据其实是投影坐标,但你以为读的是经纬度),也能发现坐标顺序问题。
我一般三层都跑一遍,尤其是处理客户给的"来路不明"的 shp 时,前三分钟多做一步校验,后面能省几小时排查。
5. 坐标拿到手却"不对":投影、单位、顺序的三个高频坑
坐标读出来了,格式也对,但数值一看就不对——这种情况大概占了 shp 相关问题的一半以上。典型表现是:X 是 120 点几,Y 是 30 点几,看着像经纬度,但一放进地图就跑到了非洲西边的海里;或者 X 是三四十万的整数,根本不像经纬度。这背后基本都是坐标系的问题,不是读取代码的问题。
很多人第一反应是"解析器写错了",然后开始怀疑字节序、怀疑偏移量,白白浪费半天。我的经验是:先看数值范围再动手改代码,这一步能省掉 90% 的无效排查。
5.1 .prj 缺失或与数据不匹配
shp 里的坐标是"裸坐标",它自己不知道代表什么。同样的数字 (116.4, 39.9),可能是经纬度,也可能是某个局部坐标系的米,全看 .prj 怎么说。所以第一件事是检查目录下有没有同名 .prj 文件。
如果有,用文本编辑器打开,内容是一段 WKT,看PROJCS还是GEOGCS开头。GEOGCS表示地理坐标系,坐标单位是度;PROJCS表示投影坐标系,单位通常是米。看AUTHORITY["EPSG","xxxx"]拿到 EPSG 编码,后面转换会用到。
如果没有 .prj,那就要靠数值范围反推:
| 数值范围 | 大概率是什么 |
|---|---|
| X 在 -180 到 180,Y 在 -90 到 90 | 经纬度,WGS84 或 CGCS2000 |
| X 在 70 到 140,Y 在 0 到 60 | 中国范围内的经纬度 |
| X 在 1e5 到 1e7 量级整数 | 投影坐标,单位米 |
| X 在 1e6 到 2e7 量级 | 可能是高斯克吕格或 Web Mercator |
如果没有 .prj 又需要转换,那就得靠数据提供方确认。这一步不要猜,猜错了整个数据集都会偏。我见过有人把西安 80 坐标系的数据当成 WGS84 去用,位置偏了上百米,做空间关联时全对不上。
还有个坑是 .prj 内容写错或者和数据不匹配。遇到过文件是从别处拷来的,.prj 是上一个项目的坐标系,数据和定义完全对不上。这时候只能靠业务常识判断——比如一个城市范围的数据,如果算出来的包围盒覆盖了大半个中国,那肯定不对。
5.2 经纬度顺序反了的判定方法
这是最阴险的一个坑。有些数据源在导出时把坐标顺序写反了,本来是 (经度, 纬度),结果存成了 (纬度, 经度)。因为大部分地区的经纬度范围有重叠(经度 30-140,纬度 20-50),你很难从单个点判断,但如果数据点都在中国境内,X 应该在 73~135,Y 应该在 3~54,只要出现 X=39.9、Y=116.4 这种组合,基本就是反了。
批量判断的方法很简单:算所有点的 X 和 Y 平均值,如果 X 的平均值明显小于 Y,而 Y 又超过 90,那必然是反了。或者直接看数据范围,如果 Y 最大值超过 90,经纬度数据不可能,直接判定为反。
修正起来一行代码的事,点对调就可以了,但要注意:只在确认数据本身是经纬度且确实写反的情况下才这么做,如果是投影坐标,对调之后数值依旧错,只会让你更迷惑。
提示:JTS 的
Coordinate里x对应第一个分量,在 EPSG:4326 这个坐标系的官方定义里第一个分量是纬度。但几乎所有数据文件、GeoJSON、实际业务里,第一个分量都是经度。这个矛盾要记住,用 EPSG:4326 做转换时 GeoTools 会按官方定义处理,容易出问题,转换时建议显式用CRS.decode("EPSG:4326", true),第二个参数true表示强制经度在前。
5.3 转换代码与精度控制
确认了源坐标系,要转成目标坐标系(比如统一成经纬度),用 GeoTools 的转换链路最短:
CoordinateReferenceSystem src = CRS.decode("EPSG:4547", true); // 源,投影坐标 CoordinateReferenceSystem dst = CRS.decode("EPSG:4326", true); // 目标,经纬度 MathTransform transform = CRS.findMathTransform(src, dst, true); Geometry converted = JTS.transform(geometry, transform);第三个参数true表示允许宽松匹配,源和目标之间没有精确转换路径时会尝试找近似路径。对于国内常用的 CGCS2000 高斯投影,这个参数很关键,因为很多 EPSG 编码之间的严格转换关系定义不全。
精度方面有几件事要注意。第一,转换会引入误差,通常在厘米级,如果业务要求毫米级,需要确认基准面转换参数,这部分 GeoTools 内置的不一定是最新的。第二,JTS.transform每次调用都会创建转换上下文,循环里调几百万次性能很差,正确做法是把MathTransform提到循环外面复用。第三,转换可能会产生无效坐标(比如超出定义域),转换后要做范围校验,剔除异常值。
另外,投影坐标转经纬度再用double存是可以的,但如果你是先转成字符串再解析回来,精度会掉。我遇到过有人把坐标格式化成%.6f存进 csv,后续再算距离时误差累积到几十米,因为小数点后 6 位大约对应 0.1 米,但累积起来就不行了。中间过程一定保持double,最后展示时才格式化。
6. 工程化细节:编码、内存与批量处理
代码能跑通只是第一步,真正上线之后会遇到一堆工程问题:中文属性乱码、文件句柄泄漏、大文件 OOM、批量处理卡死。这些问题在本地测不出来,一到生产环境就集中爆发。这一节把我遇到过的几个典型问题和解法整理出来。
6.1 DBF 中文乱码与 .cpg
.dbf 文件是 dBASE 格式,它本身没有强制编码声明,编码信息藏在文件头第 29 字节的一个语言驱动 ID 里。但很多工具在生成时压根不写这个字段,或者写了也不准,导致读出来的中文字段全是乱码。
判断方法:目录下有没有 .cpg 文件。有的话里面通常就一行,写着UTF-8或GBK或936,按它来就行。没有的话,就得靠试:用 UTF-8 读一遍,如果中文字段里出现大量�或者莫名其妙的方块,基本可以判定实际是 GBK。
用 GeoTools 的话,可以在参数里指定字符集:
params.put("charset", "GBK");如果不生效,可以在DataStore创建后通过dataStore.getFeatureSource(...)之前设置,或者直接换成ShapefileDataStore显式强转:
ShapefileDataStore sds = (ShapefileDataStore) dataStore; sds.setCharset(Charset.forName("GBK"));手写解析 .dbf 的话,就是按字节读字段,然后new String(bytes, charset)。要注意 dbf 的字段是定长的,字符字段右边用空格补齐,取出来之后要trim(),但有些字段本身就以空格结尾,trim()可能会误删,要按业务判断。
还有一个隐蔽的坑:同一个目录下的多个 shp,编码可能不一样。批量处理时不要全局设置一个字符集,按文件逐个判断,先看 .cpg,没有就试 UTF-8,出现替换字符再回退 GBK。这个逻辑我封装成一个方法,跑了上万个文件,基本没出过问题。
6.2 流式遍历与资源释放
资源释放这块,最容易被忽略的是DataStore.dispose()。前面提过一次,这里再强调:Windows 平台上,如果进程没退出,文件句柄会一直被占用,导致后续覆盖、删除、重命名全部失败。批量处理场景下,如果每个文件都漏一个句柄,跑到几百个文件时就会报"打开文件过多"。
标准写法是用 try-finally 或者 try-with-resources(GeoTools 的 DataStore 不是 AutoCloseable,得手动 finally):
DataStore ds = null; try { ds = DataStoreFinder.getDataStore(params); // 处理逻辑 } finally { if (ds != null) { ds.dispose(); } }手写解析的话,FileInputStream和BufferedInputStream都要关,用 try-with-resources 一次性搞定。另外要注意,DataInputStream包裹之后关闭外层会自动关内层,不用重复关,但如果你在中间又包了一层FileInputStream并单独持有引用,就容易漏。
内存方面,流式遍历的原则是:不要在循环里累积数据。如果最终结果要返回一个 List,那要评估数据量——一百万条线的坐标点可能有几千万个,每个Coordinate对象几十字节,加起来就是几个 G。这种场景要么改成边处理边落盘,要么用double[]原始数组代替对象数组,内存能省一个数量级。
6.3 大文件的处理策略与实测感受
单文件超过 500MB 的 shp 我处理过不少,最大的一个是 3.2GB,记录的是一条省域范围的线状数据,四百多万条。这种量级下,几个策略的差别非常明显:
| 策略 | 内存峰值 | 耗时(参考) | 适用情况 |
|---|---|---|---|
GeoTools 全量toArray() | 8G 以上 | 慢,容易 OOM | 不推荐 |
GeoToolsSimpleFeatureIterator | 800M 左右 | 中等 | 需要属性或坐标转换 |
| 手写流式解析 | 60M 左右 | 快 | 只要坐标 |
| 手写 + 多线程分片 | 200M 左右 | 最快 | 文件在本地磁盘且可随机读 |
多线程分片的思路是:先用 .shx 拿到每条记录的偏移,把记录均匀切给几个线程,每个线程用RandomAccessFile定位到自己那一区间独立解析。这个方案对磁盘 IO 要求高,机械硬盘上反而可能比单线程慢,固态盘上提速明显。实现时要注意 .shx 的偏移单位是 16 位字,要乘 2 才是字节偏移。
还有一个实践建议:如果 shp 只是中间产物,能在上游就转成更合适的格式就转。比如只取坐标这个需求,先转成二进制自定义格式或者直接入库,后续查询会快很多,也不用反复解析那个混合字节序的文件头。我在一个项目里把每天新增的 shp 转成 Parquet 存起来,查询性能提升了几十倍,代价是多了一次转换的开销,整体算下来还是划算的。
最后分享一个我一直在用的小技巧:写一个ShpInspector小工具类,接收一个路径,输出文件头信息、记录数、ShapeType 分布、坐标范围、前 5 条记录的坐标。每次拿到新的 shp 文件先跑一遍,三十秒就能判断这份数据是不是能用、有没有坐标系问题、是不是空几何。这个工具我维护了几年,每次换项目都直接拷过去,比任何文档都好使。