断面分析功能实现

dcb

2025-07-09 e53252b99e7b49b435b7a6ee3eab21ae1bd7a055

断面分析功能实现

 src/main/java/com/se/nsl/config/AsyncExecutor.java

@@ -14,7 +14,7 @@
        executor.setCorePoolSize(5);
        executor.setMaxPoolSize(10);
        executor.setQueueCapacity(25);
        executor.setThreadNamePrefix("RealTimeExecutor-");
        executor.setThreadNamePrefix("RTExecutor-");
        executor.initialize();
        return executor;
    }

 src/main/java/com/se/nsl/controller/SimuController.java

@@ -10,6 +10,7 @@
import com.se.nsl.utils.SimulateType;
import io.swagger.annotations.Api;
import io.swagger.annotations.ApiOperation;
import io.swagger.v3.oas.annotations.Parameter;
import lombok.extern.slf4j.Slf4j;
import org.gdal.ogr.Geometry;
import org.gdal.ogr.ogr;
@@ -205,24 +206,40 @@

    @ApiOperation(value = "crossSection")
    @GetMapping("/crossSection")
    public R<Object> crossSection(String serviceName, double[] startPoint, double[] endPoint) {
    public R<Object> crossSection(@Parameter(description = "方案id，示例：50") Integer id,
                                  @Parameter(description = "时间戳，示例：1751552400000") Long time,
                                  @Parameter(description = "起点坐标，示例：116.59049537485063,40.564178548127686") String startPoint,
                                  @Parameter(description = "终点坐标，示例：116.5901406492509,40.56499045715429") String endPoint) {
        if (null == id || id < 1) return clientError("id不能为空");
        Simu simu = simuService.selectById(id);
        if (simu == null) {
            return clientError("找不到对应的服务");
        }
        String serviceName = simu.getServiceName();
        if (serviceName == null) {
            return clientError("服务名不能为空");
            return fail("找不到对应的服务");
        }
        if (startPoint == null) {
            return clientError("起点不能为空");
        }
        if (startPoint.length < 2) {
            return clientError("起点至少包含x,y两个值");
        }
        if (endPoint == null) {
            return clientError("终点不能为空");
        }
        if (endPoint.length < 2) {
            return clientError("终点至少包含x,y两个值");
        String[] sp = startPoint.split(",");
        double[] s = new double[] {Double.parseDouble(sp[0]), Double.parseDouble(sp[1])};
        String[] ep = endPoint.split(",");
        double[] e = new double[] {Double.parseDouble(ep[0]), Double.parseDouble(ep[1])};
        try {
            if (time == null) {
                List<CrossSectionAnalysisResult> result = scas.crossSectionAnalysis(serviceName, s, e);
                return success(result, result.size());
            } else {
                CrossSectionAnalysisResult res = scas.crossSectionAnalysis(serviceName, s, e, time);
                return success(res, 1);
            }
        } catch (IllegalArgumentException ex) {
            return fail(ex.getMessage(), null);
        }
        List<CrossSectionAnalysisResult> result = scas.crossSectionAnalysis(serviceName, startPoint, endPoint);
        return success(result);
    }

    @ApiOperation(value = "stop")

 src/main/java/com/se/nsl/domain/vo/CrossSectionAnalysisResult.java

¶Ô±ÈÐÂÎļþ
@@ -0,0 +1,49 @@
package com.se.nsl.domain.vo;

public class CrossSectionAnalysisResult {
    private double depth;
    private double velocity; //流速
    private double flowRate; //流量
    private long time; //时间

    public CrossSectionAnalysisResult() {}

    public CrossSectionAnalysisResult(double depth,double flowRate, double velocity, long time) {
        this.depth = depth;
        this.flowRate = flowRate;
        this.velocity = velocity;
        this.time = time;
    }

    public double getDepth() {
        return depth;
    }

    public void setDepth(double depth) {
        this.depth = depth;
    }

    public double getFlowRate() {
        return flowRate;
    }

    public void setFlowRate(double flowRate) {
        this.flowRate = flowRate;
    }

    public double getVelocity() {
        return velocity;
    }

    public void setVelocity(double velocity) {
        this.velocity = velocity;
    }

    public long getTime() {
        return time;
    }

    public void setTime(long time) {
        this.time = time;
    }
}

 src/main/java/com/se/nsl/service/CrossSectionAnalysisService.java

¶Ô±ÈÐÂÎļþ
@@ -0,0 +1,317 @@
package com.se.nsl.service;

import com.se.nsl.config.PropertiesConfig;
import com.se.nsl.domain.vo.CrossSectionAnalysisResult;
import com.se.nsl.domain.vo.SimuResult;
import com.se.nsl.utils.CoordinateTransformer;
import com.se.nsl.utils.SolverTifUtil;
import com.se.nsl.utils.TimeFormatUtil;
import lombok.extern.slf4j.Slf4j;
import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;
import org.locationtech.jts.geom.*;
import org.springframework.stereotype.Service;

import javax.annotation.Resource;
import java.io.File;
import java.nio.file.Paths;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;

@Service
@Slf4j
public class CrossSectionAnalysisService {

    public static final String TIF = ".tif";
    public static final String YYYY_MM_DD_HH_MM_SS = "yyyyMMddHHmmss";
    private static final double SAMPLING_DISTANCE = 1;
    @Resource
    PropertiesConfig config;

    public List<CrossSectionAnalysisResult> crossSectionAnalysis(String serviceName, double[] startPoint, double[] endPoint) {
        List<CrossSectionAnalysisResult> results = new ArrayList<>();
        //找到指定的服务目录的地形文件
        String inPath = config.getInPath();
        File serviceNameDir = new File(inPath, serviceName);
        if (!serviceNameDir.exists()) {
            throw new IllegalArgumentException(serviceName + "不存在");
        }
        //根据起点和终点画条线，与地形构建一个多边形
        File dem = new File(serviceNameDir, "DEM.tif");
        ChannelCrossSectionExtractor extractor = new ChannelCrossSectionExtractor(dem.getAbsolutePath());
        Point s = coordinateTransform(startPoint);
        Point e = coordinateTransform(endPoint);
        extractor.extractCrossSection(s, e, SAMPLING_DISTANCE);
        File depthDir = new File(serviceNameDir, "depth");
        File[] files = depthDir.listFiles();
        if (files == null) return Collections.emptyList();
        for (File tif : files) {
            String name = tif.getName();
            if (!name.endsWith(TIF)) continue;
            CrossSectionAnalysisResult ar = getCrossSectionAnalysisResult(tif, extractor);
            if (ar != null) {
                results.add(ar);
            }
        }
        extractor.destroy();
        return results;
    }

    private CrossSectionAnalysisResult getCrossSectionAnalysisResult(File tif, ChannelCrossSectionExtractor extractor) {
        Point position = extractor.getMaxDepthPosition();
        if (position == null) return null;
        SimuResult sr = querySimuResult(tif, position);
        //查询对应点的水深
        double waterDepth = sr.getDepth();
        Polygon polygon = extractor.generateWettedPolygon(waterDepth);
        double area = calculateArea(polygon);
        //查询对应点的流速
        double v = sr.getVelocity();
        //计算出此时的水流量
        double flowRate = v * area;
        String name = tif.getName();
        String prefix = name.replace(TIF, "");
        long time = TimeFormatUtil.toMillis(prefix, YYYY_MM_DD_HH_MM_SS);
        return new CrossSectionAnalysisResult(waterDepth, flowRate, v, time);
    }

    public CrossSectionAnalysisResult crossSectionAnalysis(String serviceName, double[] startPoint, double[] endPoint, long time) {
        //找到指定的服务目录的地形文件
        String inPath = config.getInPath();
        File serviceNameDir = new File(inPath, serviceName);
        if (!serviceNameDir.exists()) {
            throw new IllegalArgumentException(serviceName + "不存在");
        }
        //根据起点和终点画条线，与地形构建一个多边形
        File dem = new File(serviceNameDir, "DEM.tif");
        ChannelCrossSectionExtractor extractor = new ChannelCrossSectionExtractor(dem.getAbsolutePath());
        Point s = coordinateTransform(startPoint);
        Point e = coordinateTransform(endPoint);
        extractor.extractCrossSection(s, e, SAMPLING_DISTANCE);
        String prefix = TimeFormatUtil.formatTime(time, YYYY_MM_DD_HH_MM_SS);
        File tif = Paths.get(inPath, serviceName, "depth", prefix + ".tif").toFile();
        CrossSectionAnalysisResult ar = getCrossSectionAnalysisResult(tif, extractor);
        extractor.destroy();
        return ar;
    }

    private Point coordinateTransform(double[] startPoint) {
        double lon = startPoint[0];
        double lat = startPoint[1];
        double[] xy = CoordinateTransformer.transform(4326, config.getEpsg(), lon, lat);
        return new Point(xy[0], xy[1]);
    }

    private SimuResult querySimuResult(File tif, Point pos) {
        double[] xy = new double[]{pos.x, pos.y};
        return SolverTifUtil.getSimuResult(tif, xy);
    }


    private double calculateArea(Geometry geometry) {
        if (geometry == null) return 0;
        return geometry.getArea();
    }


    private static class ChannelCrossSectionExtractor {
        private final GeometryFactory geometryFactory;
        private final Dataset dataset;
        private final List<CrossSectionPoint> crossSection;
        private Point maxDepthPosition; //最深位置的坐标
        private double minElevation = Double.MAX_VALUE;
        private double maxElevation = Double.MIN_VALUE;

        public ChannelCrossSectionExtractor(String tifPath) {
            geometryFactory = new GeometryFactory();
            // è¯»å–TIFF文件
            File file = new File(tifPath);
            this.dataset = gdal.Open(file.getAbsolutePath(), gdalconstConstants.GA_ReadOnly);
            crossSection = new ArrayList<>();
        }

        public void destroy() {
            if (dataset != null) {
                dataset.delete();
            }
        }

        /**
         * æå–两点之间的地形截面
         */
        public void extractCrossSection(Point startPoint, Point endPoint, double samplingDistance) {
            // è®¡ç®—两点之间的距离
            double distance = Math.sqrt(Math.pow(endPoint.x - startPoint.x, 2) + Math.pow(endPoint.y - startPoint.y, 2));
            // è®¡ç®—采样点数量
            int numPoints = (int) Math.ceil(distance / samplingDistance) + 1;

            for (int i = 0; i < numPoints; i++) {
                // æ²¿ç›´çº¿æ’值计算采样点坐标
                double fraction = (double) i / (numPoints - 1);
                double x = startPoint.x + fraction * (endPoint.x - startPoint.x);
                double y = startPoint.y + fraction * (endPoint.y - startPoint.y);
                // è®¡ç®—距离起点的水平距离
                double horizontalDistance = fraction * distance;
                // èŽ·å–该点的高程值
                double elevation = getElevation(x, y);
                minElevation = Math.min(minElevation, elevation);
                maxElevation = Math.max(maxElevation, elevation);
                crossSection.add(new CrossSectionPoint(horizontalDistance, x, y, elevation));
            }
        }

        /**
         * èŽ·å–指定位置的高程值
         */
        private double getElevation(double x, double y) {
            double[] geoTransform = dataset.GetGeoTransform();
            //计算栅格行列号
            int col = (int) ((x - geoTransform[0]) / geoTransform[1]);
            int row = (int) ((geoTransform[3] - y) / Math.abs(geoTransform[5]));
            Band band = dataset.GetRasterBand(1);
            float[] values = new float[1];
            band.ReadRaster(col, row, 1, 1, values);
            return values[0];
        }

        /**
         * ç”Ÿæˆè¿‡æ°´æ–­é¢å¤šè¾¹å½¢
         */
        private Polygon generateWettedPolygon(double waterLevel) {
            if (crossSection.isEmpty()) {
                return null;
            }
            // åˆ›å»ºåœ°å½¢æˆªé¢çº¿
            Coordinate[] terrainCoords = new Coordinate[crossSection.size()];
            for (int i = 0; i < crossSection.size(); i++) {
                CrossSectionPoint point = crossSection.get(i);
                double elevation = point.getElevation();
                terrainCoords[i] = new Coordinate(point.getDistance(), elevation);
//                log.info("{},{}", point.getDistance(), point.getElevation());
            }
            LineString terrainLine = geometryFactory.createLineString(terrainCoords);
//            log.info(terrainLine.toString());
            // åˆ›å»ºæ°´ä½çº¿
            CrossSectionPoint first = crossSection.get(0);
            CrossSectionPoint last = crossSection.get(crossSection.size() - 1);
            double minDistance = first.getDistance();
            double maxDistance = last.getDistance();
            Coordinate[] waterCoords = new Coordinate[2];
            if (maxElevation - minElevation > waterLevel) {
                waterLevel = minElevation + waterLevel;
            } else {
                waterLevel = Math.min(first.getElevation(), last.getElevation());
            }
//            log.info("waterLevle:{}", waterLevel);
            waterCoords[0] = new Coordinate(minDistance, waterLevel);
            waterCoords[1] = new Coordinate(maxDistance, waterLevel);
            LineString waterLine = geometryFactory.createLineString(waterCoords);
//            log.info(waterLine.toString());
            // è®¡ç®—水位线与地形的交点
            Geometry intersections = terrainLine.intersection(waterLine);
            if (intersections.isEmpty()) {
                return null; // æ°´ä½ä½ŽäºŽåœ°å½¢ï¼Œæ— æ°´
            }
            // èŽ·å–最外侧的两个交点作为河岸
            Coordinate leftBank;
            Coordinate rightBank;
            if (intersections instanceof org.locationtech.jts.geom.Point) {
                return null; // ç‰¹æ®Šæƒ…况：水位刚好与地形相切于一点
            } else if (intersections instanceof MultiPoint) {
                MultiPoint multiPoint = (MultiPoint) intersections;
                int num = multiPoint.getNumGeometries();
                if (num >= 2) {
                    leftBank = multiPoint.getGeometryN(0).getCoordinate();
                    rightBank = multiPoint.getGeometryN(num - 1).getCoordinate();
                } else {
                    return null;
                }
            } else {
                return null;
            }

            // æå–水面以下的地形点
            double finalWaterLevel = waterLevel;
            List<Coordinate> wettedCoords = crossSection.stream()
                    .filter(c -> c.getElevation() < finalWaterLevel)
                    .map(s -> new Coordinate(s.getDistance(), s.getElevation()))
                    .collect(Collectors.toList());
            wettedCoords.add(0, leftBank);
            wettedCoords.add(rightBank);
            // é—­åˆå¤šè¾¹å½¢
            if (wettedCoords.size() >= 3) {
                // ç¡®ä¿å¤šè¾¹å½¢é—­åˆï¼ˆé¦–尾相连）
                if (!wettedCoords.get(0).equals2D(wettedCoords.get(wettedCoords.size() - 1))) {
                    wettedCoords.add(new Coordinate(wettedCoords.get(0)));
                }
                return geometryFactory.createPolygon(wettedCoords.toArray(new Coordinate[0]));
            } else {
                return null; // æ— æ³•å½¢æˆæœ‰æ•ˆå¤šè¾¹å½¢
            }
        }

        public Point getMaxDepthPosition() {
            if (maxDepthPosition != null) {
                return maxDepthPosition;
            } else {
                for (CrossSectionPoint csp : crossSection) {
                    double elevation = csp.getElevation();
                    if (Math.abs(elevation - minElevation) < 1e-15) {
                        maxDepthPosition = new Point(csp.getX(), csp.getY());
                        return maxDepthPosition;
                    }
                }
            }
            return null;
        }
    }

    /**
     * æˆªé¢ç‚¹ç±»
     */
    private static class CrossSectionPoint {
        private final double distance; // è·ç¦»èµ·ç‚¹çš„水平距离
        private final double x;
        private final double y;
        private final double elevation; // é«˜ç¨‹

        public CrossSectionPoint(double distance, double x, double y, double elevation) {
            this.distance = distance;
            this.x = x;
            this.y = y;
            this.elevation = elevation;
        }

        public double getDistance() {
            return distance;
        }

        public double getX() {
            return x;
        }

        public double getY() {
            return y;
        }

        public double getElevation() {
            return elevation;
        }
    }

    private static class Point {
        public double x;
        public double y;

        public Point(double x, double y) {
            this.x = x;
            this.y = y;
        }

    }

}

 src/main/java/com/se/nsl/utils/SolverTifUtil.java

@@ -60,7 +60,7 @@
        return readPixelValue(cr.dataset, cr.col, cr.row, band);
    }

    private static ColumnRow getColumnRow(File tifFile, double x, double y) {
    public static ColumnRow getColumnRow(File tifFile, double x, double y) {
        Dataset dataset = gdal.Open(tifFile.getAbsolutePath(), gdalconstConstants.GA_ReadOnly);
        // èŽ·å–地理变换参数（6元素数组）
        // [0]: å·¦ä¸Šè§’X坐标, [1]: åƒå…ƒå®½åº¦, [2]: X方向旋转,

 src/main/resources/application-dev.yml

@@ -149,8 +149,8 @@
  rainfallSite: beijing
  epsg: 4548
#  saveFrames: 3
  #生成帧数的间隔时间，单位是分钟，设置为5表示每隔5分钟生成一帧
  saveFrameInterval: 20
  #生成帧数的间隔时间，单位是秒，设置为30表示每隔30秒生成一帧
  saveFrameInterval: 600
  saveFilter: 0.015
  evaporation: 0.27
  # åœŸåœ°åˆ©ç”¨ï¼š1-Cropland,2-Forest,3-Shrub,4-Grassland,5-Water,6-Snow/Ice,7-Barren,8-Impervious,9-Wetland
@@ -172,9 +172,10 @@
  #tif中的高程文件名要含有dem，土地利用要含有landuse，站点要含有station，不限制大小写，文字尽量用英文字母表达
  keyDitch: D:\other\simu\CudaUWSolver-2.2.1\KeyDitch
realtime-simulate-config:
  url: http://192.168.56.106:9522/ylclyPacket/getData
  token: YjhhYjAwOWFhMjk1MTM1ZDA0NGU3ZWZlMDQzMzUzZDE1MGJmY2Q4ZWEyYjliNjQzZjcwMjhlNDY0ZjAxNWZjOTZmNzMwYmNmZDA2YmVmNTIzNjU0ZDgzODRjYTUxYTM1
  realTimeInterval: 5
  url: http://192.168.56.106:9533/ylclyPacket/getData
  token: NmNmNWZlMThmNTExYzEyMjEwMzNlMGViYjFiN2Y2MmEwMzkwNjA2MmFlZjYzMzU0MzE4YTIzMDUyZTM2MzU5ZWJjMDllNTYwZDNiN2JjZTU5YTFhNjg5Y2IwZGRlODRh
  #实时模拟每帧的间隔，单位秒
  frameInterval: 5
  #请求雨量计数据时，时间范围相差多少分钟
  requestOffsetMinutes: 15
  #往前偏移多少天，为了方便实时模拟以前的降雨情况,设置为0表示，从当前时间开始计算