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.ArrayList;
|
import java.util.List;
|
import java.util.Vector;
|
|
/**
|
* 栅格分析服务
|
* @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> analysisPoint(Geometry point, int 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());
|
|
processPoint(entity, pub, point, size);
|
rs.add(entity);
|
}
|
|
return rs;
|
}
|
|
/**
|
* 设置错误
|
*/
|
private void setError(AnalysisResultEntity entity, String info) {
|
entity.setCode(StaticData.I500);
|
entity.setInfo(info);
|
}
|
|
/**
|
* 处理点
|
*/
|
private void processPoint(AnalysisResultEntity entity, PublishEntity pub, Geometry point, int size) {
|
List<MetaEntity> metas = publishService.selectMetasByPubid(pub.getId());
|
if (null == metas || metas.isEmpty()) {
|
setError(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, "源数据不存在");
|
return;
|
}
|
|
processPoint(entity, filePath, point, size);
|
}
|
|
/**
|
* 处理点
|
*/
|
public void processPoint(AnalysisResultEntity entity, String filePath, Geometry point, int size) {
|
Dataset ds = null;
|
try {
|
ds = gdal.Open(filePath);
|
if (null == ds) {
|
throw new Exception("打开栅格数据失败");
|
}
|
|
double x = point.GetX(), y = point.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 xPixel = (int) Math.floor((x - minX) / pixelWidth);
|
int yPixel = (int) Math.floor((maxY - y) / Math.abs(pixelHeight));
|
|
int bandCount = ds.getRasterCount();
|
|
List<double[]> sum = new ArrayList<>();
|
for (int bandIndex = 1; bandIndex <= bandCount; bandIndex++) {
|
double[] pixelValues = new double[size * size];
|
Band band = ds.GetRasterBand(bandIndex);
|
band.ReadRaster(xPixel, yPixel, size, size, pixelValues);
|
|
sum.add(pixelValues);
|
}
|
|
//return average(calculateAverage(sum));
|
} catch (Exception ex) {
|
setError(entity, ex.getMessage());
|
log.error(ex.getMessage(), ex);
|
} finally {
|
// gdal.GDALDestroyDriverManager()
|
if (null != ds) {
|
ds.delete();
|
}
|
}
|
}
|
|
/**
|
* 分析线
|
*/
|
public List<AnalysisResultEntity> analysisPolyline(Geometry polyline) {
|
List<AnalysisResultEntity> rs = new ArrayList<>();
|
//
|
|
return rs;
|
}
|
|
/**
|
* 分析面
|
*/
|
public List<AnalysisResultEntity> analysisPolygon(Geometry polygon) {
|
List<AnalysisResultEntity> rs = new ArrayList<>();
|
//
|
|
return rs;
|
}
|
|
|
public List<Double> processLine(String imagePath, String geometryString) {
|
// 打开栅格图像数据集
|
Dataset dataset = gdal.Open(imagePath);
|
if (dataset == null) {
|
throw new RuntimeException("Failed to open raster dataset.");
|
}
|
|
Geometry geometry = Geometry.CreateFromWkt(geometryString);
|
|
// 栅格像素范围
|
double[] geoTransform = dataset.GetGeoTransform();
|
double minX = geoTransform[0];
|
double pixelWidth = geoTransform[1];
|
double rotationX = geoTransform[2];
|
double maxY = geoTransform[3];
|
double rotationY = geoTransform[4];
|
double pixelHeight = geoTransform[5];
|
|
// geometry像素范围
|
double[] bounds = new double[6];
|
geometry.GetEnvelope(bounds);
|
double startX = bounds[0];
|
double endX = bounds[1];
|
double startY = bounds[2];
|
double endY = bounds[3];
|
|
int xStartPixel = (int) Math.floor((startX - minX) / pixelWidth);
|
int ySstartPixel = (int) Math.floor((maxY - startY) / Math.abs(pixelHeight));
|
int xEndPixel = (int) Math.floor((endX - minX) / pixelWidth);
|
int yEndPixel = (int) Math.floor((maxY - endY) / Math.abs(pixelHeight));
|
|
List<Double> values = new ArrayList<>();
|
for (int y = Math.min(ySstartPixel, yEndPixel); y <= Math.max(ySstartPixel, yEndPixel); y++) {
|
for (int x = Math.min(xStartPixel, xEndPixel); x <= Math.max(xStartPixel, xEndPixel); x++) {
|
Geometry point = new Geometry(ogr.wkbPoint);
|
point.AddPoint(minX + pixelWidth * x, maxY - pixelHeight * y);
|
if (geometry.Contains(point)) {
|
values.add(bandValueHandle(dataset, x, y));
|
}
|
}
|
}
|
|
return values;
|
}
|
|
public double processPolygon(String imagePath, String geometryString) {
|
// 打开栅格图像数据集
|
Dataset dataset = gdal.Open(imagePath);
|
|
if (dataset == null) {
|
throw new RuntimeException("Failed to open raster dataset.");
|
}
|
|
Geometry geometry = Geometry.CreateFromWkt(geometryString);
|
|
// 栅格像素范围
|
double[] geoTransform = dataset.GetGeoTransform();
|
double minX = geoTransform[0];
|
double pixelWidth = geoTransform[1];
|
double rotationX = geoTransform[2];
|
double maxY = geoTransform[3];
|
double rotationY = geoTransform[4];
|
double pixelHeight = geoTransform[5];
|
|
// geometry像素范围
|
double[] env = new double[6];
|
geometry.GetEnvelope(env);
|
int xMinPixel = (int) Math.floor((env[0] - minX) / pixelWidth);
|
int yMinPixel = (int) Math.floor((maxY - env[3]) / Math.abs(pixelHeight));
|
int xMaxPixel = (int) Math.ceil((env[1] - minX) / pixelWidth);
|
int yMaxPixel = (int) Math.ceil((maxY - env[2]) / Math.abs(pixelHeight));
|
|
int bandCount = dataset.getRasterCount();
|
int geometryWidth = Math.abs(xMaxPixel - xMinPixel);
|
int geometryHeight = Math.abs(yMaxPixel - yMinPixel);
|
|
List<double[]> sum = new ArrayList<>();
|
for (int bandIndex = 1; bandIndex <= bandCount; bandIndex++) {
|
Band band = dataset.GetRasterBand(bandIndex);
|
|
double[] pixelValues = new double[geometryWidth * geometryHeight];
|
band.ReadRaster(xMinPixel, yMinPixel, geometryWidth, geometryHeight, pixelValues);
|
|
// 多波段归一化处理
|
sum.add(pixelValues);
|
}
|
|
return average(calculateAverage(sum));
|
}
|
|
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 double[] 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);
|
}
|
|
return calculateAverage(sum);
|
}
|
|
public double 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] = average(pixelValues);
|
}
|
return average(bandSum);
|
}
|
|
public static double average(double[] values) {
|
double sum = 0;
|
for (double value : values) {
|
sum += value;
|
}
|
return sum / values.length;
|
}
|
|
public static double[] calculateAverage(List<double[]> list) {
|
if (list == null || list.isEmpty()) {
|
return new double[0];
|
}
|
|
int arrayLength = list.get(0).length;
|
double[] outArray = new double[arrayLength];
|
|
for (int i = 0; i < arrayLength; i++) {
|
double sum = 0;
|
|
for (double[] array : list) {
|
sum += array[i];
|
}
|
|
outArray[i] = sum / list.size();
|
}
|
|
return outArray;
|
}
|
|
public static double bandValueHandle(Dataset dataset, int x, int y) {
|
double[] bandSum = new double[dataset.getRasterCount()];
|
for (int bandIndex = 1; bandIndex <= dataset.getRasterCount(); bandIndex++) {
|
double[] pixelValues = new double[1];
|
Band band = dataset.GetRasterBand(bandIndex);
|
band.ReadRaster(x, y, 1, 1, pixelValues);
|
|
bandSum[bandIndex - 1] = pixelValues[0];
|
}
|
|
return average(bandSum);
|
}
|
}
|
