package com.moon.server.service.data; import com.alibaba.fastjson.JSONObject; 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.GeoHelper; import com.moon.server.helper.PathHelper; import com.moon.server.helper.WebHelper; import com.moon.server.service.all.WebSocketService; 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.awt.geom.Point2D; import java.io.File; import java.io.IOException; 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); /** * 测试 * * POINT (26.72 -48.58) * LINESTRING(18.44 -46.69,20.14 -45.93) * POLYGON ((5.11 -41.64,7.68 -43.59,4.73 -43.3,5.11 -41.64)) * POLYGON ((42.6 -40.3,49.1 -40.4,49.1 -65.7,46.3 -65.6,42.6 -40.3)) */ public List test(Geometry geo, Integer size) { List rs = new ArrayList<>(); AnalysisResultEntity entity = new AnalysisResultEntity(); entity.setLayerName("Test"); openRaster(entity, "D:\\Moon\\data\\DOM\\Lunar_LRO_LOLA_ClrShade_Global_128ppd_v04_2.tif", geo, size); rs.add(entity); return rs; } /** * 使用WKT查询分析,结果以消息推送 */ public void analysisForPost(Geometry geo, Integer size, String token) { List pubs = publishService.selectRaster(); if (null == pubs || pubs.isEmpty()) { return; } // for (PublishEntity pub : pubs) { pubs.parallelStream().forEach(pub -> { try { AnalysisResultEntity entity = new AnalysisResultEntity(token); entity.setLayerName(pub.getName()); List metas = publishService.selectMetasByPubid(pub.getId()); if (null == metas || metas.isEmpty()) { setError(entity, "找不到发布数据"); postInfo(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, "源数据不存在"); postInfo(entity); return; } openRaster(entity, filePath, geo, size); postInfo(entity); } catch (Exception ex) { log.error(ex.getMessage(), ex); } }); } /** * 推送消息 */ private void postInfo(AnalysisResultEntity entity) throws IOException { JSONObject map = new JSONObject(); map.put("analysisForPost", entity); String json = JSONObject.toJSONString(map); // System.out.println(json) WebSocketService.broadCastInfo(json); } /** * 分析方法 */ 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 { // filePath = "D:\\Moon\\data\\DOM\\test.tif" 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, size); 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 geo, int size) { double[] transform = ds.GetGeoTransform(); double pixelWidth = transform[1], pixelHeight = Math.abs(transform[5]); double buffer = Math.max(pixelWidth, pixelHeight) * size; geo = geo.Buffer(buffer); analysisPolygon(entity, ds, geo); } /** * 分析线 */ public void analysisPolyline(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) { double[] transform = ds.GetGeoTransform(); int xSize = ds.getRasterXSize(), ySize = ds.getRasterYSize(); double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = transform[5]; double len = 0; int bandCount = ds.getRasterCount(); double lineDis = getPolylineDistance(geo); double distance = lineDis / Math.max(geo.GetPointCount(), size); List points = getPointsByDistance(geo, distance); for (int i = 0, c = points.size(); i < c; i++) { double[] xy = points.get(i); if (i > 0) { double[] p = points.get(i - 1); len += GeoHelper.getDistance(p[0], p[1], xy[0], xy[1]); } int xPixel = (int) Math.floor((xy[0] - minX) / pixelWidth); int yPixel = (int) Math.floor((maxY - xy[1]) / Math.abs(pixelHeight)); if (xPixel < 1 || xPixel > xSize || yPixel < 1 || yPixel > ySize) { continue; } List vals = new ArrayList<>(); for (int j = 1; j <= bandCount; j++) { double[] pixelValues = new double[1]; ds.GetRasterBand(j).ReadRaster(xPixel, yPixel, 1, 1, pixelValues); if (!Double.isNaN(pixelValues[0])) { vals.add(WebHelper.round(pixelValues[0], 3)); } } if (vals.size() > 0) { entity.addPoint(xy[0], xy[1], len, vals); } } } /** * 分析面 */ public void analysisPolygon(AnalysisResultEntity entity, Dataset ds, Geometry geo) { double[] transform = ds.GetGeoTransform(); int xSize = ds.getRasterXSize(), ySize = ds.getRasterYSize(); double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = Math.abs(transform[5]), rotationX = transform[2], rotationY = transform[4]; double[] env = new double[4]; geo.GetEnvelope(env); entity.addPoint(geo.Centroid().GetX(), geo.Centroid().GetY(), 0, null); int xMinPixel = Math.max((int) Math.floor((env[0] - minX) / pixelWidth), 1); int yMinPixel = Math.max((int) Math.floor((maxY - env[3]) / pixelHeight), 1); int xMaxPixel = Math.min((int) Math.floor((env[1] - minX) / pixelWidth), xSize); int yMaxPixel = Math.min((int) Math.floor((maxY - env[2]) / pixelHeight), ySize); if (xMaxPixel < 1 || yMaxPixel < 1 || xMaxPixel - xMinPixel < 0 || yMaxPixel - yMinPixel < 0) { setError(entity, "查询范围无效"); return; } int bandCount = ds.getRasterCount(); int width = xMaxPixel - xMinPixel; int height = yMaxPixel - yMinPixel; if (width * height > StaticData.I64 * StaticData.I64) { readRasterForBlocks(entity, ds, bandCount, xMinPixel, yMinPixel, width, height); return; } for (int i = 1; i <= bandCount; i++) { double[] pixelValues = new double[width * height]; ds.GetRasterBand(i).ReadRaster(xMinPixel, yMinPixel, width, height, pixelValues); setBandVals(entity, pixelValues); } } /** * 按照块读取栅格数据 */ private void readRasterForBlocks(AnalysisResultEntity entity, Dataset ds, int bandCount, int xMinPixel, int yMinPixel, int width, int height) { List xList = getSamples(xMinPixel, width - 1); List yList = getSamples(yMinPixel, height - 1); double[] pixelValues = new double[1]; for (int i = 1; i <= bandCount; i++) { List list = new ArrayList<>(); for (Integer x : xList) { for (Integer y : yList) { ds.GetRasterBand(i).ReadRaster(x, y, 1, 1, pixelValues); if (!Double.isNaN(pixelValues[0])) { list.add(pixelValues[0]); } } } setBandVals(entity, list); } } /** * 获取抽样列表 */ private List getSamples(int start, int strip) { List list = new ArrayList<>(); double avg = 1.0 * strip / (StaticData.I64 - 1); for (int i = 0; i < StaticData.I64; i++) { list.add(start + (int) Math.ceil(avg * i)); } return list; } /** * 设置错误 */ private void setError(AnalysisResultEntity entity, String info) { entity.setCode(StaticData.I500); entity.setInfo(info); } /** * 设置Band值 */ private void setBandVals(AnalysisResultEntity entity, double[] pixelValues) { if (null == pixelValues || pixelValues.length == 0) { return; } List list = new ArrayList<>(); for (double val : pixelValues) { if (!Double.isNaN(val)) { list.add(val); } } setBandVals(entity, list); } /** * 设置Band值 */ private void setBandVals(AnalysisResultEntity entity, List list) { if (null == list || list.isEmpty()) { return; } double min = Collections.min(list); double max = Collections.max(list); double avg = list.stream().mapToDouble(Double::valueOf).average().getAsDouble(); entity.addBandVals(min, avg, max); } /** * 获取折线距离 */ private double getPolylineDistance(Geometry geo) { double dis = 0; for (int i = 0, c = geo.GetPointCount() - 1; i < c; i++) { double[] xy1 = geo.GetPoint(i); double[] xy2 = geo.GetPoint(i + 1); dis += GeoHelper.getDistance(xy1[0], xy1[1], xy2[0], xy2[1]); } return dis; } /** * 根据距离获取节点 */ private List getPointsByDistance(Geometry geo, double distance) { List list = new ArrayList<>(); int c = geo.GetPointCount() - 1; for (int i = 0; i < c; i++) { double[] xy1 = geo.GetPoint(i); double[] xy2 = geo.GetPoint(i + 1); list.add(xy1); double lineDis = GeoHelper.getDistance(xy1[0], xy1[1], xy2[0], xy2[1]); if (lineDis < distance) { continue; } double d = distance; double angle = GeoHelper.getAngle(xy1[0], xy1[1], xy2[0], xy2[1]); while (d + distance * StaticData.D05 < lineDis) { Point2D point = GeoHelper.getPointByDistanceAndAngle(xy1[0], xy1[1], angle, d); list.add(new double[]{point.getX(), point.getY()}); d += distance; } } list.add(geo.GetPoint(c)); return list; } /** * 分析点 2 */ public void analysisPoint2(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) { double x = geo.GetX(), y = geo.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)); 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); } } /** * 分析线 2 */ public void analysisPolyline2(AnalysisResultEntity entity, Dataset ds, Geometry geo, int size) { double[] transform = ds.GetGeoTransform(); // double rotationX = transform[2]; double rotationY = transform[4] double minX = transform[0], pixelWidth = transform[1], maxY = transform[3], pixelHeight = Math.abs(transform[5]); double buffer = Math.max(pixelWidth, pixelHeight) * size; double[] env = new double[4]; geo = geo.Buffer(buffer); geo.GetEnvelope(env); int xMinPixel = (int) Math.floor((env[0] - minX) / pixelWidth); int yMinPixel = (int) Math.floor((maxY - env[3]) / pixelHeight); int xMaxPixel = (int) Math.floor((env[1] - minX) / pixelWidth); int yMaxPixel = (int) Math.floor((maxY - env[2]) / 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 (geo.Intersects(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) } private 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); } } private 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)); } }