概述

在Armap工具箱‘分析工具->疊加分析’，不得不說，非常好用，本文給你講講如何在geotools中實現。

關鍵點

要實現類似的功能有兩個關鍵點：
1、已經計算過的兩個資料不能重複計算；
2、需要保留兩個shp圖形的屬性。
這兩點在後面的程式碼裡面會有相對比較詳細的註釋的。

實現結果

實現程式碼

package com.lzugis.test;

import com.alibaba.fastjson.JSONObject;
import com.vividsolutions.jts.geom.Geometry;
import
 
 com.vividsolutions.jts.geom.MultiPolygon;
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.data.simple.SimpleFeatureCollection;
import
 
 org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import
 
 org.opengis.feature.type.AttributeDescriptor;
import org.opengis.feature.type.AttributeType;
import org.opengis.feature.type.Name;

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;

public class ShapeIntersect {

    public static void main(String[] args){
        double start = System.currentTimeMillis();

        String inputPath1 = "D:\\data\\province.shp",
                inputPath2 = "D:\\data\\test.shp",
                outputPath = "D:\\data\\province_test.shp";

        try {
            File inputFile1 = new File(inputPath1);
            File inputFile2 = new File(inputPath2);

            ShapefileDataStore shpDataStore1 = new ShapefileDataStore(inputFile1.toURL());
            ShapefileDataStore shpDataStore2 = new ShapefileDataStore(inputFile2.toURL());

            //屬性編碼
            Charset charset = Charset.forName("GBK");

            shpDataStore1.setCharset(charset);
            shpDataStore2.setCharset(charset);

            String typeName1 = shpDataStore1.getTypeNames()[0];
            String typeName2 = shpDataStore2.getTypeNames()[0];

            SimpleFeatureSource featureSource1 = shpDataStore1.getFeatureSource(typeName1);
            SimpleFeatureSource featureSource2 = shpDataStore2.getFeatureSource(typeName2);

            SimpleFeatureCollection featureCollection1 = featureSource1.getFeatures();
            SimpleFeatureCollection featureCollection2 = featureSource2.getFeatures();

            /**
             * mapFields記錄的是兩個圖層的屬性名稱，
             *          在處理第二個圖層的時候，如果已經有了這個名稱，
             *          會在欄位後面加‘_1’予以區分
             * fields1為圖層1的欄位
             * fields2為圖層2的欄位
             */
            Map<String, Class> mapFields = new HashMap();
            List<Map> fields1 = new ArrayList(),
                fields2 = new ArrayList();

            SimpleFeatureType featureType1 = featureCollection1.getSchema();
            List<AttributeDescriptor> attrList1 = featureType1.getAttributeDescriptors();
            for(int i=0;i<attrList1.size();i++){
                AttributeDescriptor attr = attrList1.get(i);
                String name = attr.getName().toString();
                Class type = attr.getType().getBinding();
                if(name != "the_geom"){
                    mapFields.put(name, type);
                    Map map = new HashMap();
                    map.put("fieldShp", name);
                    map.put("fieldNew", name);
                    fields1.add(map);
                }
            }
            SimpleFeatureType featureType2 = featureCollection2.getSchema();
            List<AttributeDescriptor> attrList2 = featureType2.getAttributeDescriptors();
            for(int j=0;j<attrList2.size();j++){
                AttributeDescriptor attr = attrList1.get(j);
                String name = attr.getName().toString();
                Class type = attr.getType().getBinding();
                if(name != "the_geom"){
                    String _name = name;
                    if(mapFields.containsKey(name)){
                        _name = _name+"_1";
                    }
                    mapFields.put(_name, type);
                    Map map = new HashMap();
                    map.put("fieldShp", name);
                    map.put("fieldNew", _name);
                    fields2.add(map);
                }
            }

            SimpleFeatureIterator itertor1 = featureCollection1.features();

            //建立輸出檔案
            File outputFile = new File(outputPath);
            Map<String, Serializable> params = new HashMap<String, Serializable>();
            params.put( ShapefileDataStoreFactory.URLP.key, outputFile.toURI().toURL() );
            ShapefileDataStore ds = (ShapefileDataStore) new ShapefileDataStoreFactory().createNewDataStore(params);
            //定義圖形資訊和屬性資訊
            SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
            tb.setCRS(DefaultGeographicCRS.WGS84);
            tb.setName("shapefile");
            tb.add("the_geom", MultiPolygon.class);
            for(String key:mapFields.keySet()){
                tb.add(key, mapFields.get(key));
            }
            ds.createSchema(tb.buildFeatureType());
            //設定編碼
            ds.setCharset(charset);
            //設定Writer
            FeatureWriter<SimpleFeatureType, SimpleFeature> writer = ds.getFeatureWriter(ds.getTypeNames()[0], Transaction.AUTO_COMMIT);

            //記錄已經參與過計算的資料
            Map hasDone = new HashMap();
            //開始計算
            while (itertor1.hasNext()) {
                SimpleFeature feature1 = itertor1.next();
                Geometry geom1 = (Geometry) feature1.getDefaultGeometry();
                String id1 = feature1.getID();

                SimpleFeatureIterator itertor2 = featureCollection2.features();
                while (itertor2.hasNext()){
                    SimpleFeature feature2 = itertor2.next();
                    Geometry geom2 = (Geometry) feature2.getDefaultGeometry();
                    String id2 = feature2.getID();
                    //判斷是否已經參與了計算，需要考慮1∩2和2∩1兩種情況
                    boolean isDone1 = hasDone.containsKey(id1+"-"+id2),
                        isDone2 = hasDone.containsKey(id2+"-"+id1),
                        isIntersect = geom1.intersects(geom2);
                    if(!isDone1 && !isDone2 && isIntersect){
                        Geometry geomOut = geom1.intersection(geom2);
                        SimpleFeature featureOut = writer.next();
                        featureOut.setAttribute("the_geom", geomOut);
                        for(int i=0;i<fields1.size();i++){
                            Map map = fields1.get(i);
                            String fieldShp = map.get("fieldShp").toString(),
                                    fieldNew = map.get("fieldNew").toString();
                            featureOut.setAttribute(fieldNew, feature1.getAttribute(fieldShp));
                        }
                        for(int i=0;i<fields2.size();i++){
                            Map map = fields2.get(i);
                            String fieldShp = map.get("fieldShp").toString(),
                                    fieldNew = map.get("fieldNew").toString();
                            featureOut.setAttribute(fieldNew, feature2.getAttribute(fieldShp));
                        }
                        writer.write();
                    }
                    hasDone.put(id1+"-"+id2, true);
                    hasDone.put(id2+"-"+id1, true);
                }
                itertor2.close();
            }

            writer.close();
            ds.dispose();

            itertor1.close();
        } catch (Exception e){
            e.printStackTrace();
        }

        double end = System.currentTimeMillis();
        System.out.println("共耗時"+(end-start)+"MS");
    }
}

技術部落格
CSDN：http://blog.csdn.NET/gisshixisheng
線上教程
https://edu.csdn.net/course/detail/799
https://edu.csdn.net/course/detail/7471
聯絡方式

“GIS講堂”知識星球開通了，在星球，我將提供一對一的問答服務，你問我答，期待與你相見。