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;
|
}
|
|
}
|
|
}
|
