package com.moon.server.service.data;
|
|
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.PathHelper;
|
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.io.File;
|
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> 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;
|
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 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));
|
|
List<double[]> list = new ArrayList<>();
|
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);
|
|
list.add(pixelValues);
|
}
|
|
processResult(entity, list);
|
}
|
|
/**
|
* 分析线
|
*/
|
public void analysisPolyline(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);
|
}
|
|
/**
|
* 分析面
|
*/
|
public void analysisPolygon(AnalysisResultEntity entity, Dataset ds, Geometry geo) {
|
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[] env = new double[4];
|
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 geoWidth = Math.abs(xMaxPixel - xMinPixel);
|
int geoHeight = Math.abs(yMaxPixel - yMinPixel);
|
|
List<double[]> list = new ArrayList<>();
|
for (int i = 1; i <= bandCount; i++) {
|
Band band = ds.GetRasterBand(i);
|
double[] pixelValues = new double[geoWidth * geoHeight];
|
band.ReadRaster(xMinPixel, yMinPixel, geoWidth, geoHeight, pixelValues);
|
|
list.add(pixelValues);
|
}
|
|
processResult(entity, list);
|
}
|
|
/**
|
* 设置错误
|
*/
|
private void setError(AnalysisResultEntity entity, String info) {
|
entity.setCode(StaticData.I500);
|
entity.setInfo(info);
|
}
|
|
/**
|
* 处理结果
|
*/
|
private void processResult(AnalysisResultEntity entity, List<double[]> list) {
|
if (null == list || list.isEmpty()) {
|
return;
|
}
|
|
List<Double> rs = new ArrayList<>();
|
for (double[] ds : list) {
|
if (null != ds && ds.length > 0) {
|
for (double d : ds) {
|
rs.add(d);
|
}
|
}
|
}
|
|
processResult(rs, entity);
|
}
|
|
/**
|
* 处理结果
|
*/
|
private void processResult(List<Double> rs, AnalysisResultEntity entity) {
|
if (null == rs || rs.isEmpty()) {
|
return;
|
}
|
|
double avg = rs.stream().mapToDouble(Double::valueOf).average().getAsDouble();
|
entity.setMin(Collections.min(rs));
|
entity.setMax(Collections.max(rs));
|
entity.setAvg(avg);
|
}
|
|
public 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);
|
}
|
}
|
|
public static 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));
|
}
|
|
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.");
|
}
|
|
Geometry geometry = Geometry.CreateFromWkt(geometryString);
|
|
Dataset clippedDataset = clipRaster(dataset, geometry);
|
|
int width = clippedDataset.GetRasterXSize();
|
int height = clippedDataset.GetRasterYSize();
|
int bandCount = clippedDataset.getRasterCount();
|
|
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];
|
}
|
}
|
}
|
