修改栅格分析的线分析功能，添加GeoTools工具类

13693261870

2023-09-14 3b74a18403fa5004e0b3820b0c171d90a615ca21

修改栅格分析的线分析功能，添加GeoTools工具类

 pom.xml

@@ -70,18 +70,6 @@
            <artifactId>mybatis-plus-boot-starter</artifactId>
            <version>3.5.1</version>
        </dependency>
        <!--支持PGIS的geometry类型：https://www.91mszl.com/zhangwuji/article/details/1350-->
        <!--<dependency>
            <groupId>com.eyougo</groupId>
            <artifactId>mybatis-typehandlers-postgis</artifactId>
            <version>1.0</version>
        </dependency>-->
        <!--mybatis-->
        <!--<dependency>
            <groupId>org.mybatis.spring.boot</groupId>
            <artifactId>mybatis-spring-boot-starter</artifactId>
            <version>2.2.0</version>
        </dependency>-->
        <!--redis-->
        <dependency>
            <groupId>org.springframework.boot</groupId>
@@ -113,7 +101,7 @@
            <groupId>org.springframework.boot</groupId>
            <artifactId>spring-boot-starter-aop</artifactId>
        </dependency>
        <!--页面热部署-->
        <!--devtools-->
        <dependency>
            <groupId>org.springframework.boot</groupId>
            <artifactId>spring-boot-devtools</artifactId>
@@ -126,7 +114,7 @@
            <artifactId>fastjson</artifactId>
            <version>1.2.47</version>
        </dependency>
        <!--整合druid数据源-->
        <!--druid-->
        <dependency>
            <groupId>com.alibaba</groupId>
            <artifactId>druid</artifactId>
@@ -265,6 +253,28 @@
            <artifactId>sqlite-jdbc</artifactId>
            <version>3.36.0.3</version>
        </dependency>
        <!--geotools-->
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-main</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-epsg-hsql</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-api</artifactId>
            <version>20.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-opengis</artifactId>
            <version>22.0</version>
            <scope>compile</scope>
        </dependency>
    </dependencies>

    <build>
@@ -280,4 +290,29 @@
            </plugin>
        </plugins>
    </build>

    <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>
    </repositories>
</project>

 src/main/java/com/moon/server/helper/GdalHelper.java

@@ -53,6 +53,7 @@
        // æ³¨å†Œæ‰€æœ‰çš„驱动
        gdal.AllRegister();
        ogr.RegisterAll();
        GeoHelper.initSr();
    }

    /**

 src/main/java/com/moon/server/helper/GeoHelper.java

¶Ô±ÈÐÂÎļþ
@@ -0,0 +1,112 @@
package com.moon.server.helper;

import com.moon.server.entity.all.StaticData;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import org.gdal.osr.SpatialReference;
import org.geotools.geometry.DirectPosition2D;
import org.geotools.referencing.CRS;
import org.geotools.referencing.GeodeticCalculator;
import org.opengis.referencing.crs.CoordinateReferenceSystem;

import java.awt.geom.Point2D;

/**
 * Geo帮助类
 * @author WWW
 * @date 2023-09-14
 */
public class GeoHelper {
    public static SpatialReference sr104903;

    public final static double MOON_RADIUS = 1738000;

    public static CoordinateReferenceSystem crs104903;

    private final static Log log = LogFactory.getLog(GeoHelper.class);

    /**
     * åˆå§‹åŒ–坐标系
     */
    public static void initSr() {
        try {
            sr104903 = new SpatialReference(StaticData.MOON_2000_WKT);

            crs104903 = CRS.parseWKT(StaticData.MOON_2000_WKT);
        } catch (Exception ex) {
            log.error(ex.getMessage(), ex);
        }
    }

