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 analysisPoint(Geometry point, int size) { List rs = new ArrayList<>(); List 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 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 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 analysisPolyline(Geometry polyline) { List rs = new ArrayList<>(); // return rs; } /** * 分析面 */ public List analysisPolygon(Geometry polygon) { List rs = new ArrayList<>(); // return rs; } public List 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 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 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 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 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 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); } }