java调用GDAL及JTS实现生成泰森多边形（Voronoi图）的一种方法

一、关于泰森多边形

1.泰森多边形的特性

2.本文的目的

二、实现思路

1.gdal和jts库的maven坐标

2.jts生成泰森多边形的关键代码

3.使用GDAL读取源文件信息的关键代码

4.使用GDAL将生成的泰森多边形写入文件

三、实现结果

1.实现的效果

2.完整代码示例

一、关于泰森多边形

泰森多边形，又称Voronoi图，是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。

1.泰森多边形的特性

每个泰森多边形内仅含有一个离散点数据；
泰森多边形内的点到相应离散点的距离最近；
位于泰森多边形边上的点到其两边的离散点的距离相等。

2.本文的目的

泰森多边形在地理信息系统（GIS）领域有着广泛的应用，一般可用于点插值，也可以在不直接计算距离的情况进行最邻近分析。目前常见的桌面GIS软件基本都有此功能。本文不讨论泰森多边形的实现算法，仅仅从应用开发的角度出发介绍如何使用已有的java矢量数据读写库、几何处理库来实现泰森多边形的生成。

二、实现思路

GDAL是一种常用的地理空间栅格及矢量数据的读写库其由C/C++编写而成，存在java绑定库，可以被java语言调用。GDAL内置了部分几何处理和空间分析的算法，经作者了解，其暂未内置泰森多边形算法。GDAL在矢量数据的功能方面，支持读写geopackage、shapefile、kml、geojson、gml、xlsx等多种格式的数据。

JTS是一个java语言开发的几何图形处理库，具有较丰富的几何图形处理能力。经了解，JTS内置了泰森多边形算法，即VoronoiDiagramBuilder。基于上述分析，使用java调用GDAL和JTS库实现泰森多边形的生成在技术上是可行的。

1.gdal和jts库的maven坐标

    <dependency><groupId>org.gdal</groupId><artifactId>gdal</artifactId><version>3.2.0</version></dependency><dependency><groupId>org.locationtech.jts</groupId><artifactId>jts-core</artifactId><version>1.18.2</version></dependency>

2.jts生成泰森多边形的关键代码

VoronoiDiagramBuilder voronoiBuilder = new VoronoiDiagramBuilder();
voronoiBuilder.setSites(points);
voronoiBuilder.setClipEnvelope(envelope);
Geometry voronoiPolygons = voronoiBuilder.getDiagram(geometryFactory);

其中points为坐标类的Coordinate列表， envelope为泰森多边形的空间范围定义。

points、envelope均需要从源文件中读取，这时可以使用gdal读取有关信息。

3.使用GDAL读取源文件信息的关键代码

Geometry ogrGeo = f.GetGeometryRef();
Coordinate coordinate = new Coordinate(ogrGeo.GetX(), ogrGeo.GetY());

上面的代码是读取矢量文件坐标点的gdal代码，f为gdal的Feature类。此处的Geometry为gdal的Geometry类，非JTS的Geometry类。下面的代码是读取矢量文件范围的gdal代码。

double[] extent = layer.GetExtent();
Envelope envelope = new Envelope(extent[0],extent[1],extent[2],extent[3]);

4.使用GDAL将生成的泰森多边形写入文件

通过gdal读取矢量文件中调用jts获取泰森多边形需要的参数后，jts生成了Geometry类（org.locationtech.jts.geom.Geometry），逐个读取Geometry的面状图形，将面状图形的坐标点生成wkt文本，利用gdal的根据wkt文本创建几何的功能，将坐标串转换为gdal的几何对象，最后调用gdal写矢量图层的方法将数据生成指定格式的文件即可。以下是转换数据的核心代码。

String wkt = voronoiPolygons.getGeometryN(i).toText();
Geometry geometryOut = Geometry.CreateFromWkt(wkt);

三、实现结果

1.实现的效果

通过本文的示例代码，可以实现最基础的泰森多边形的生成，从而应用于所需要的应用场景中。以下是QGIS中查看代码生成的泰森多边形的实现效果，可以看到，除了生成了多边形外，还继承了点的属性数据，与QGIS自带的泰森多边形功能生成的多边形几何形状一致。详细的示例代码请见下一节。

2.完整代码示例

