6.gis空间插值之java opengis+goetools+wcontour

需要引入的gis库：org.geotools、wcontour、org.opengis,直接上代码:

   //创建测试数据
    public List<InterpolData> createInterpolData() {
        List<InterpolData> listData = new ArrayList<>();

        double lon = 105.99;
        double lat = 32.148;
        double pValue = 60;
        for (int i = 0; i < 50; i++) {
            InterpolData interp = new InterpolData();
            lon += Math.random();
            lat += Math.random();
            pValue += Math.random() * 50;
            interp.setLon(lon);
            interp.setLat(lat);
            interp.setPValue(pValue);
            listData.add(interp);
        }
        return listData;
    }


//生成等值线或等值面
    public FeatureCollection createequiSurface() {
        SimpleFeatureCollection fc = createFeatureCl('');
        Object geom = fc.features().next().getAttribute("the_geom");
        WKTReader2 wktReader = new WKTReader2();
        Geometry oldGeom = wktReader.read(geom.toString());
        FeatureCollection pfeatureCol = null;
        double _undefData = -9999.0;
        List<PolyLine> cPolylineList;
        List<Polygon> cPolygonList;
        int width = 100, height = 100;
        double[] _X = new double[width];
        double[] _Y = new double[height];
        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;
        int nc =12;
        double[] dataInterval = {62, 65, 72, 75, 78, 81, 85, 90, 96, 100, 105, 110.9};

        double[][] trainData = trainData(createInterpolData());
        // IDW插值格网点
        _gridData = Interpolate.interpolation_IDW_Neighbor(trainData,
                _X, _Y, 12, _undefData);
        int[][] S1 = new int[_gridData.length][_gridData[0].length];
        // 获取轮廓
        List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y,
                S1, _undefData);
        // 生成等值线
        cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc,
                dataInterval, _undefData, _borders, S1);
        // 平滑
        cPolylineList = Contour.smoothLines(cPolylineList);
       //等值面
            cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList,
                    _borders, dataInterval);
            pfeatureCol = getFeatureCollection(cPolygonList);
        //等值线
           // pfeatureCol = getFeaturesBycPolyline(cPolylineList);
        
     if (inParam.getClip()) {
            pfeatureCol = clipFeatureCollection(fc, (SimpleFeatureCollection) pfeatureCol);
            pfeatureCol = new ClippedFeatureCollection((SimpleFeatureCollection) pfeatureCol, oldGeom, true);
        }

        return pfeatureCol;
    }
    /**
     * 将等值线转换成FeatureCollection
     */
    public static FeatureCollection getFeaturesBycPolyline(List<PolyLine> cPolylineList) throws SchemaException {
        if (cPolylineList == null || cPolylineList.size() == 0) {
            return null;
        }
        FeatureCollection cs = null;
        List<Map<String, Object>> values = new ArrayList<Map<String, Object>>();
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
        LineString line;
        for (PolyLine pPolyline : cPolylineList) {
            Coordinate[] coordinates = new Coordinate[pPolyline.PointList.size()];
            for (int i = 0, len = pPolyline.PointList.size(); i < len; i++) {
                PointD ptd = pPolyline.PointList.get(i);
                coordinates[i] = new Coordinate(ptd.X, ptd.Y);
            }

            line = geometryFactory.createLineString(coordinates);

            Map<String, Object> map = new HashMap<String, Object>();
            map.put("the_geom", line);
            map.put("value", pPolyline.Value);
            values.add(map);
        }
        cs = creatFeatureCollection(
                "polygons",
                "the_geom:LineString,value:double",
                values);
        return cs;
    }

    /**
     * 将生成的等值面转换成FeatureCollection
     */
    public static FeatureCollection getFeatureCollection(List<Polygon> cPolygonList) {

        if (cPolygonList == null || cPolygonList.size() == 0) {
            return null;
        }
        FeatureCollection cs = null;
        List<Map<String, Object>> values = new ArrayList<Map<String, Object>>();
        try {
            for (Polygon pPolygon : cPolygonList) {
                //外圈
                LinearRing mainRing;
                Coordinate[] coordinates = new Coordinate[pPolygon.OutLine.PointList.size()];
                for (int i = 0, len = pPolygon.OutLine.PointList.size(); i < len; i++) {
                    PointD ptd = pPolygon.OutLine.PointList.get(i);
                    coordinates[i] = new Coordinate(ptd.X, ptd.Y);
                }
                mainRing = geometryFactory.createLinearRing(coordinates);
                //孔洞
                LinearRing[] holeRing = new LinearRing[pPolygon.HoleLines.size()];
                for (int i = 0; i < pPolygon.HoleLines.size(); i++) {
                    PolyLine hole = pPolygon.HoleLines.get(i);
                    Coordinate[] coordinates_h = new Coordinate[hole.PointList.size()];
                    for (int j = 0, len = hole.PointList.size(); j < len; j++) {
                        PointD ptd = hole.PointList.get(j);
                        coordinates_h[j] = new Coordinate(ptd.X, ptd.Y);
                    }
                    holeRing[i] = geometryFactory.createLinearRing(coordinates_h);
                }
                org.locationtech.jts.geom.Polygon geo = geometryFactory.createPolygon(mainRing, holeRing);
                Map<String, Object> map = new HashMap<String, Object>();
                map.put("the_geom", geo);
                map.put("value", pPolygon.LowValue);
                values.add(map);
            }

            cs = creatFeatureCollection(
                    "polygons",
                    "the_geom:Polygon,value:double",
                    values);
        } catch (Exception e) {
            e.printStackTrace();
        }
        return cs;
    }

    /**
     *
     * 构建特定的FeatureCollection
     */
    public static FeatureCollection creatFeatureCollection(String typeName, String typeSpec, List<Map<String, Object>> values) throws SchemaException, SchemaException {
        SimpleFeatureType type = DataUtilities.createType(typeName, typeSpec);
        SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(type);
        DefaultFeatureCollection collection = new DefaultFeatureCollection();
        for (Map feat : values) {
            featureBuilder.reset();
            featureBuilder.add(feat.get("the_geom"));
            System.out.println(feat.get("value"));
            featureBuilder.add(feat.get("value"));
            SimpleFeature feature = featureBuilder.buildFeature(null);
            collection.add(feature);
        }
        return collection;
    }

    /**
     * 
     * 两要素进行裁剪
     */
    public FeatureCollection clipFeatureCollection(FeatureCollection fc, SimpleFeatureCollection gs) {

        FeatureCollection cs = null;
        try {
            List<Map<String, Object>> values = new ArrayList<Map<String, Object>>();
            FeatureIterator contourFeatureIterator = gs.features();
            FeatureIterator dataFeatureIterator = fc.features();
            while (dataFeatureIterator.hasNext()) {
                Feature dataFeature = dataFeatureIterator.next();
                Object dataGeometry = dataFeature.getProperty(
                        "the_geom").getValue();
                if (dataGeometry instanceof MultiPolygon) {
                    MultiPolygon p = (MultiPolygon) dataGeometry;
                    while (contourFeatureIterator.hasNext()) {
                        Feature contourFeature = contourFeatureIterator.next();
                        Geometry contourGeometry = (Geometry) contourFeature
                                .getProperty("the_geom").getValue();
                        double v = (Double) contourFeature.getProperty("value")
                                .getValue();
                        if (p.intersects(contourGeometry)) {
                            Geometry geo = p.intersection(contourGeometry);
                            Map<String, Object> map = new HashMap<String, Object>();
                            map.put("the_geom", geo);
                            map.put("value", v);
                            values.add(map);
                        }
                    }
                } else {
                    Geometry p = (Geometry) dataGeometry;
                    while (contourFeatureIterator.hasNext()) {
                        Feature contourFeature = contourFeatureIterator.next();
                        Geometry contourGeometry = (Geometry) contourFeature
                                .getProperty("the_geom").getValue();
                        double v = (Double) contourFeature.getProperty("value")
                                .getValue();
                        if (p.intersects(contourGeometry)) {
                            Geometry geo = p.intersection(contourGeometry);
                            Map<String, Object> map = new HashMap<String, Object>();
                            map.put("the_geom", geo);
                            map.put("value", v);
                            values.add(map);
                        }
                    }
                }
            }
            contourFeatureIterator.close();
            dataFeatureIterator.close();
            cs = creatFeatureCollection(
                    "MultiPolygons",
                    "geom:MultiPolygon:srid=4326,value:double",
                    values);

        } catch (SchemaException e) {
            e.printStackTrace();
        }
        return cs;
    }

    /**
     * 根据条件获取裁剪数据的要素
     */
    public SimpleFeatureCollection createFeatureCl(String stWhere) throws IOException, CQLException {
        ShapefileDataStore store = new ShapefileDataStore(new File(filePath)
                .toURI().toURL());
        store.setCharset(Charset.forName(encode));
        FeatureSource<SimpleFeatureType, SimpleFeature> features = store
                .getFeatureSource(store.getTypeNames()[0]);
        if (!StringUtils.isEmpty(stWhere)) {
            Filter filter = ff.equals(ff.property("NAME"), ff.literal(stWhere));
            return (SimpleFeatureCollection) features.getFeatures(filter);
        }
        return (SimpleFeatureCollection) features.getFeatures();
    }
   /**
     * 转换插值数据类型
     */
    public double[][] trainData(List<InterpolData> inData) {
        double[][] douData = new double[inData.size()][3];
        for (int i = 0; i < inData.size(); i++) {
            douData[i][0] = inData.get(i).getLon();
            douData[i][1] = inData.get(i).getLat();
            douData[i][2] = inData.get(i).getPValue();
        }
        return douData;
    }