java实现反距离权重插值算法生成geojson矢量数据

一、java代码

package cn.cepec.talroad.air.common;

import cn.cepec.talroad.air.common.geoutil.FeaureUtil;
import cn.cepec.talroad.air.common.geoutil.GeoJSONUtil;
import com.vividsolutions.jts.geom.Geometry;
import net.sf.json.JSONObject;
import org.geotools.data.FeatureSource;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.geotools.geojson.feature.FeatureJSON;
import org.opengis.feature.Feature;
import org.opengis.feature.simple.SimpleFeature;
import wcontour.Contour;
import wcontour.global.Border;
import wcontour.global.PointD;
import wcontour.global.PolyLine;
import wcontour.global.Polygon;
import wcontour.Interpolate;

import java.io.File;
import java.io.StringWriter;
import java.nio.charset.Charset;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class InterpolationUtils {
    
    // 解析Shapefile文件
    public static String calEquiSurface(double[][] trainData, double[] dataInterval, int[] size, String boundryFile,
                                 boolean isclip) {
        String geojsonpogylon = "";
        try {
            double _undefData = -9999.0;
            // 多边形集合
            SimpleFeatureCollection polygonCollection = null;
            // 多边形线集合
            List cPolylineList = new ArrayList();
            // 多边形List
            List cPolygonList = new ArrayList();

            int width = size[0], height = size[1];
            double[] _X = new double[width];
            double[] _Y = new double[height];
            File file = new File(boundryFile);

            // 源shape文件
            ShapefileDataStore shpDataStore = null;
            shpDataStore = new ShapefileDataStore(file.toURL());
            // 设置编码
            Charset charset = Charset.forName("GBK");
            shpDataStore.setCharset(charset);
            // 获取文件后缀名
            String typeName = shpDataStore.getTypeNames()[0];
            // SimpleFeatureSource相当于AE中的featureClass，设置属性
            SimpleFeatureSource featureSource = null;
            featureSource = shpDataStore.getFeatureSource();
            SimpleFeatureCollection fc = featureSource.getFeatures();
            // 最大最小经纬度
            double minX = fc.getBounds().getMinX();
            double minY = fc.getBounds().getMinY();
            double maxX = fc.getBounds().getMaxX();
            double maxY = fc.getBounds().getMaxY();

            Interpolate.createGridXY_Num(minX, minY, maxX, maxY, _X, _Y);
            double[][] _gridData = new double[width][height];

            // 数据间隔长度
            int nc = dataInterval.length;
            // IDW插值 （训练数据（离散数据阵列）、宽（栅格X阵）、高（栅格Y阵）、默认数（最近邻居数））
            _gridData = Interpolate.interpolation_IDW_Neighbor(trainData, _X, _Y, 12, _undefData);

            int[][] S1 = new int[_gridData.length][_gridData[0].length];
            // 边界
            List _borders = Contour.tracingBorders(_gridData, _X, _Y, S1, _undefData);
            // 生成等值线（IDW插值、宽、高、数据间隔长度、数据间隔、默认值、边界、IDW长度）
            cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc, dataInterval, _undefData, _borders, S1);
            // 平滑
            cPolylineList = Contour.smoothLines(cPolylineList);

            cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList, _borders, dataInterval);
            // 多边形Json
            geojsonpogylon = getPolygonGeoJson(cPolygonList);

            if (isclip) {
                // 读取GeoJSON字符串 返回 SimpleFeatureCollection 要素集合
                polygonCollection = GeoJSONUtil.readGeoJsonByString(geojsonpogylon);
                // 裁剪等值面
                FeatureSource dc = clipFeatureCollection(fc, polygonCollection);
                // 得到多边形GeoJson
                geojsonpogylon = getPolygonGeoJson(dc.getFeatures());
            }
        } catch (Exception e) {
            e.printStackTrace();
        }
        return geojsonpogylon;
    }


    
    public static String calEquiSurfaceLine(double[][] trainData, double[] dataInterval, int[] size, String boundryFile,
                                     boolean isclip) {
        String geojsonline = "";
        try {
            double _undefData = -9999.0;
            SimpleFeatureCollection polylineCollection = null;
            List cPolylineList = new ArrayList();
            List cPolygonList = new ArrayList();

            int width = size[0],
                    height = size[1];
            double[] _X = new double[width];
            double[] _Y = new double[height];

            File file = new File(boundryFile);
            ShapefileDataStore shpDataStore = null;

            shpDataStore = new ShapefileDataStore(file.toURL());
            //设置编码
            Charset charset = Charset.forName("GBK");
            shpDataStore.setCharset(charset);
            String typeName = shpDataStore.getTypeNames()[0];
            SimpleFeatureSource featureSource = null;
            featureSource = shpDataStore.getFeatureSource(typeName);
            SimpleFeatureCollection fc = featureSource.getFeatures();

            double minX = fc.getBounds().getMinX();
            double minY = fc.getBounds().getMinY();
            double maxX = fc.getBounds().getMaxX();
            double maxY = fc.getBounds().getMaxY();

            Interpolate.createGridXY_Num(minX, minY, maxX, maxY, _X, _Y);
            double[][] _gridData = new double[width][height];

            int nc = dataInterval.length;

            _gridData = Interpolate.interpolation_IDW_Neighbor(trainData,
                    _X, _Y, 12, _undefData);// IDW插值

            int[][] S1 = new int[_gridData.length][_gridData[0].length];
            
            List _borders = Contour.tracingBorders(_gridData, _X, _Y,
                    S1, _undefData);

            
            cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc,
                    dataInterval, _undefData, _borders, S1);// 生成等值线

            cPolylineList = Contour.smoothLines(cPolylineList);// 平滑


            geojsonline = getPolylineGeoJson(cPolylineList);

            if (isclip) {
                polylineCollection = GeoJSONUtil.readGeoJsonByString(geojsonline);
                FeatureSource dc = clipFeatureLineCollection(fc, polylineCollection);
                geojsonline = getPolylineGeoJson(dc.getFeatures());
            }

        } catch (Exception e) {
            e.printStackTrace();
        }
        return geojsonline;
    }

    // 等值面
    private static String getPolygonGeoJson(List cPolygonList) {
        String geo = null;
        String geometry = " { "type":"Feature","geometry":";
        String properties = ","properties":{ "hvalue":";

        String head = "{"type": "FeatureCollection"," + ""features": [";
        String end = "  ] }";
        if (cPolygonList == null || cPolygonList.size() == 0) {
            return null;
        }
        try {
            for (Polygon pPolygon : cPolygonList) {

                List