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读取坐标点并赋值给points
        Feature 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;

    }
}