JAVA后端——依据离散点/格点生成GEOJSON以渲染色斑图

JAVA后端——依据离散点/格点生成GEOJSON以渲染色斑图

纯干货，一文既懂，相关解释都放进了代码的注释里，不多赘述

所需依赖

1.org.geotools相关依赖
2.jts相关依赖
3.netcdfAll相关依赖
4.其他常规JAVA依赖

代码

对geojson的源文件进行解析，以某一个要素为例，需要先获取对应的文件（一般为栅格文件），解析并取出对应的数据，形成栅格数据或进行特殊处理降噪后的离散点数据，再进行插值。
以下代码是以第一个时间维度为例做的解析，实际要根据需求修改代码

public static void main(String[] args) throws IOException {
        //生成geojosn的对应文件路径
        String filePath = "filePath";
        //读取文件
        NetcdfFile ncFile = NetcdfFile.open(filePath);

        // 打印变量信息
        NcFileUtils.printVariableInfo(ncFile, "x");

        // 获取维度信息
        int[] shape = NcFileUtils.getArrayShape(ncFile, "x");
        log.info("原始维度: {}", java.util.Arrays.toString(shape));

        // 读取经纬度（一维数组）
        float[] latArray = NcFileUtils.formatArray(NcFileUtils.readVariableTo1D(ncFile, "lat"));
        float[] lonArray = NcFileUtils.formatArray(NcFileUtils.readVariableTo1D(ncFile, "lon"));
        log.info("经度数量: {}, 纬度数量: {}", lonArray.length, latArray.length);

        // 根据维度选择合适的读取方法
        // 读取三维数据（保持完整结构）
        float[][][] data3D = NcFileUtils.formatArray3D(NcFileUtils.readVariableTo3D(ncFile, "x"));
        log.info("三维数据形状: {} x {} x {}", data3D.length, data3D[0].length, data3D[0][0].length);

        // 经纬度范围（根据边界数组计算的最值）
        double xMin = HUNAN_LON_MIN;   // 最小经度（西边界）
        double xMax = HUNAN_LON_MAX;   // 最大经度（东边界）
        double yMin = HUNAN_LAT_MIN;    // 最小纬度（南边界）
        double yMax = HUNAN_LAT_MAX;    // 最大纬度（北边界）

        //1: 3.14日 9: 3.22日
        float[][] data1 = data3D[0];
        List<Double> latList = new ArrayList<>();
        List<Double> lonList = new ArrayList<>();
        List<Double> valueList = new ArrayList<>();

        for (int i = 0; i < data1.length; i++) {

            double lat = latArray[i];
            for (int j = 0; j < data1[i].length; j++) {
                double lon = lonArray[j];
                if (lat < yMin || lat > yMax || lon < xMin || lon > xMax) {
                    continue;
                }
                double value = data1[i][j];
                latList.add(lat);
                lonList.add(lon);
                valueList.add(value);
//                log.info("经度: {}, 纬度: {}, 值: {}", lon, lat, value);
            }

        }
        try {
            ncFile.close();
        } catch (IOException e) {
            throw new RuntimeException(e);
        }

        double[] xs = new double[lonList.size()];
        double[] ys = new double[latList.size()];
        double[] values = new double[valueList.size()];
        for (int i = 0; i < lonList.size(); i++) {
            xs[i] = BigDecimalUtil.roundToDecimalPlaces(lonList.get(i), 3);
            ys[i] = BigDecimalUtil.roundToDecimalPlaces(latList.get(i), 3);
            values[i] = BigDecimalUtil.roundToDecimalPlaces(valueList.get(i), 3);
        }
        log.info("xs size:{} ys size:{}", xs.length, ys.length);

是否需要对最后的geojson进行裁剪

        ShpBoundaryUtil.BoundaryResult boundaryResult = ShpBoundaryUtil.getBoundaryByFeatures(HUNAN_SHP_PATH);
        double[] clipYs = boundaryResult.getLatitudes();
        double[] clipXs = boundaryResult.getLongitudes();        

预期的栅格数据量大小，根据实际情况调整

        int xNum = 300;   // 经度方向格点数（更密集的网格）
        int yNum = 300;   // 纬度方向格点数

		// 无效值过滤
        double undefValue = -9999.0;

注意：将对数据进行插值，以下是几种算法的特点
/**
* IDW 临近点算法、固定最近N离散点数的插值，点约近权重越大
*/
IDW_NEIGHBOR,

/**
 * IDW 范围查找（kdTree 实现）检索以该点为中心、指定半径范围内的所有离散点，再用 IDW 计算，幂次衰减：权重随距离幂次（p）单调下降，无严格边界（即是受所有数据影响）
 * 属于 IDW_RADIUS中大数据量的优化方案
 */
IDW_RADIUS_KD_TREE,

/**
 * Cressman kdTree 指定半径范围内的所有离散点，非线性衰减：距离越近，权重增速越快，有明确边界（受范围内数据影响） 
 * Cressman中大数据量的优化方案
 */
CRESSMAN_KD_TREE,

插值并生成结果

        // 插值算法（从离散点插值到网格）
        InterpolateAlgorithm algorithm = InterpolateAlgorithm.IDW_NEIGHBOR;
//        InterpolateAlgorithm algorithm = null;

        // 定义颜色区间（适合气温范围）
        InterpolateColor.Interval interval = InterpolateColor.PRE_24H;

        // 输出类型：POLYGON（等值面）、POLYLINE（等值线）、IMAGE（图像）
        ContourOutType type = ContourOutType.POLYGON;
        // 调用 toJSON 方法
        JSONObject result = ContourUtil.toJSON(
                xs,
                ys,
                values,
                clipXs,
                clipYs,
                xMin,
                xMax,
                yMin,
                yMax,
                xNum,
                yNum,
                undefValue,
                algorithm,
                interval,
                type
        );

        // 输出结果
        System.out.println("生成的离散点插值 GeoJSON: ");
        System.out.println(result.toJSONString()
                .substring(0, Math.min(result.toJSONString().length(), 100)) +" ...");

        // 保存到文件
        ContourUtil.saveToFile("xxxxxxxx\\output.json", result);

若是离散数据，按照此方法执行，获取对应值范围的色卡，即geojosn的色卡，并处理离散点数据，数据分级以便于生成geojson

   /**
     * 离散点数据轮廓处理
     *
     * @param xs         经度
     * @param ys         纬度
     * @param values     观测值
     * @param clipXs     裁剪边界经度
     * @param clipYs     裁剪边界纬度
     * @param xMin       最小经度
     * @param xMax       最大经度
     * @param yMin       最小纬度
     * @param yMax       最大纬度
     * @param xNum       径向格点数
     * @param yNum       纬向格点数
     * @param undefValue 无效值
     * @param algorithm  插值算法
     * @param interval   轮廓色标值区间
     * @param type       输出类型
     * @return 结果
     */
    public JSONObject toJSON(
        double[] xs,
        double[] ys,
        double[] values,
        double[] clipXs,
        double[] clipYs,
        double xMin,
        double xMax,
        double yMin,
        double yMax,
        int xNum,
        int yNum,
        double undefValue,
        InterpolateAlgorithm algorithm,
        InterpolateColor.Interval interval,
        ContourOutType type
    ) {
        // 轮廓区间值
        double[] contourValues = interval.getValues();
        // 处理离散数据
        double[][] discreteData = processDiscreteValues(xs, ys, values, contourValues);

        return toJSON(
            discreteData,
            clipXs,
            clipYs,
            xMin,
            xMax,
            yMin,
            yMax,
            xNum,
            yNum,
            undefValue,
            algorithm,
            interval,
            type
        );
    }

    public double[][] processDiscreteValues(double[] xs, double[] ys, double[] values, double[] intervalValues) {
        double[][] discreteData = new double[xs.length][3];
        for (int i = 0; i < xs.length; i++) {
            discreteData[i][0] = xs[i];
            discreteData[i][1] = ys[i];
            if (values[i] <= intervalValues[0]) {
                discreteData[i][2] = intervalValues[0] - 1;
                continue;
            }
            if (values[i] >= intervalValues[intervalValues.length - 1]) {
                discreteData[i][2] = 1.5 * intervalValues[intervalValues.length - 1];
                continue;
            }
            for (int j = 0; j < intervalValues.length - 1; j++) {
                if (values[i] >= intervalValues[j] && values[i] < intervalValues[j + 1]) {
                    discreteData[i][2] = 0.5 * (intervalValues[j] + intervalValues[j + 1]);
                    break;
                }
            }
        }

        return discreteData;
    }

以分级后的数据生成插值获得网格数据、计算平滑轮廓线，是最后的geojson更加自然，最后输出等值面json

    /**
     * 离散点数据轮廓处理
     *
     * @param clipXs     裁剪边界经度
     * @param clipYs     裁剪边界纬度
     * @param xMin       最小经度
     * @param xMax       最大经度
     * @param yMin       最小纬度
     * @param yMax       最大纬度
     * @param xNum       径向格点数
     * @param yNum       纬向格点数
     * @param undefValue 无效值
     * @param algorithm  插值算法
     * @param interval   轮廓色标值区间
     * @param type       输出类型
     * @return 结果
     */
    public JSONObject toJSON(
        double[][] gridData,
        double[] clipXs,
        double[] clipYs,
        double xMin,
        double xMax,
        double yMin,
        double yMax,
        int xNum,
        int yNum,
        double undefValue,
        InterpolateAlgorithm algorithm,
        InterpolateColor.Interval interval,
        ContourOutType type
    ) {
        // 轮廓区间值
        double[] contourValues = interval.getValues();
        // 创建栅格坐标
        double[] gridX = new double[xNum];
        double[] gridY = new double[yNum];
        Interpolate.createGridXY_Num(xMin, yMin, xMax, yMax, gridX, gridY);

        long start = System.currentTimeMillis();
        if (algorithm != null) {
            // 插值获得网格数据
            gridData = interpolate(gridData, gridX, gridY, undefValue, algorithm);
        }
        log.info("插值花费时间{}", (System.currentTimeMillis() - start));
        // 追踪计算网格边界
        int[][] temp = new int[gridData.length][gridData[0].length];
        List<Border> borders = Contour.tracingBorders(gridData, gridX, gridY, temp, undefValue);

        // 追踪计算轮廓线
        List<PolyLine> contourLines = Contour.tracingContourLines(
            gridData,
            gridX,
            gridY,
            contourValues.length,
            contourValues,
            undefValue,
            borders,
            temp
        );

        // 平滑处理轮廓线
        List<PolyLine> lines = Contour.smoothLines(contourLines);
        if (ContourOutType.POLYLINE.equals(type)) {
            return toPolyLineJson(lines, clipXs, clipYs, interval);
        }

        // 追踪计算轮廓面
        List<Polygon> polygons = Contour.tracingPolygons(gridData, lines, borders, contourValues);
        // 输出等值面 GeoJSON
        if (ContourOutType.POLYGON.equals(type)) {
            return toPolygonJson(polygons, clipXs, clipYs, interval);
        }
        return new JSONObject();
    }

插值方法

    /**
     * 根据算法将离散点数据插值成网格数据
     *
     * @param discreteData 离散点数据
     * @param gridX        网格 X 坐标集合（eg: Lon List）
     * @param gridY        网格 Y 坐标集合（eg: Lat List）
     * @param undefValue   默认值
     * @param algorithm    算法 {@link InterpolateAlgorithm}
     * @return 插值后的网格数据
     */
    public double[][] interpolate(
        double[][] discreteData,
        double[] gridX,
        double[] gridY,
        double undefValue,
        InterpolateAlgorithm algorithm
    ) {
        // 定义网格数据
        double[][] gridData = new double[gridX.length][gridY.length];

        // 根据指定算法进行插值处理
        if (InterpolateAlgorithm.IDW_NEIGHBOR.equals(algorithm)) {
            gridData = Interpolate.interpolation_IDW_Neighbor(discreteData, gridX, gridY, 4, undefValue);
        } else if (InterpolateAlgorithm.IDW_RADIUS.equals(algorithm)) {
            gridData = Interpolate.interpolation_IDW_Radius(discreteData, gridX, gridY, 4, 2, undefValue);
        } else if (InterpolateAlgorithm.IDW_RADIUS_KD_TREE.equals(algorithm)) {
            gridData = Interpolate.idw_Radius_kdTree(discreteData, gridX, gridY, 4, 2, undefValue);
        } else if (InterpolateAlgorithm.CRESSMAN.equals(algorithm)) {
            List<Double> list = Arrays.asList(10D, 8D, 6D, 4D, 2D);
            gridData = Interpolate.cressman(discreteData, gridX, gridY, undefValue, list);
        } else if (InterpolateAlgorithm.CRESSMAN_KD_TREE.equals(algorithm)) {
            List<Double> list = Arrays.asList(10D, 8D, 6D, 4D, 2D);
            gridData = Interpolate.cressman_kdTree(discreteData, gridX, gridY, undefValue, list);
        }

        return gridData;
    }

根据interval拼geojosn的方法

    /**
     * 将 Polygons 根据数值区间转换为 GeoJSON
     *
     * @param polygons 
     * @param interval 
     * @return GeoJSON
     */
    private static JSONObject toPolygonJson(
        List<Polygon> polygons,
        double[] clipXs,
        double[] clipYs,
        InterpolateColor.Interval interval
    ) {
        List<PointD> clips = getClipPointList(clipXs, clipYs);
        if (!clips.isEmpty()) {
            polygons = Contour.clipPolygons(polygons, clips);
        }
        String features = polygons
            .stream()
            .map(polygon -> {
                List<String> lines = new ArrayList<>();
                lines.add("\"type\":\"Feature\"");
                List<String> coordinates = new ArrayList<>();
                PointD pointFirst = polygon.OutLine.PointList.get(0);
                Double firstX = NumberUtil.round(pointFirst.X, 6).doubleValue();
                Double firstY = NumberUtil.round(pointFirst.Y, 6).doubleValue();
                PointD pointLast = polygon.OutLine.PointList.get(polygon.OutLine.PointList.size() - 1);
                Double lastX = NumberUtil.round(pointLast.X, 6).doubleValue();
                Double lastY = NumberUtil.round(pointLast.Y, 6).doubleValue();
                double absX = Math.abs(NumberUtil.sub(firstX, lastX));
                if (absX > 0) {
                    pointLast.X = firstX;
                }
                double absY = Math.abs(NumberUtil.sub(firstY, lastY));
                if (absY > 0) {
                    pointLast.Y = firstY;
                }
                coordinates.add(
                    "[" +
                    polygon.OutLine.PointList
                        .stream()
                        .map(point ->
                            "[" +
                            NumberUtil.round(point.X, 6).doubleValue() +
                            "," +
                            NumberUtil.round(point.Y, 6).doubleValue() +
                            "]"
                        )
                        .collect(Collectors.joining(",")) +
                    "]"
                );
                if (polygon.HasHoles()) {
                    coordinates.add(
                        polygon.HoleLines
                            .stream()
                            .map(line -> {
                                PointD pointFirst1 = line.PointList.get(0);
                                Double firstX1 = NumberUtil.round(pointFirst1.X, 6).doubleValue();
                                Double firstY1 = NumberUtil.round(pointFirst1.Y, 6).doubleValue();
                                PointD pointLast1 = line.PointList.get(line.PointList.size() - 1);
                                Double lastX1 = NumberUtil.round(pointLast1.X, 6).doubleValue();
                                Double lastY1 = NumberUtil.round(pointLast1.Y, 6).doubleValue();
                                double absX1 = Math.abs(NumberUtil.sub(firstX1, lastX1));
                                if (absX1 > 0 && absX1 < 0.0001) {
                                    pointLast1.X = firstX1;
                                }
                                double absY1 = Math.abs(NumberUtil.sub(firstY1, lastY1));
                                if (absY1 > 0 && absY1 < 0.0001) {
                                    pointLast1.Y = firstY1;
                                }
                                return (
                                    "[" +
                                    line.PointList
                                        .stream()
                                        .map(point ->
                                            (
                                                "[" +
                                                NumberUtil.round(point.X, 6) +
                                                "," +
                                                NumberUtil.round(point.Y, 6) +
                                                "]"
                                            )
                                        )
                                        .collect(Collectors.joining(",")) +
                                    "]"
                                );
                            })
                            .collect(Collectors.joining(","))
                    );
                }
                lines.add(
                    "\"geometry\":{\"type\":\"Polygon\",\"coordinates\":[" + String.join(",", coordinates) + "]}"
                );

                String color = getRangeColor(interval, polygon.LowValue, polygon.IsHighCenter);
                lines.add("\"properties\":{\"color\":\"" + color + "\"}");

                return "{" + String.join(",", lines) + "}";
            })
            .collect(Collectors.joining(","));
        String result = "{\"type\":\"FeatureCollection\",\"features\":[" + features + "]}";
        return JSONObject.parseObject(result);
    }

网格数据不过多解释，代码如下

下集更精彩~