目录
一、关于泰森多边形
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;}
}