    /**
     * èŽ·å–距离1
     */
    public static double getDistance1(double lon1, double lat1, double lon2, double lat2) {
        double radLat1 = Math.toRadians(lat1);
        double radLat2 = Math.toRadians(lat2);
        double a = radLat1 - radLat2;
        double b = Math.toRadians(lon1) - Math.toRadians(lon2);
        double s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2) + Math.cos(radLat1) * Math.cos(radLat2) * Math.pow(Math.sin(b / 2), 2)));

        return s * MOON_RADIUS;
    }

    /**
     * èŽ·å–距离2
     */
    public static double getDistance2(double x1, double y1, double x2, double y2) {
        GeodeticCalculator geodeticCalculator = new GeodeticCalculator(crs104903);
        geodeticCalculator.setStartingGeographicPoint(x1, y1);
        geodeticCalculator.setDestinationGeographicPoint(x2, y2);

        return geodeticCalculator.getOrthodromicDistance();
    }

    /**
     * èŽ·å–方向角
     */
    public static double getBearing(double x1, double y1, double x2, double y2) {
        return (90 - Math.toDegrees(Math.atan2(x2 - x1, y2 - y1)) + 360) % 360;
    }

    /**
     * èŽ·å–角度1
     */
    public static double getAngle1(double x1, double y1, double x2, double y2) {
        DirectPosition2D p1 = new DirectPosition2D(crs104903, x1, y1);
        DirectPosition2D p2 = new DirectPosition2D(crs104903, x2, y2);

        GeodeticCalculator gc = new GeodeticCalculator();
        gc.setStartingGeographicPoint(p1);
        gc.setDestinationGeographicPoint(p2);

        return gc.getAzimuth();
    }

    /**
     * èŽ·å–角度2
     */
    public static double getAngle2(double x1, double y1, double x2, double y2) {
        DirectPosition2D p1 = new DirectPosition2D(x1, y1);
        DirectPosition2D p2 = new DirectPosition2D(x2, y2);

        GeodeticCalculator gc = new GeodeticCalculator();
        gc.setStartingGeographicPoint(p1);
        gc.setDestinationGeographicPoint(p2);

        return gc.getAzimuth();
    }

    /**
     * æ ¹æ®è·ç¦»å’Œè§’度获取目标点
     */
    private static Point2D getPointByDistanceAndAngle(double x, double y, double angle, double distance) {
        DirectPosition2D p1 = new DirectPosition2D(x, y);

        GeodeticCalculator gc = new GeodeticCalculator();
        gc.setStartingGeographicPoint(p1);
        gc.setDirection(angle, distance);

        return gc.getDestinationGeographicPoint();
    }
}

 src/main/java/com/moon/server/service/data/RasterAnalysisService.java

@@ -4,6 +4,7 @@
import com.moon.server.entity.data.AnalysisResultEntity;
import com.moon.server.entity.data.MetaEntity;
import com.moon.server.entity.data.PublishEntity;
import com.moon.server.helper.GeoHelper;
import com.moon.server.helper.PathHelper;
import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
@@ -13,7 +14,10 @@
import org.gdal.gdal.gdal;
import org.gdal.ogr.Geometry;
import org.gdal.ogr.ogr;
import org.gdal.osr.SpatialReference;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.springframework.stereotype.Service;
import org.geotools.referencing.CRS;

