美文网首页
2021-10-28 实现等值面的一些问题

2021-10-28 实现等值面的一些问题

作者: MrSwilder | 来源:发表于2021-10-28 16:41 被阅读0次

pom依赖

<?xml version="1.0" encoding="UTF-8"?>
<project xmlns="http://maven.apache.org/POM/4.0.0"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 https://maven.apache.org/xsd/maven-4.0.0.xsd">
    <modelVersion>4.0.0</modelVersion>
    <parent>
        <groupId>org.springframework.boot</groupId>
        <artifactId>spring-boot-starter-parent</artifactId>
        <version>2.2.1.RELEASE</version>
        <relativePath /> <!-- lookup parent from repository -->
    </parent>
    <groupId>com.summit</groupId>
    <artifactId>summit-tools-gis</artifactId>
    <version>0.0.1-SNAPSHOT</version>
    <name>summit-tools-gis</name>
    <description>gis工具</description>

    <properties>
        <java.version>1.8</java.version>
        <geotools.version>26.0</geotools.version>
        <swagger.version>2.7.0</swagger.version>
        <mybatisplus.version>3.3.0</mybatisplus.version>
        <druid.version>1.1.16</druid.version>
    </properties>

    <dependencies>
        <dependency>
            <groupId>org.springframework.boot</groupId>
            <artifactId>spring-boot-starter-web</artifactId>
        </dependency>
        <!-- druid数据库连接池 -->
        <dependency>
            <groupId>com.alibaba</groupId>
            <artifactId>druid-spring-boot-starter</artifactId>
            <version>${druid.version}</version>
        </dependency>
        <dependency>
            <groupId>com.baomidou</groupId>
            <artifactId>mybatis-plus-boot-starter</artifactId>
            <version>${mybatisplus.version}</version>
            <exclusions>
                <exclusion>
                    <groupId>com.baomidou</groupId>
                    <artifactId>mybatis-plus-generator</artifactId>
                </exclusion>
            </exclusions>
        </dependency>
        <dependency>
            <groupId>org.springframework.boot</groupId>
            <artifactId>spring-boot-devtools</artifactId>
            <scope>runtime</scope>
            <optional>true</optional>
        </dependency>
        <!-- pg数据库驱动 -->
        <dependency>
            <groupId>org.postgresql</groupId>
            <artifactId>postgresql</artifactId>
            <scope>runtime</scope>
        </dependency>

        <dependency>
            <groupId>org.projectlombok</groupId>
            <artifactId>lombok</artifactId>
            <optional>true</optional>
        </dependency>
        <dependency>
            <groupId>org.springframework.boot</groupId>
            <artifactId>spring-boot-starter-test</artifactId>
            <scope>test</scope>
            <exclusions>
                <exclusion>
                    <groupId>org.junit.vintage</groupId>
                    <artifactId>junit-vintage-engine</artifactId>
                </exclusion>
            </exclusions>
        </dependency>
        <!-- fastjson -->
        <dependency>
            <groupId>com.alibaba</groupId>
            <artifactId>fastjson</artifactId>
            <version>1.2.62</version>
        </dependency>

        <!--&lt;!&ndash;IDW 插值 &ndash;&gt;-->
        <dependency>
            <groupId>com.wContour</groupId>
            <artifactId>wContour</artifactId>
            <version>2</version>
            <scope>system</scope>
            <systemPath>${project.basedir}/src/main/resources/lib/wContour.jar</systemPath>
        </dependency>
        <!--swagger2 -->
        <dependency>
            <groupId>io.springfox</groupId>
            <artifactId>springfox-swagger2</artifactId>
            <version>${swagger.version}</version>
        </dependency>
        <dependency>
            <groupId>io.springfox</groupId>
            <artifactId>springfox-swagger-ui</artifactId>
            <version>${swagger.version}</version>
        </dependency>

        <!-- GEOTOOLS相关的依赖 -->
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-main</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-shapefile</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-epsg-hsql</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-swing</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-geojson</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-coverage</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-metadata</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-referencing</artifactId>
            <version>${geotools.version}</version>
        </dependency>

        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-grid</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <!-- https://mvnrepository.com/artifact/org.locationtech.jts/jts-core -->
        <dependency>
            <groupId>org.locationtech.jts</groupId>
            <artifactId>jts-core</artifactId>
            <version>1.18.2</version>
        </dependency>

        <!--<dependency>-->
            <!--<groupId>com.vividsolutions</groupId>-->
            <!--<artifactId>jts</artifactId>-->
            <!--<version>1.13</version>-->
        <!--</dependency>-->
        <!-- https://mvnrepository.com/artifact/wContour/wContour -->
        <!--<dependency>-->
            <!--<groupId>wContour</groupId>-->
            <!--<artifactId>wContour</artifactId>-->
            <!--<version>1.6.1</version>-->
        <!--</dependency>-->




    </dependencies>
    <repositories>
        <repository>
            <id>osgeo</id>
            <name>OSGeo Release Repository</name>
            <url>https://repo.osgeo.org/repository/release/</url>
            <snapshots><enabled>false</enabled></snapshots>
            <releases><enabled>true</enabled></releases>
        </repository>
        <repository>
            <id>osgeo-snapshot</id>
            <name>OSGeo Snapshot Repository</name>
            <url>https://repo.osgeo.org/repository/snapshot/</url>
            <snapshots><enabled>true</enabled></snapshots>
            <releases><enabled>false</enabled></releases>
        </repository>
        <repository>
            <id>maven</id>
            <name>maven</name>
            <url>http://mvnrepository.com/</url>
            <snapshots><enabled>true</enabled></snapshots>
            <releases><enabled>false</enabled></releases>
        </repository>

    </repositories>
    <build>
        <plugins>
            <plugin>
                <groupId>org.springframework.boot</groupId>
                <artifactId>spring-boot-maven-plugin</artifactId>
                <configuration>
                    <includeSystemScope>true</includeSystemScope>
                </configuration>
            </plugin>
        </plugins>
    </build>

