PM20230220027_Service.git

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<AnalysisResultEntity> analysis(Geometry geo, Integer size) {
        List<AnalysisResultEntity> rs = new ArrayList<>();
 
        List<PublishEntity> pubs = publishService.selectRaster();
        if (null == pubs || pubs.isEmpty()) {
            return rs;
        }
 
        for (PublishEntity pub : pubs) {
            AnalysisResultEntity entity = new AnalysisResultEntity();
            entity.setLayerName(pub.getName());
 
            List<MetaEntity> 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<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);
    }
 
    /**
     * 分析线
     */
    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[] bounds = new double[4];
        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<Double> 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);
    }
 
    /**
     * 分析面
     */
    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[2]) / Math.abs(pixelHeight));
        int xMaxPixel = (int) Math.floor((env[1] - minX) / pixelWidth);
        int yMaxPixel = (int) Math.floor((maxY - env[3]) / 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) {
        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<double[]> 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<String> 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];
        }
    }
}