package com.hjzx.util;import org.gdal.ogr.*;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.triangulate.VoronoiDiagramBuilder;import java.io.File;
import java.util.*;public class Voronoi {public static void main(String[] args) {String vectorPath = "E:\\随机点.gpkg";String outPath = "E:\\泰森多边形.shp";boolean success = createVoronoiByVector(vectorPath, outPath, "ESRI Shapefile");}/*** 生成泰森多边形* @param vectorPath 矢量数据的绝对路径* @param outPath 要输出的泰森多边形绝对路径* @param outDriverName 要输出的数据格式对应的驱动字符串* @return 是否生成成功*/public static boolean createVoronoiByVector(String vectorPath, String outPath, String outDriverName) {ogr.RegisterAll();DataSource dataSource = ogr.Open(vectorPath);System.out.printf("输入的矢量文件路径：%s\t", vectorPath);if (dataSource == null) {System.out.println("打开数据源失败");return false;}Driver driver = dataSource.GetDriver();if (driver == null) {System.out.println("打开驱动失败");return false;}System.out.printf("打开驱动成功：驱动名称：%s\n", driver.getName());//这里简化了操作，只读取了第1个图层Layer layer = dataSource.GetLayer(0);if (layer == null) {System.out.printf("打开矢量图层失败，文件路径：%s", vectorPath);return false;}if (ogrConstants.wkbPoint != layer.GetGeomType()) {System.out.printf("矢量图层不是点状几何图形，无法生成泰森多边形，文件路径：%s", vectorPath);return false;}//gdal获取矢量图层范围double[] extent = layer.GetExtent();// 创建JTS几何工厂GeometryFactory geometryFactory = new GeometryFactory();// 创建一组点List<Coordinate> points = new ArrayList<>();//gdal读取坐标点并赋值给pointsFeature feature;//使用Map集合存储坐标与要素的对应关系Map<Coordinate, Feature> gdalKeyValue = new HashMap<>();while ((feature = layer.GetNextFeature()) != null) {Geometry ogrGeo = feature.GetGeometryRef();//将gdal矢量要素的点转换到JTS中去Coordinate coordinate = new Coordinate(ogrGeo.GetX(), ogrGeo.GetY());points.add(coordinate);gdalKeyValue.put(coordinate, feature);}System.out.printf("输入矢量数据要素转换为坐标点列表成功,要素数量：%d\n", points.size());Envelope envelope = new Envelope(extent[0], extent[1], extent[2], extent[3]);VoronoiDiagramBuilder voronoiBuilder = new VoronoiDiagramBuilder();voronoiBuilder.setSites(points);voronoiBuilder.setClipEnvelope(envelope);org.locationtech.jts.geom.Geometry voronoiPolygons = voronoiBuilder.getDiagram(geometryFactory);int geometryCount = voronoiPolygons.getNumGeometries();System.out.printf("泰森多边形几何图形创建成功,要素数量：%d\n", geometryCount);System.out.println("开始写入泰森多边形矢量文件！");//获取数据数据源的驱动Driver outDriver = ogr.GetDriverByName(outDriverName);if (outDriver == null) {System.out.printf("打开输出文件的驱动失败，驱动字符串：%s\t", outDriverName);return false;}DataSource outDataSource = outDriver.CreateDataSource(outPath);Vector<String> options = new Vector<>();options.add("ENCODING=UTF-8");/*将文件名作为图层名*/File file = new File(outPath);String name = file.getName();name = name.substring(0, name.lastIndexOf("."));//创建输出矢量图层，沿用源坐标系Layer layerOut = outDataSource.CreateLayer(name, layer.GetSpatialRef(), ogr.wkbMultiPolygon, options);//沿用源文件字段FeatureDefn featureDefn = layer.GetLayerDefn();System.out.printf("输入源文件要素定义数量：%d\n", featureDefn.GetFieldCount());FeatureDefn featureDefnOut = new FeatureDefn();//给输出图层创建字段for (int i = 0; i < featureDefn.GetFieldCount(); i++) {featureDefnOut.AddFieldDefn(featureDefn.GetFieldDefn(i));layerOut.CreateField(featureDefn.GetFieldDefn(i));}for (int i = 0; i < geometryCount; i++) {//取出泰森多边形中的每个面对应的Coordinate对象Coordinate coordinate = (Coordinate) voronoiPolygons.getGeometryN(i).getUserData();String wkt = voronoiPolygons.getGeometryN(i).toText();Geometry geometryOut = Geometry.CreateFromWkt(wkt);//根据Coordinate对象到hashMap中获取矢量要素Feature featureOut = gdalKeyValue.get(coordinate);//改变要素的几何图形featureOut.SetGeometry(geometryOut);//创建要素到要输出的图层layerOut.CreateFeature(featureOut);}//保存数据outDataSource.SyncToDisk();//销毁输出数据源layerOut.delete();outDataSource.delete();//销毁输入数据源layer.delete();dataSource.delete();System.out.printf("创建泰森多边形文件成功，文件路径：%s\n", outPath);return true;}
}