</project>

核心代码

/**  
 * All rights Reserved, Designed By www.summit.com.cn
 */
package com.summit.gis.modules.iso.utils;

import java.io.File;
import java.io.IOException;
import java.io.StringWriter;
import java.nio.charset.Charset;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import com.summit.gis.modules.iso.vo.Station;
import org.geotools.data.DataUtilities;
import org.geotools.data.FeatureSource;
import org.geotools.data.collection.ListFeatureCollection;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.geotools.feature.SchemaException;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.geojson.feature.FeatureJSON;
import org.json.simple.JSONArray;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.opengis.feature.Feature;
import org.opengis.feature.simple.SimpleFeature;

import com.alibaba.fastjson.JSONObject;
import com.summit.gis.modules.iso.enums.GeometryEnum;
import com.summit.gis.modules.iso.vo.InterpolationVO;
//import com.vividsolutions.jts.geom.Coordinate;
//import com.vividsolutions.jts.geom.Geometry;

import org.opengis.feature.simple.SimpleFeatureType;
import wContour.Contour;
import wContour.Global.PointD;
import wContour.Global.PolyLine;
import wContour.Contour;
import wContour.Global.PointD;
import wContour.Global.PolyLine;
import wContour.Interpolate;

import wContour.Contour;
import wContour.Interpolate;
import wContour.Global.Border;
import wContour.Global.PointD;
import wContour.Global.PolyLine;
import wContour.Global.Polygon;

/**
 * @ClassName: IsoSurfaceUtil
 * @Description: 生成等值面的工具
 * @author: zhangchaojia
 * @date: 2019年12月2日 上午10:13:31
 * 
 */
public class IsoSurfaceUtil {

    private static final double UNDEF_DATA = -9999.0;

    private static final double EPSILON = 0.001;

