❤️利用geoTools计算shp面积❤️

1
2
3
4
5
6
7
8
9
10
11
<dependency>
    <groupId>org.geotools</groupId>
    <artifactId>gt-shapefile</artifactId>
    <version>18.4</version>
</dependency>
<dependency>
    <groupId>org.geotools</groupId>
    <artifactId>gt-epsg-hsql</artifactId>
    <version>18.4</version>
</dependency>

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

import com.vividsolutions.jts.geom.Geometry;
import org.geotools.data.Query;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.shapefile.ShapefileDataStoreFactory;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.opengis.feature.Property;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import java.io.File;
import java.io.IOException;
import java.net.MalformedURLException;
import java.nio.charset.Charset;

/**
 * @Author: LEAVES
 * @Version 1.0
 * @Date: 2021年09月13日  10时55分08秒
 * @Description:    利用geoTools计算shp面积
 * <p>
 * 参考：https://docs.geotools.org/latest/tutorials/geometry/geometrycrs.html#epsg-codes
 */
public class CalculateShpArea {

    private static Logger log = LoggerFactory.getLogger(GeoToolsUtil.class);

    private final static String EPSG_CODE = "EPSG:3857";

    /**
     * 利用geoTools计算shp面积
     *
     * @param shpFilePath
     * @return
     */
    public static double getShpArea(String shpFilePath) {
        //此时计算出的面积单位是平方米
        double shpArea = getArea(shpFilePath);
        //换算成亩
        shpArea = shpArea / 666.7;
        return shpArea;
    }

    /**
     * 计算shp面积
     *
     * @param shpFilePath
     * @return
     */
    public static double getArea(String shpFilePath) {
        //定义面积
        double shpTotalArea = 0;
        //读取shp文件
        ShapefileDataStore shpStore = buildDataStore(shpFilePath);
        try {
            //获取shp文件坐标系
            CoordinateReferenceSystem srcCRS = CRS.parseWKT(getCRSWkt(shpFilePath));
            SimpleFeatureSource source = shpStore.getFeatureSource();
            SimpleFeatureType schema = source.getSchema();
            Query query = new Query(schema.getTypeName());
            FeatureCollection<SimpleFeatureType, SimpleFeature> collection = source.getFeatures(query);
            try (FeatureIterator<SimpleFeature> features = collection.features()) {
                while (features.hasNext()) {
                    SimpleFeature feature = features.next();
                    log.info(feature.getID() + ": ");
                    for (Property attribute : feature.getProperties()) {
                        log.info("t" + attribute.getName() + ":" + attribute.getValue());
                    }
                    Geometry geometry = (Geometry) feature.getDefaultGeometry();
                    geometry.setSRID(4326);
                    CoordinateReferenceSystem crsTarget = CRS.decode(EPSG_CODE);
                    // 投影转换
                    MathTransform transform = CRS.findMathTransform(srcCRS, crsTarget);
                    Geometry res = JTS.transform(geometry, transform);
                    //可以打印看下是否是个面（起点必须和结束点重合）
                    for (Coordinate coordinate : res.getCoordinates()) {
                        log.info("coordinate = " + coordinate);
                    }
                    //累加每一个面的面积
                    shpTotalArea = shpTotalArea + res.getArea();
                }
            }
        } catch (IOException e) {
            e.printStackTrace();
        } catch (FactoryException e) {
            e.printStackTrace();
        } catch (TransformException e) {
            e.printStackTrace();
        }
        return shpTotalArea;
    }

    /**
     * 构建ShapeDataStore对象。
     *
     * @param shpFilePath shape文件路径。
     * @return
     */
    public static ShapefileDataStore buildDataStore(String shpFilePath) {
        ShapefileDataStoreFactory factory = new ShapefileDataStoreFactory();
        try {
            ShapefileDataStore dataStore = (ShapefileDataStore) factory
                    .createDataStore(new File(shpFilePath).toURI().toURL());
            if (dataStore != null) {
                //这儿根据shp数据实际编码来自己确定
                dataStore.setCharset(Charset.forName("UTF-8"));
//                dataStore.setCharset(Charset.forName("GBK"));
            }
            return dataStore;
        } catch (MalformedURLException e) {
            e.printStackTrace();
        } catch (IOException e) {
            e.printStackTrace();
        }
        return null;
    }

    /**
     * 读取WKT格式的坐标系信息
     * <p>
     * shape文件并不是一个文件而是一堆文件，坐标系信息存储在*.prj文件，*.prj文件中存储的直接就是WKT格式的坐标系信息。
     * <p>
     * 也可以直接读取*.prj来获取对应的WKT格式数据
     *
     * @param shpFilePath
     * @return
     */
    public static String getCRSWkt(String shpFilePath) {
        ShapefileDataStore dataStore = buildDataStore(shpFilePath);
        String wkt = null;
        try {
            wkt = dataStore.getSchema().getCoordinateReferenceSystem().toWKT();
        } catch (IOException e) {
            e.printStackTrace();
        }
        return wkt;
    }
}