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 analysis(Geometry geo, Integer 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()); List 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); 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 point, int size) { 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 bandCount = ds.getRasterCount(); int xPixel = (int) Math.floor((x - minX) / pixelWidth); int yPixel = (int) Math.floor((maxY - y) / Math.abs(pixelHeight)); List 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 polyline) { 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]; double[] env = new double[4]; polyline.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.floor((env[1] - minX) / pixelWidth); int yMaxPixel = (int) Math.floor((maxY - env[2]) / Math.abs(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 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 (polyline.Contains(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 polygon) { 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]; double[] env = new double[4]; polygon.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.floor((env[1] - minX) / pixelWidth); int yMaxPixel = (int) Math.floor((maxY - env[2]) / Math.abs(pixelHeight)); int bandCount = ds.getRasterCount(); int geoWidth = Math.abs(xMaxPixel - xMinPixel); int geoHeight = Math.abs(yMaxPixel - yMinPixel); List 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 list) { if (null == list || list.isEmpty()) { return; } List 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 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 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 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]; } } }