    /**
     * 
     * <p>Title: calculateIsoSurfaceByIDW</p>
     * 
     * @param vo
     * @return
     * @throws IOException
     */
    public static String calculateIsoSurfaceByIDW(InterpolationVO vo) throws IOException, SchemaException {
        boolean validate = paramsValidate(vo);
        if (validate) {
//          File file = new File(vo.getBoundaryFilePath());
            String geoJson;
            Integer[] gridSize = vo.getGridSize();
            int width = gridSize[0];
            int height = gridSize[1];

            double[] X = new double[width];
            double[] Y = new double[height];
            //------------根据geojson生成simpleFeatureCollection--
            JSONObject boundary=vo.getBoundary();
            FeatureJSON featureJSON=new FeatureJSON();
            FeatureCollection featureCollection=featureJSON.readFeatureCollection(boundary.toString());
            // 获取边界
            double minX = featureCollection.getBounds().getMinX();
            double minY = featureCollection.getBounds().getMinY();
            double maxX = featureCollection.getBounds().getMaxX()+1;
            double maxY = featureCollection.getBounds().getMaxY()+1;

            // 创建网格X / Y坐标
            Interpolate.CreateGridXY_Num(minX, minY, maxX, maxY, X, Y);

            double[] dataInterval = vo.getDataInterval();
            int nc = dataInterval.length;
            //获取站点数据
            List<Station> list=vo.getStationList();
            double[][] data =new double[3][list.size()];
            if(list.size()<3){
                return null;
            }

            for(int i=0;i<list.size();i++){
                Station station=list.get(i);
                data[0][i]=station.getX();
                data[1][i]=station.getY();
                data[2][i]=station.getValue();
            }
            // IDW 插值
            double[][] gridData = Interpolate.Interpolation_IDW_Neighbor(data, X, Y, dataInterval.length, UNDEF_DATA);
            int[][] S1 = new int[gridData.length][gridData[0].length];

            // 用未定义的数据跟踪网格数据的轮廓边界。 网格数据从左到右，从下到上。 网格数据数组：第一个尺寸为Y，第二个尺寸为X。
            List<Border> borders = Contour.tracingBorders(gridData, X, Y, S1, UNDEF_DATA);

            // 轮廓边界
            List<PolyLine> polylineList = Contour.tracingContourLines(gridData, X, Y, nc, dataInterval, UNDEF_DATA,
                    borders, S1);

            // 平滑曲线
            polylineList = Contour.smoothLines(polylineList);
            GeometryEnum type = vo.getType();
            SimpleFeatureCollection contourSimpleFeatureCollection,clipSimpleFeatureCollection;
            switch (type) {
            case Line:

                geoJson = polylineToGeoJson(polylineList);
                //获取等值面要素
                contourSimpleFeatureCollection=getSimpleCollectionByGeoJSON(geoJson);
                // 根据边界裁剪等值面
                clipSimpleFeatureCollection=clipPolylineFeatureCollection(featureCollection,contourSimpleFeatureCollection);
                return SimpleCollectionToGeoJSON(clipSimpleFeatureCollection);

            case Polygon:
                // 根据等值线描述出等值面数据
                List<Polygon> polygonList = Contour.tracingPolygons(gridData, polylineList, borders, dataInterval);
                geoJson = polygonToGeoJson(polygonList);
                //获取等值面要素
                contourSimpleFeatureCollection=getSimpleCollectionByGeoJSON(geoJson);
                // 根据边界裁剪等值面
                clipSimpleFeatureCollection=clipPolygonFeatureCollection(featureCollection,contourSimpleFeatureCollection);

                return SimpleCollectionToGeoJSON(clipSimpleFeatureCollection);
            default:
                geoJson = null;
                break;
            }


        }
        return null;
    }

    /**
     * 通过geojson获取要素类
     * @param geojsom
     * @return
     * @throws IOException
     */
    private static SimpleFeatureCollection getSimpleCollectionByGeoJSON(String geojsom) throws IOException {
        FeatureJSON featureJSON=new FeatureJSON();
        FeatureCollection featureCollection=featureJSON.readFeatureCollection(geojsom);
        List<SimpleFeature> list=new ArrayList<>();
        FeatureIterator featureIterator=featureCollection.features();
        while (featureIterator.hasNext()){
            SimpleFeature feature=(SimpleFeature) featureIterator.next();
            list.add(feature);
        }
        SimpleFeatureCollection simpleFeatureCollection=new ListFeatureCollection((SimpleFeatureType) featureCollection.getSchema(),list);
        return simpleFeatureCollection;
    }


    /**
     * 要素类转geojson
     * @param simpleFeatureCollection
     * @return
     * @throws IOException
     */
    private static String SimpleCollectionToGeoJSON(SimpleFeatureCollection simpleFeatureCollection) throws IOException {
        FeatureCollection featureCollection=(FeatureCollection) simpleFeatureCollection;
        FeatureJSON featureJSON=new FeatureJSON();
        StringWriter writer = new StringWriter();
        featureJSON.writeFeatureCollection(featureCollection,writer);
        return writer.toString();

    }

