geotools中泰森多边形的生成

概述

本文讲述如何在geotools中生成泰森多边形,并shp输出。

泰森多边形

1、定义

泰森多边形又叫冯洛诺伊图(Voronoi diagram),得名于Georgy Voronoi,是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。

2、建立步骤

建立泰森多边形算法的关键是对离散数据点合理地连成三角网,即构建Delaunay三角网。建立泰森多边形的步骤为: 1)离散点自动构建三角网,即构建Delaunay三角网。对离散点和形成的三角形编号,记录每个三角形是由哪三个离散点构成的。 2)找出与每个离散点相邻的所有三角形的编号,并记录下来。这只要在已构建的三角网中找出具有一个相同顶点的所有三角形即可。 3)对与每个离散点相邻的三角形按顺时针或逆时针方向排序,以便下一步连接生成泰森多边形。设离散点为o。找出以o为顶点的一个三角形,设为A;取三角形A除o以外的另一顶点,设为a,则另一个顶点也可找出,即为f;则下一个三角形必然是以of为边的,即为三角形F;三角形F的另一顶点为e,则下一三角形是以oe为边的;如此重复进行,直到回到oa边。 4)计算每个三角形的外接圆圆心,并记录之。 5)根据每个离散点的相邻三角形,连接这些相邻三角形的外接圆圆心,即得到泰森多边形。对于三角网边缘的泰森多边形,可作垂直平分线与图廓相交,与图廓一起构成泰森多边形。

3、特征

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

geotools中的生成

1、创建测试点

	    int xmin = 0, xmax=180;
        int ymin = 0, ymax=90;
        Random random = new Random();
        List<Geometry> geomsPoints = new ArrayList<Geometry>();
	    for(int i=0;i<100;i++){
	    	int x = random.nextInt(xmax)%(xmax-xmin+1) + xmin,
	            y = random.nextInt(ymax)%(ymax-ymin+1) + ymin;
	        Coordinate coord = new Coordinate(x,y,i);
	        coords.add(coord);
	        clipEnvelpoe.expandToInclude(coord);
	        geomsPoints.add(new GeometryFactory().createPoint(coord));
	    }

2、生成泰森多边形

	    voronoiDiagramBuilder.setSites(coords);
	    voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe);
	    Geometry geom = voronoiDiagramBuilder.getDiagram(JTSFactoryFinder.getGeometryFactory());
	    List<Geometry> geoms = new ArrayList<Geometry>();
	    for(int i=0;i<geom.getNumGeometries();i++){
	        geoms.add(geom.getGeometryN(i));
	    }

完整代码如下:

package com.lzugis.geotools;

import java.io.File;
import java.io.Serializable;
import java.nio.charset.Charset;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Random;

import org.geotools.data.FeatureWriter;
import org.geotools.data.Transaction;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.shapefile.ShapefileDataStoreFactory;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Point;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.triangulate.VoronoiDiagramBuilder;

public class TsdbxTest {
	static TsdbxTest tsdbx = new TsdbxTest();
	public void voronoiTest(){
	    VoronoiDiagramBuilder voronoiDiagramBuilder = new VoronoiDiagramBuilder();
	    List<Coordinate> coords = new ArrayList<Coordinate>();
	    Envelope clipEnvelpoe = new Envelope();
	    int xmin = 0, xmax=180;
        int ymin = 0, ymax=90;
        Random random = new Random();
        List<Geometry> geomsPoints = new ArrayList<Geometry>();
	    for(int i=0;i<100;i++){
	    	int x = random.nextInt(xmax)%(xmax-xmin+1) + xmin,
	            y = random.nextInt(ymax)%(ymax-ymin+1) + ymin;
	        Coordinate coord = new Coordinate(x,y,i);
	        coords.add(coord);
	        clipEnvelpoe.expandToInclude(coord);
	        geomsPoints.add(new GeometryFactory().createPoint(coord));
	    }
	    String pointpath = "d:/data/tsdbxpt.shp";
	    tsdbx.writeShape(pointpath,"Point", geomsPoints);

	    voronoiDiagramBuilder.setSites(coords);
	    voronoiDiagramBuilder.setClipEnvelope(clipEnvelpoe);
	    Geometry geom = voronoiDiagramBuilder.getDiagram(JTSFactoryFinder.getGeometryFactory());
	    List<Geometry> geoms = new ArrayList<Geometry>();
	    for(int i=0;i<geom.getNumGeometries();i++){
	        geoms.add(geom.getGeometryN(i));
	    }
	    
	    String polygonpath = "d:/data/tsdbx.shp";
	    tsdbx.writeShape(polygonpath,"Polygon", geoms);
	}
	/**
	 * 
	 * @param filepath
	 * @param geoType
	 * @param geoms
	 */
	public void writeShape(String filepath, String geoType, List<Geometry> geoms) {
		try {
			//创建shape文件对象
			File file = new File(filepath);
			Map<String, Serializable> params = new HashMap<String, Serializable>();
			params.put( ShapefileDataStoreFactory.URLP.key, file.toURI().toURL() );
			ShapefileDataStore ds = (ShapefileDataStore) new ShapefileDataStoreFactory().createNewDataStore(params);
			//定义图形信息和属性信息
			SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
			tb.setCRS(DefaultGeographicCRS.WGS84);
			tb.setName("shapefile");
			if(geoType=="Point"){
				tb.add("the_geom", Point.class);
			}
			else{
				tb.add("the_geom", Polygon.class);
			}
			
			ds.createSchema(tb.buildFeatureType());
			//设置编码
            Charset charset = Charset.forName("GBK");
            ds.setCharset(charset);
			//设置Writer
			FeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0], Transaction.AUTO_COMMIT);
			for(int i=0,len=geoms.size();i<len;i++){
				//写下一条
				SimpleFeature feature = writer.next();
				Geometry geom = geoms.get(i);
				feature.setAttribute("the_geom", geom);
			}
			writer.write();
			writer.close();
			ds.dispose();
		} 
		catch (Exception e) {
			e.printStackTrace();
		}
	}
	
	public static void main(String[] args){
		long start = System.currentTimeMillis();
		tsdbx.voronoiTest();
		System.out.println("共耗时"+(System.currentTimeMillis() - start)+"ms");
	}
}

