屏蔽NoData

13693261870

2023-12-04 45a773850057dd90c4292d30715b4f9dfdc86740

 src/main/java/com/moon/server/service/data/RasterAnalysisService.java

@@ -1,10 +1,15 @@
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.StringHelper;
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;
@@ -13,10 +18,14 @@
import org.gdal.gdal.gdal;
import org.gdal.ogr.Geometry;
import org.gdal.ogr.ogr;
import org.springframework.beans.factory.annotation.Value;
import org.springframework.stereotype.Service;

import javax.annotation.Resource;
import java.awt.geom.Point2D;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.*;

/**
@@ -26,6 +35,9 @@
 */
@Service
public class RasterAnalysisService {
    @Value("${sys.gdal_path}")
    String gdalPath;

    @Resource
    PathHelper pathHelper;

@@ -33,6 +45,102 @@
    PublishService publishService;

    private final static Log log = LogFactory.getLog(RasterAnalysisService.class);

    /**
     * 测试
     * <p>
     * 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<AnalysisResultEntity> test(Geometry geo, Integer size) {
        List<AnalysisResultEntity> 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<PublishEntity> 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());
                entity.setUnit(pub.getUnit());

                List<MetaEntity> metas = publishService.selectMetasByPubid(pub.getId());
                if (null == metas || metas.isEmpty()) {
                    setError(entity, "找不到发布数据");
                    postInfo(entity);
                    return;
                }

                String filePath = getFilePath(metas);
                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 String getFilePath(List<MetaEntity> metas) throws IOException {
        String basePath = pathHelper.getConfig().getUploadPath() + File.separator;
        if (metas.size() == 1) {
            return basePath + metas.get(0).getPath();
        }

        String tempPath = pathHelper.getTempPath() + File.separator, guid = StringHelper.getGuid();
        String fileList = tempPath + guid + ".txt";
        FileWriter fw = new FileWriter(fileList);
        for (int i = 0, c = metas.size(); i < c; i++) {
            String str = (i > 0 ? "\r\n" : "") + basePath + metas.get(i).getPath();
            fw.write(str.toCharArray());
        }
        fw.close();

        String vrt = tempPath + guid + ".vrt";
        String cmd = String.format("\"%s\\gdalbuildvrt.exe\" -input_file_list \"%s\" \"%s\"", gdalPath, fileList, vrt);
        WebHelper.exec(cmd);

        return vrt;
    }

    /**
     * 推送消息
     */
    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);
    }

    /**
     * 分析方法
@@ -46,26 +154,30 @@
        }

        for (PublishEntity pub : pubs) {
            AnalysisResultEntity entity = new AnalysisResultEntity();
            entity.setLayerName(pub.getName());
            try {
                AnalysisResultEntity entity = new AnalysisResultEntity();
                entity.setLayerName(pub.getName());
                entity.setUnit(pub.getUnit());

            List<MetaEntity> metas = publishService.selectMetasByPubid(pub.getId());
            if (null == metas || metas.isEmpty()) {
                setError(entity, "找不到发布数据");
                continue;
                List<MetaEntity> metas = publishService.selectMetasByPubid(pub.getId());
                if (null == metas || metas.isEmpty()) {
                    setError(entity, "找不到发布数据");
                    continue;
                }

                String filePath = getFilePath(metas);
                File file = new File(filePath);
                if (!file.exists() || file.isDirectory()) {
                    setError(entity, "源数据不存在");
                    continue;
                }

                openRaster(entity, filePath, geo, size);
                rs.add(entity);
            } catch (Exception ex) {
                log.error(ex.getMessage(), ex);
            }

            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;
    }

@@ -75,6 +187,7 @@
    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("打开栅格数据失败");
@@ -85,7 +198,7 @@
                    analysisPoint(entity, ds, geo, size);
                    break;
                case ogr.wkbLineString:
                    analysisPolyline(entity, ds, geo);
                    analysisPolyline(entity, ds, geo, size);
                    break;
                default:
                    analysisPolygon(entity, ds, geo);
@@ -104,8 +217,217 @@
    /**
     * 分析点
     */
    public void analysisPoint(AnalysisResultEntity entity, Dataset ds, Geometry point, int size) {
        double x = point.GetX(), y = point.GetY();
    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<double[]> 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<Double> 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<Integer> xList = getSamples(xMinPixel, width - 1);
        List<Integer> yList = getSamples(yMinPixel, height - 1);

        double[] pixelValues = new double[1];
        for (int i = 1; i <= bandCount; i++) {
            List<Double> 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<Integer> getSamples(int start, int strip) {
        List<Integer> 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<Double> list = new ArrayList<>();
        for (double val : pixelValues) {
            if (!Double.isNaN(val) && val > Integer.MIN_VALUE) {
                list.add(val);
            }
        }
        setBandVals(entity, list);
    }

    /**
     * 设置Band值
     */
    private void setBandVals(AnalysisResultEntity entity, List<Double> 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<double[]> getPointsByDistance(Geometry geo, double distance) {
        List<double[]> 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];
@@ -114,32 +436,30 @@
        int xPixel = (int) Math.floor((x - minX) / pixelWidth);
        int yPixel = (int) Math.floor((maxY - y) / Math.abs(pixelHeight));

        List<double[]> 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);
    }

    /**
     * 分析线
     * 分析线 2
     */
    public void analysisPolyline(AnalysisResultEntity entity, Dataset ds, Geometry polyline) {
    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 = transform[5];
        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];
        polyline.GetEnvelope(env);
        geo = geo.Buffer(buffer);
        geo.GetEnvelope(env);

        int xMinPixel = (int) Math.floor((env[0] - minX) / pixelWidth);
        int yMinPixel = (int) Math.floor((maxY - env[3]) / Math.abs(pixelHeight));
        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]) / Math.abs(pixelHeight));
        int yMaxPixel = (int) Math.floor((maxY - env[2]) / pixelHeight);

        int bandCount = ds.getRasterCount();
        int xMin = Math.min(xMinPixel, xMaxPixel);
@@ -152,7 +472,7 @@
            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)) {
                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);
@@ -161,84 +481,10 @@
                }
            }
        }

        processResult(list, entity);
        // 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<double[]> 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<double[]> list) {
        if (null == list || list.isEmpty()) {
            return;
        }

        List<Double> 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<Double> 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) {
    private void processClippedDataByLine(String imagePath, String geometryString) {
        // 注册GDAL驱动
        gdal.AllRegister();
        // 打开栅格图像数据集
@@ -266,7 +512,7 @@
        }
    }

    public static Dataset clipRaster(Dataset dataset, Geometry geometry) {
    private Dataset clipRaster(Dataset dataset, Geometry geometry) {
        Vector<String> warpOptions = new Vector<>();
        warpOptions.add("-crop_to_cutline");
        warpOptions.add("-cutline");
@@ -274,34 +520,5 @@
        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];
        }
    }
}