    /**
     * Contour面转geojson
     * @param polygonList
     * @return
     */
    private static String polygonToGeoJson(List<Polygon> polygonList) {
        String geo;
        String geometry = "{ \"type\":\"Feature\",\"geometry\":";
        String properties = ",\"properties\":{ \"hvalue\":";

        String head = "{\"type\": \"FeatureCollection\"," + "\"features\":[";
        String end = "]}";

        if (polygonList == null || polygonList.isEmpty()) {
            return null;
        }

        StringBuilder geoBuilder = new StringBuilder();

        for (Polygon pPolygon : polygonList) {
            List<Object> ptsTotal = new ArrayList<>();
            List<Object> pts = new ArrayList<>();
            pPolygon.OutLine.PointList.stream().map(point -> {
                List<Double> pt = new ArrayList<>();
                pt.add(point.X);
                pt.add(point.Y);
                return pt;
            }).forEachOrdered(pts::add);
            ptsTotal.add(pts);

            if (pPolygon.HasHoles()) {
                // 有孔洞的多面
                pPolygon.HoleLines.stream().map(cptLine -> {
                    List<Object> cpts = new ArrayList<>();

                    cptLine.PointList.stream().map(ccptD -> {
                        List<Double> pt = new ArrayList<>();
                        pt.add(ccptD.X);
                        pt.add(ccptD.Y);
                        return pt;
                    }).forEachOrdered(cpts::add);
                    return cpts;
                }).filter(cpts -> !cpts.isEmpty()).forEachOrdered(ptsTotal::add);
            }

            JSONObject js = new JSONObject();

            js.put("type", "Polygon");
            js.put("coordinates", ptsTotal);

            double hv = pPolygon.HighValue;
            double lv = pPolygon.LowValue;

//          if (Math.abs(hv - lv) < EPSILON) {
//              // 是否顺时针
//              if (pPolygon.IsClockWise) {
//                  // 是否居中
//                  if (!pPolygon.IsHighCenter) {
//                      hv = hv + 0.1;
//                      lv = lv + 0.1;
//                  }
//              } else {
//                  if (!pPolygon.IsHighCenter) {
//                      hv = hv - 0.1;
//                      lv = lv - 0.1;
//                  }
//              }
//          } else {
//              if (!pPolygon.IsClockWise) {
//                  lv = lv + 0.1;
//              } else {
//                  if (pPolygon.IsHighCenter) {
//                      hv = hv - 0.1;
//                  }
//              }
//
//          }
            if (Math.abs(hv - lv) < EPSILON) {
                if(pPolygon.IsHighCenter){
                    hv=hv+0.1;
                    lv=lv+0.1;
                }else{
                    hv=hv-0.1;
                    lv=lv-0.1;
                }
            }else{
                hv=hv-0.1;
                lv=lv+0.1;
            }

            geoBuilder.insert(0, geometry + js.toString() + properties + hv + ", \"lvalue\":" + lv + "}}" + ",");
        }
        geo = geoBuilder.toString();

        if (geo.contains(",")) {
            geo = geo.substring(0, geo.lastIndexOf(","));
        }

        geo = head + geo + end;

        return geo;
    }



    /**
     * 
     * <p>Title: polylineToGeoJson</p>
     * <p>Description: polyline 转换为GeoJSON字符串</p>
     * 
     * @param polylineList
     * @return
     */
    private static String polylineToGeoJson(List<PolyLine> polylineList) {
        StringBuilder geo = new StringBuilder();
        String geometry = " { \"type\":\"Feature\",\"geometry\":";
        String properties = ",\"properties\":{ \"value\":";

        String head = "{\"type\": \"FeatureCollection\"," + "\"features\": [";
        String end = "  ] }";

        if (polylineList.isEmpty()) {
            return null;
        }

        for (PolyLine pPolyline : polylineList) {
            List<Object> ptsTotal = new ArrayList<>();

            for (PointD ptD : pPolyline.PointList) {
                List<Double> pt = new ArrayList<>();
                pt.add(ptD.X);
                pt.add(ptD.Y);
                ptsTotal.add(pt);
            }

            JSONObject js = new JSONObject();
            js.put("type", "LineString");
            js.put("coordinates", ptsTotal);

            geo.insert(0, geometry + js.toString() + properties + pPolyline.Value + "} }" + ",");
        }

        if (geo.toString().contains(",")) {
            geo = new StringBuilder(geo.substring(0, geo.lastIndexOf(",")));
        }

        geo = new StringBuilder(head + geo + end);
        return geo.toString();
    }



    /**
     *
     * <p>Title: paramsValidate</p>
     * <p> Description: 参数校验</p>
     *
     * @param vo
     * @return
     */
    private static boolean paramsValidate(InterpolationVO vo) {

        return true;
    }