import javax.annotation.Resource;
import java.io.File;
@@ -143,6 +147,58 @@
     * åˆ†æžçº¿
     */
    public void analysisPolyline(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) {
        /*Double len = geo.Length();


        double[] transform = ds.GetGeoTransform();
        // double rotationX = transform[2]; double rotationY = transform[4]
        double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = Math.abs(transform[5]);

        double buffer = Math.max(pixelWidth, pixelHeight) * size;
        double[] env = new double[4];
        geo = geo.Buffer(buffer);
        geo.GetEnvelope(env);

        int xMinPixel = (int) Math.floor((env[0] - minX) / pixelWidth);
        int yMinPixel = (int) Math.floor((maxY - env[3]) / pixelHeight);
        int xMaxPixel = (int) Math.floor((env[1] - minX) / pixelWidth);
        int yMaxPixel = (int) Math.floor((maxY - env[2]) / pixelHeight);

        int bandCount = ds.getRasterCount();
        int xMin = Math.min(xMinPixel, xMaxPixel);
        int xMax = Math.max(xMinPixel, xMaxPixel);
        int yMin = Math.min(yMinPixel, yMaxPixel);
        int yMax = Math.max(yMinPixel, yMaxPixel);

        List<Double> list = new ArrayList<>();
        for (int y = yMin; y <= yMax; y++) {
            for (int x = xMin; x <= xMax; x++) {
                Geometry point = new Geometry(ogr.wkbPoint);
                point.AddPoint(minX + pixelWidth * x, maxY - pixelHeight * y);
                if (geo.Intersects(point)) {
                    for (int i = 1; i <= bandCount; i++) {
                        double[] values = new double[1];
                        ds.GetRasterBand(i).ReadRaster(x, y, 1, 1, values);
                        list.add(values[0]);
                    }
                }
            }
        }

        processResult(list, entity);*/

        //double dis = GeoHelper.getDistance1()

        double lineDis = getPolylineDistance(geo);
        double distance = lineDis / Math.max(geo.GetPointCount(), size);

        List<double[]> list = getPointsByDistance(geo, distance);
    }

    /**
     * åˆ†æžçº¿
     */
    public void analysisPolyline2(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) {
        double[] transform = ds.GetGeoTransform();
        // double rotationX = transform[2]; double rotationY = transform[4]
        double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = Math.abs(transform[5]);
@@ -292,32 +348,34 @@
        return gdal.Warp("", new Dataset[]{dataset}, new WarpOptions(warpOptions));
    }

    public void processClippedDataByPolygon(String imagePath, String geometryString) {
        // æ³¨å†ŒGDAL驱动
        gdal.AllRegister();
        // æ‰“开栅格图像数据集
        Dataset dataset = gdal.Open(imagePath);
        if (dataset == null) {
            throw new RuntimeException("Failed to open raster dataset.");
    /**
     * èŽ·å–折线距离
     */
    private double getPolylineDistance(Geometry geo) {
        double dis = 0;
        for (int i = 0, c = geo.GetPointCount() - 1; i < c; i++) {
            double[] xy1 = geo.GetPoint(i);
            double[] xy2 = geo.GetPoint(i + 1);

            dis += GeoHelper.getDistance2(xy1[0], xy1[1], xy2[0], xy2[1]);
        }

        Geometry geometry = Geometry.CreateFromWkt(geometryString);
        return dis;
    }

        Dataset clippedDataset = clipRaster(dataset, geometry);
    /**
     * æ ¹æ®è·ç¦»èŽ·å–节点
     */
    private List<double[]> getPointsByDistance(Geometry geo, double distance) {
        List<double[]> list = new ArrayList<>();
        for (int i = 0, c = geo.GetPointCount() - 1; i < c; i++) {
            double[] xy1 = geo.GetPoint(i);
            double[] xy2 = geo.GetPoint(i + 1);

        int width = clippedDataset.GetRasterXSize();
        int height = clippedDataset.GetRasterYSize();
        int bandCount = clippedDataset.getRasterCount();
            double dis = GeoHelper.getDistance2(xy1[0], xy1[1], xy2[0], xy2[1]);

        double[] bandSum = new double[bandCount];
        for (int bandIndex = 1; bandIndex <= bandCount; bandIndex++) {
            Band band = clippedDataset.GetRasterBand(bandIndex);

            double[] pixelValues = new double[width * height];
            band.ReadRaster(0, 0, width, height, pixelValues);

            // å¤šæ³¢æ®µå½’一化处理
            bandSum[bandIndex - 1] = pixelValues[0];
        }

        return list;
    }
}