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 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 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); } catch (Exception ex) { setError(entity, ex.getMessage()); log.error(ex.getMessage(), ex); } finally { // gdal.GDALDestroyDriverManager() if (null != ds) { ds.delete(); } } } /** * 处理结果 */ 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) { double min = Collections.min(rs); double max = Collections.max(rs); double avg = rs.stream().mapToDouble(Double::valueOf).average().getAsDouble(); entity.setMin(min); entity.setMax(max); entity.setAvg(avg); } /** * 分析线 */ public List analysisPolyline(Geometry polyline) { 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()); processPolyline(entity, pub, polyline); rs.add(entity); } return rs; } /** * 处理线 */ private void processPolyline(AnalysisResultEntity entity, PublishEntity pub, Geometry polyline) { 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; } processPolyline(entity, filePath, polyline); } /** * 处理线 */ public void processPolyline(AnalysisResultEntity entity, String filePath, Geometry polyline) { Dataset ds = null; try { ds = gdal.Open(filePath); if (null == ds) { throw new Exception("打开栅格数据失败"); } 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[] bounds = new double[6]; polyline.GetEnvelope(bounds); double startX = bounds[0], endX = bounds[1], startY = bounds[2], endY = bounds[3]; int bandCount = ds.getRasterCount(); int xStartPixel = (int) Math.floor((startX - minX) / pixelWidth); int yStartPixel = (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 list = new ArrayList<>(); for (int y = Math.min(yStartPixel, yEndPixel); y <= Math.max(yStartPixel, 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 (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); } catch (Exception ex) { setError(entity, ex.getMessage()); log.error(ex.getMessage(), ex); } finally { // gdal.GDALDestroyDriverManager() if (null != ds) { ds.delete(); } } } /** * 分析面 */ public List analysisPolygon(Geometry polygon) { List rs = new ArrayList<>(); // return rs; } 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 0; } 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]; } } }