    /**
     * 裁剪图形
     * @param fc
     * @param gs
     * @return
     * @throws SchemaException
     */
    private static SimpleFeatureCollection clipPolygonFeatureCollection(FeatureCollection fc,
                                                SimpleFeatureCollection gs) throws SchemaException {
        SimpleFeatureCollection simpleFeatureCollection = null;
        final SimpleFeatureType TYPE= DataUtilities.createType("polygons", "the_geom:MultiPolygon:srid=4326,lvalue:double,hvalue:double");
        try {
            List<SimpleFeature> list=new ArrayList<>();
            SimpleFeatureIterator contourFeatureIterator = gs.features();
            FeatureIterator dataFeatureIterator = fc.features();
            while (dataFeatureIterator.hasNext()) {
                SimpleFeature dataFeature =(SimpleFeature) dataFeatureIterator.next();
                Geometry dataGeometry = (Geometry)dataFeature.getDefaultGeometry();
                while (contourFeatureIterator.hasNext()) {
                    SimpleFeature contourFeature = contourFeatureIterator.next();
                    Geometry contourGeometry = (Geometry) contourFeature.getDefaultGeometry();
                    double lv = (Double) contourFeature.getProperty("lvalue")
                            .getValue();
                    double hv = (Double) contourFeature.getProperty("hvalue")
                            .getValue();
                    try{
                        if(dataGeometry.getGeometryType()=="MultiPolygon"){
                            for(int i=0;i<dataGeometry.getNumGeometries();i++){
                                Geometry geom=dataGeometry.getGeometryN(i);
                                if (geom.intersects(contourGeometry)) {
                                    Geometry geo = geom
                                            .intersection(contourGeometry);

                                    SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
                                    featureBuilder.add(geo);
                                    featureBuilder.add(lv);
                                    featureBuilder.add(hv);
                                    SimpleFeature feature =featureBuilder.buildFeature(null);
                                    list.add(feature);
//                      }

                                }
                            }
                        }else {
                            if (dataGeometry.intersects(contourGeometry)) {
                                Geometry geo = dataGeometry
                                        .intersection(contourGeometry);

                                SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);
                                featureBuilder.add(geo);
                                featureBuilder.add(lv);
                                featureBuilder.add(hv);
                                SimpleFeature feature = featureBuilder.buildFeature(null);
                                list.add(feature);
                            }
                        }


                    }catch (Exception e){
                        System.out.println(e.getMessage());
                    }finally {
                        continue;
                    }


                }

            }

            contourFeatureIterator.close();
            dataFeatureIterator.close();
            simpleFeatureCollection=new ListFeatureCollection(TYPE,list);

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

        return simpleFeatureCollection;
    }

    /**
     * 裁剪线
     * @param fc
     * @param gs
     * @return
     * @throws SchemaException
     */
    private static SimpleFeatureCollection clipPolylineFeatureCollection(FeatureCollection fc,
                                                                        SimpleFeatureCollection gs) throws SchemaException {
        SimpleFeatureCollection simpleFeatureCollection = null;
        final SimpleFeatureType TYPE= DataUtilities.createType("polygons", "the_geom:LineString:srid=4326,value:double");
        try {
            List<SimpleFeature> list=new ArrayList<>();
            SimpleFeatureIterator contourFeatureIterator = gs.features();
            FeatureIterator dataFeatureIterator = fc.features();
            while (dataFeatureIterator.hasNext()) {
                SimpleFeature dataFeature =(SimpleFeature) dataFeatureIterator.next();
                Geometry dataGeometry = (Geometry)dataFeature.getDefaultGeometry();
                while (contourFeatureIterator.hasNext()) {
                    SimpleFeature contourFeature = contourFeatureIterator.next();
                    Geometry contourGeometry = (Geometry) contourFeature.getDefaultGeometry();

                    double value = (Double) contourFeature.getProperty("value")
                            .getValue();
                    if (dataGeometry.intersects(contourGeometry)) {
                        Geometry geo = dataGeometry
                                .intersection(contourGeometry);
                        SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);

                        featureBuilder.add(geo);
                        featureBuilder.add(value);
                        SimpleFeature feature =featureBuilder.buildFeature(null);
                        list.add(feature);
                    }
                }
            }
            contourFeatureIterator.close();
            dataFeatureIterator.close();
            simpleFeatureCollection=new ListFeatureCollection(TYPE,list);

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

        return simpleFeatureCollection;
    }
}

! wcontour依赖包不要用maven下载的

相关文章

  • 2021-10-28 实现等值面的一些问题

    pom依赖 核心代码 ! wcontour依赖包不要用maven下载的

  • Geotools中实现NC转等值面

    概述： 前面的文章有实现IDW插值并生成等值面的，本文在前文基础上实现气象NC数据生成等值面。 效果： Arcgi...

  • geotools等值线生成

    概述 前文中，提到了等值面的生成,后面有人经常会问等值线的生成，本文在前文的基础上做了一点修改，完成了等值线的ge...

  • 多表查询

    等值连接 表别名 多表等值连接 自然连接、USING子句、ON子句 自关联 非等值连接 等值连接   其中sele...

  • （六）Hive函数大全

    一、关系运算： 1. 等值比较: = 2. 等值比较:<=> 3. 不等值比较: <>和!= 4. 小于比较: <...

  • element ui checkbox 拿到想要的数据

    需求：从对象中拿到label id等值,以便传值给后面的元素（对象）