package com.moon.server.service.data;
|
|
import com.alibaba.fastjson.JSONObject;
|
import com.moon.server.entity.all.StaticData;
|
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 com.moon.server.helper.WebHelper;
|
import com.moon.server.service.all.WebSocketService;
|
import org.apache.commons.logging.Log;
|
import org.apache.commons.logging.LogFactory;
|
import org.gdal.gdal.Band;
|
import org.gdal.gdal.Dataset;
|
import org.gdal.gdal.WarpOptions;
|
import org.gdal.gdal.gdal;
|
import org.gdal.ogr.Geometry;
|
import org.gdal.ogr.ogr;
|
import org.springframework.stereotype.Service;
|
|
import javax.annotation.Resource;
|
import java.awt.geom.Point2D;
|
import java.io.File;
|
import java.io.IOException;
|
import java.util.*;
|
|
/**
|
* 栅格分析服务
|
* @author LJW
|
* @date 2023-09-01 6:16
|
*/
|
@Service
|
public class RasterAnalysisService {
|
@Resource
|
PathHelper pathHelper;
|
|
@Resource
|
PublishService publishService;
|
|
private final static Log log = LogFactory.getLog(RasterAnalysisService.class);
|
|
/**
|
* 测试
|
*/
|
public List<AnalysisResultEntity> test(Geometry geo, Integer size) {
|
List<AnalysisResultEntity> rs = new ArrayList<>();
|
AnalysisResultEntity entity = new AnalysisResultEntity();
|
entity.setLayerName("Test");
|
|
openRaster(entity, "D:\\Moon\\data\\DOM\\Lunar_LRO_LOLA_ClrShade_Global_128ppd_v04_2.tif", geo, size);
|
rs.add(entity);
|
|
return rs;
|
}
|
|
/**
|
* 使用WKT查询分析，结果以消息推送
|
*/
|
public void analysisForPost(Geometry geo, Integer size, String token) {
|
List<PublishEntity> pubs = publishService.selectRaster();
|
if (null == pubs || pubs.isEmpty()) {
|
return;
|
}
|
|
// for (PublishEntity pub : pubs) {
|
pubs.parallelStream().forEach(pub -> {
|
try {
|
AnalysisResultEntity entity = new AnalysisResultEntity(token);
|
entity.setLayerName(pub.getName());
|
|
List<MetaEntity> metas = publishService.selectMetasByPubid(pub.getId());
|
if (null == metas || metas.isEmpty()) {
|
setError(entity, "找不到发布数据");
|
postInfo(entity);
|
return;
|
}
|
|
String filePath = pathHelper.getConfig().getUploadPath() + File.separator + metas.get(0).getPath();
|
File file = new File(filePath);
|
if (!file.exists() || file.isDirectory()) {
|
setError(entity, "源数据不存在");
|
postInfo(entity);
|
return;
|
}
|
|
openRaster(entity, filePath, geo, size);
|
postInfo(entity);
|
} catch (Exception ex) {
|
log.error(ex.getMessage(), ex);
|
}
|
});
|
}
|
|
/**
|
* 推送消息
|
*/
|
private void postInfo(AnalysisResultEntity entity) throws IOException {
|
JSONObject map = new JSONObject();
|
map.put("analysisForPost", entity);
|
|
String json = JSONObject.toJSONString(map);
|
// System.out.println(json)
|
WebSocketService.broadCastInfo(json);
|
}
|
|
/**
|
* 分析方法
|
*/
|
public List<AnalysisResultEntity> analysis(Geometry geo, Integer size) {
|
List<AnalysisResultEntity> rs = new ArrayList<>();
|
|
List<PublishEntity> pubs = publishService.selectRaster();
|
if (null == pubs || pubs.isEmpty()) {
|
return rs;
|
}
|
|
for (PublishEntity pub : pubs) {
|
AnalysisResultEntity entity = new AnalysisResultEntity();
|
entity.setLayerName(pub.getName());
|
|
List<MetaEntity> metas = publishService.selectMetasByPubid(pub.getId());
|
if (null == metas || metas.isEmpty()) {
|
setError(entity, "找不到发布数据");
|
continue;
|
}
|
|
String filePath = pathHelper.getConfig().getUploadPath() + File.separator + metas.get(0).getPath();
|
File file = new File(filePath);
|
if (!file.exists() || file.isDirectory()) {
|
setError(entity, "源数据不存在");
|
continue;
|
}
|
|
openRaster(entity, filePath, geo, size);
|
rs.add(entity);
|
}
|
|
return rs;
|
}
|
|
/**
|
* 打开栅格数据
|
*/
|
private void openRaster(AnalysisResultEntity entity, String filePath, Geometry geo, int size) {
|
Dataset ds = null;
|
filePath = "D:\\Moon\\data\\DOM\\test.tif";
|
try {
|
ds = gdal.Open(filePath);
|
if (null == ds) {
|
throw new Exception("打开栅格数据失败");
|
}
|
|
switch (geo.GetGeometryType()) {
|
case ogr.wkbPoint:
|
analysisPoint(entity, ds, geo, size);
|
break;
|
case ogr.wkbLineString:
|
analysisPolyline(entity, ds, geo, size);
|
break;
|
default:
|
analysisPolygon(entity, ds, geo);
|
}
|
} catch (Exception ex) {
|
setError(entity, ex.getMessage());
|
log.error(ex.getMessage(), ex);
|
} finally {
|
// gdal.GDALDestroyDriverManager()
|
if (null != ds) {
|
ds.delete();
|
}
|
}
|
}
|
|
/**
|
* 分析点
|
*/
|
public void analysisPoint(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) {
|
double[] transform = ds.GetGeoTransform();
|
double pixelWidth = transform[1], pixelHeight = Math.abs(transform[5]);
|
|
double buffer = Math.max(pixelWidth, pixelHeight) * size;
|
geo = geo.Buffer(buffer);
|
|
analysisPolygon(entity, ds, geo);
|
}
|
|
/**
|
* 分析线
|
*/
|
public void analysisPolyline(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) {
|
double[] transform = ds.GetGeoTransform();
|
int xSize = ds.getRasterXSize(), ySize = ds.getRasterYSize();
|
double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = transform[5];
|
|
double len = 0;
|
int bandCount = ds.getRasterCount();
|
double lineDis = getPolylineDistance(geo);
|
double distance = lineDis / Math.max(geo.GetPointCount(), size);
|
|
List<double[]> points = getPointsByDistance(geo, distance);
|
for (int i = 0, c = points.size(); i < c; i++) {
|
double[] xy = points.get(i);
|
if (i > 0) {
|
double[] p = points.get(i - 1);
|
len += GeoHelper.getDistance(p[0], p[1], xy[0], xy[1]);
|
}
|
|
int xPixel = (int) Math.floor((xy[0] - minX) / pixelWidth);
|
int yPixel = (int) Math.floor((maxY - xy[1]) / Math.abs(pixelHeight));
|
if (xPixel < 1 || xPixel > xSize || yPixel < 1 || yPixel > ySize) {
|
continue;
|
}
|
|
List<Double> vals = new ArrayList<>();
|
for (int j = 1; j <= bandCount; j++) {
|
double[] pixelValues = new double[1];
|
ds.GetRasterBand(j).ReadRaster(xPixel, yPixel, 1, 1, pixelValues);
|
if (!Double.isNaN(pixelValues[0])) {
|
vals.add(WebHelper.round(pixelValues[0], 3));
|
}
|
}
|
if (vals.size() > 0) {
|
entity.addPoint(xy[0], xy[1], len, vals);
|
}
|
}
|
}
|
|
/**
|
* 分析面
|
*/
|
public void analysisPolygon(AnalysisResultEntity entity, Dataset ds, Geometry geo) {
|
double[] transform = ds.GetGeoTransform();
|
int xSize = ds.getRasterXSize(), ySize = ds.getRasterYSize();
|
double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = Math.abs(transform[5]), rotationX = transform[2], rotationY = transform[4];
|
double[] env = new double[4];
|
geo.GetEnvelope(env);
|
entity.addPoint(geo.Centroid().GetX(), geo.Centroid().GetY(), 0, null);
|
|
int xMinPixel = Math.max((int) Math.floor((env[0] - minX) / pixelWidth), 1);
|
int yMinPixel = Math.max((int) Math.floor((maxY - env[3]) / pixelHeight), 1);
|
int xMaxPixel = Math.min((int) Math.floor((env[1] - minX) / pixelWidth), xSize);
|
int yMaxPixel = Math.min((int) Math.floor((maxY - env[2]) / pixelHeight), ySize);
|
if (xMaxPixel < 1 || yMaxPixel < 1 || xMaxPixel - xMinPixel < 0 || yMaxPixel - yMinPixel < 0) {
|
setError(entity, "查询范围无效");
|
return;
|
}
|
|
int bandCount = ds.getRasterCount();
|
int width = xMaxPixel - xMinPixel;
|
int height = yMaxPixel - yMinPixel;
|
if (width * height > StaticData.I64 * StaticData.I64) {
|
readRasterForBlocks(entity, ds, bandCount, xMinPixel, yMinPixel, width, height);
|
return;
|
}
|
|
for (int i = 1; i <= bandCount; i++) {
|
double[] pixelValues = new double[width * height];
|
ds.GetRasterBand(i).ReadRaster(xMinPixel, yMinPixel, width, height, pixelValues);
|
setBandVals(entity, pixelValues);
|
}
|
}
|
|
/**
|
* 按照块读取栅格数据
|
*/
|
private void readRasterForBlocks(AnalysisResultEntity entity, Dataset ds, int bandCount, int xMinPixel, int yMinPixel, int width, int height) {
|
List<Integer> xList = getSamples(xMinPixel, width - 1);
|
List<Integer> yList = getSamples(yMinPixel, height - 1);
|
|
double[] pixelValues = new double[1];
|
for (int i = 1; i <= bandCount; i++) {
|
List<Double> list = new ArrayList<>();
|
for (Integer x : xList) {
|
for (Integer y : yList) {
|
ds.GetRasterBand(i).ReadRaster(x, y, 1, 1, pixelValues);
|
if (!Double.isNaN(pixelValues[0])) {
|
list.add(pixelValues[0]);
|
}
|
}
|
}
|
setBandVals(entity, list);
|
}
|
}
|
|
/**
|
* 获取抽样列表
|
*/
|
private List<Integer> getSamples(int start, int strip) {
|
List<Integer> list = new ArrayList<>();
|
|
double avg = 1.0 * strip / (StaticData.I64 - 1);
|
for (int i = 0; i < StaticData.I64; i++) {
|
list.add(start + (int) Math.ceil(avg * i));
|
}
|
|
return list;
|
}
|
|
/**
|
* 设置错误
|
*/
|
private void setError(AnalysisResultEntity entity, String info) {
|
entity.setCode(StaticData.I500);
|
entity.setInfo(info);
|
}
|
|
/**
|
* 设置Band值
|
*/
|
private void setBandVals(AnalysisResultEntity entity, double[] pixelValues) {
|
if (null == pixelValues || pixelValues.length == 0) {
|
return;
|
}
|
|
List<Double> list = new ArrayList<>();
|
for (double val : pixelValues) {
|
if (!Double.isNaN(val)) {
|
list.add(val);
|
}
|
}
|
setBandVals(entity, list);
|
}
|
|
/**
|
* 设置Band值
|
*/
|
private void setBandVals(AnalysisResultEntity entity, List<Double> list) {
|
if (null == list || list.isEmpty()) {
|
return;
|
}
|
|
double min = Collections.min(list);
|
double max = Collections.max(list);
|
double avg = list.stream().mapToDouble(Double::valueOf).average().getAsDouble();
|
entity.addBandVals(min, avg, max);
|
}
|
|
/**
|
* 获取折线距离
|
*/
|
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.getDistance(xy1[0], xy1[1], xy2[0], xy2[1]);
|
}
|
|
return dis;
|
}
|
|
/**
|
* 根据距离获取节点
|
*/
|
private List<double[]> getPointsByDistance(Geometry geo, double distance) {
|
List<double[]> list = new ArrayList<>();
|
|
int c = geo.GetPointCount() - 1;
|
for (int i = 0; i < c; i++) {
|
double[] xy1 = geo.GetPoint(i);
|
double[] xy2 = geo.GetPoint(i + 1);
|
list.add(xy1);
|
|
double lineDis = GeoHelper.getDistance(xy1[0], xy1[1], xy2[0], xy2[1]);
|
if (lineDis < distance) {
|
continue;
|
}
|
|
double d = distance;
|
double angle = GeoHelper.getAngle(xy1[0], xy1[1], xy2[0], xy2[1]);
|
while (d + distance * StaticData.D05 < lineDis) {
|
Point2D point = GeoHelper.getPointByDistanceAndAngle(xy1[0], xy1[1], angle, d);
|
list.add(new double[]{point.getX(), point.getY()});
|
d += distance;
|
}
|
}
|
list.add(geo.GetPoint(c));
|
|
return list;
|
}
|
|
/**
|
* 分析点 2
|
*/
|
public void analysisPoint2(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) {
|
double x = geo.GetX(), y = geo.GetY();
|
double[] transform = ds.GetGeoTransform();
|
// double rotationX = transform[2]; double rotationY = transform[4]
|
double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = transform[5];
|
|
int bandCount = ds.getRasterCount();
|
int xPixel = (int) Math.floor((x - minX) / pixelWidth);
|
int yPixel = (int) Math.floor((maxY - y) / Math.abs(pixelHeight));
|
|
for (int i = 1; i <= bandCount; i++) {
|
double[] pixelValues = new double[size * size];
|
Band band = ds.GetRasterBand(i);
|
band.ReadRaster(xPixel, yPixel, size, size, pixelValues);
|
}
|
}
|
|
/**
|
* 分析线 2
|
*/
|
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]);
|
|
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)
|
}
|
|
private void processClippedDataByLine(String imagePath, String geometryString) {
|
// 注册GDAL驱动
|
gdal.AllRegister();
|
// 打开栅格图像数据集
|
Dataset dataset = gdal.Open(imagePath);
|
if (dataset == null) {
|
throw new RuntimeException("Failed to open raster dataset.");
|
}
|
|
Geometry geometry = Geometry.CreateFromWkt(geometryString);
|
|
Dataset clippedDataset = clipRaster(dataset, geometry);
|
|
int width = clippedDataset.GetRasterXSize();
|
int height = clippedDataset.GetRasterYSize();
|
int bandCount = clippedDataset.getRasterCount();
|
|
List<double[]> sum = new ArrayList<>();
|
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);
|
|
sum.add(pixelValues);
|
}
|
}
|
|
private Dataset clipRaster(Dataset dataset, Geometry geometry) {
|
Vector<String> warpOptions = new Vector<>();
|
warpOptions.add("-crop_to_cutline");
|
warpOptions.add("-cutline");
|
warpOptions.add(geometry.ExportToWkt());
|
warpOptions.add("-dstalpha");
|
|
return gdal.Warp("", new Dataset[]{dataset}, new WarpOptions(warpOptions));
|
}
|
}
|
