package com.moon.server.service.data;
|
|
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 java.util.ArrayList;
|
import java.util.List;
|
import java.util.Vector;
|
|
/**
|
* 栅格分析服务
|
* @author LJW
|
* @date 2023-9-1 6:16
|
*/
|
@Service
|
public class RasterAnalysisService {
|
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 double processPoint(String imagePath, String geometryString, int size) {
|
// 打开栅格图像数据集
|
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];
|
|
double x = geometry.GetX();
|
double y = geometry.GetY();
|
|
int xPixel = (int) Math.floor((x - minX) / pixelWidth);
|
int yPixel = (int) Math.floor((maxY - y) / Math.abs(pixelHeight));
|
|
int bandCount = dataset.getRasterCount();
|
|
List<double[]> sum = new ArrayList<>();
|
for (int bandIndex = 1; bandIndex <= bandCount; bandIndex++) {
|
double[] pixelValues = new double[size * size];
|
Band band = dataset.GetRasterBand(bandIndex);
|
band.ReadRaster(xPixel, yPixel, size, size, pixelValues);
|
|
sum.add(pixelValues);
|
}
|
|
// 释放资源
|
gdal.GDALDestroyDriverManager();
|
|
return average(calculateAverage(sum));
|
}
|
|
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 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);
|
}
|
}
|