参考文献:

1、百度百科

2、geotools user guide

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏服务端技术杂谈

算出实体箱子之间可以存多少水

出题 ? 下午,群里面有同学出了一个简单的算法题,意思是一个房间内,有多个1立方米的箱子,多个箱子可以垂直落在一起,问:剩下的空间可以存多少立方的水(如图)。 ...

3347
来自专栏C语言及其他语言

[每日一题]方程的根

今天的每日一题是大家小学、初中、高中、大学都需要会的一种数学题,但只要我们会了代码,一切都只要输入数据就行,答案秒出,是不是简单了很多呢 题目描述 求方程 的...

2513
来自专栏深度学习自然语言处理

一分钟搞懂的算法之BPE算法

昨天总结实验数据分析的时候发现一个机器翻译的其中的一个脚本,其中用到的算法就是BPE算法,刚开始感觉很高大上的,因为总是听到带上算法帽子的东西就觉得666。等自...

2.2K8
来自专栏PPV课数据科学社区

【学习】《R实战》读书笔记(第四章)

读书会是一种在于拓展视野、宏观思维、知识交流、提升生活的活动。PPV课R语言读书会以“学习、分享、进步”为宗旨,通过成员协作完成R语言专业书籍的精读和分享,达到...

2975
来自专栏美团技术团队

函数式编程在Redux/React中的应用

本文简述了软件复杂度问题及应对策略:抽象和组合;展示了抽象和组合在函数式编程中的应用;并展示了Redux/React在解决前端状态管理的复杂度方面对上述理论的实...

3669
来自专栏数说工作室

统计师的Python日记【第十天:数据聚合】

本文是【统计师的Python日记】第10天的日记 回顾一下: 第1天学习了Python的基本页面、操作,以及几种主要的容器类型。 第2天学习了python的函数...

5298
来自专栏C语言及其他语言

[每日一题]杨辉三角(1231)

题目描述 还记得中学时候学过的杨辉三角吗?具体的定义这里不再描述,你可以参考以下的图形: 1 1 1 1 2 1 1 3 3 1 1 4 6 4 1 ...

3437
来自专栏趣谈编程

冒泡排序

面试官: 写一个冒泡排序吧 冒泡排序是一个比较经典和简单的排序算法,今天我们从从算法本身,时间复杂度以及稳定性方面来看看冒泡排序,这些方面也是研究其他排序算法的...

20110
来自专栏java学习

每日一练(2017/5/22)

Java基础 | 数据库 | Android | 学习视频 | 学习资料下载 课前导读 ●回复"每日一练"获取以前的题目! ●答案公布时间:为每期发布题目的第二...

3428
来自专栏颇忒脱的技术博客

Elasticsearch中将Doc根据A字段排序获得第一个Doc的B字段值的方法

最近遇到这样一个需求,要通过Elasticsearch将Doc根据A字段降序,然后获得B字段的值,最终根据B字段的值再去做Pipeline Aggregatio...

1782

扫码关注云+社区

领取腾讯云代金券