geotools中泰森多邊形的生成

概述

本文講述如何在geotools中生成泰森多邊形，並shp輸出。

泰森多邊形

一、定義

泰森多邊形又叫馮洛諾伊圖（Voronoi diagram），得名於Georgy Voronoi，是由一組由鏈接兩鄰點直線的垂直平分線組成的連續多邊形組成。

二、創建步驟

創建泰森多邊形算法的關鍵是對離散數據點合理地連成三角網，即構建Delaunay三角網。創建泰森多邊形的步驟爲：
1）離散點自動構建三角網，即構建Delaunay三角網。對離散點和造成的三角形編號，記錄每一個三角形是由哪三個離散點構成的。
2）找出與每一個離散點相鄰的全部三角形的編號，並記錄下來。這隻要在已構建的三角網中找出具備一個相同頂點的全部三角形便可。
3）對與每一個離散點相鄰的三角形按順時針或逆時針方向排序，以便下一步鏈接生成泰森多邊形。設離散點爲o。找出以o爲頂點的一個三角形，設爲A；取三角形A除o之外的另外一頂點，設爲a，則另外一個頂點也可找出，即爲f；則下一個三角形必然是以of爲邊的，即爲三角形F；三角形F的另外一頂點爲e，則下一三角形是以oe爲邊的；如此重複進行，直到回到oa邊。
4）計算每一個三角形的外接圓圓心，並記錄之。
5）根據每一個離散點的相鄰三角形，鏈接這些相鄰三角形的外接圓圓心，即獲得泰森多邊形。對於三角網邊緣的泰森多邊形，可做垂直平分線與圖廓相交，與圖廓一塊兒構成泰森多邊形。

三、特徵

1）每一個泰森多邊形內僅含有一個離散點數據；
2）泰森多邊形內的點到相應離散點的距離最近；
3）位於泰森多邊形邊上的點到其兩邊的離散點的距離相等。

geotools中的生成

一、建立測試點

	    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));
	    }

二、生成泰森多邊形

	    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");
	}
}

參考文獻：

一、百度百科

二、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/

聯繫方式

q       q:1004740957

e-mail:niujp08@qq.com

公衆號:lzugis15

Q Q 羣:452117357(webgis)

             337469080(Android)