geotools中泰森多邊形的生成

本文轉載自查看原文 2017-07-23 09:17 3251 ——【GIS二次開發】/ ● 3S技術

概述

本文講述如何在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

---------------------------------------------------------------------------------------------------------------

技術博客

CSDN：http://blog.csdn.NET/gisshixisheng

博客園：http://www.cnblogs.com/lzugis/

在線教程

http://edu.csdn.Net/course/detail/799

Github

https://github.com/lzugis/

聯系方式