geotools等值线生成

概述

前文中,提到了等值面的生成,后面有人经常会问等值线的生成,本文在前文的基础上做了一点修改,完成了等值线的geotools生成。

效果

实现代码

package com.lzugis.geotools;

import com.amazonaws.util.json.JSONObject;
import com.lzugis.CommonMethod;
import com.lzugis.geotools.utils.FeaureUtil;
import com.lzugis.geotools.utils.GeoJSONUtil;
import com.vividsolutions.jts.geom.Geometry;
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.*;

/**
 * Created by admin on 2017/8/29.
 */
public class EquiSurfaceLine {
    private static String rootPath = System.getProperty("user.dir");

    /**
     * 生成等值面
     *
     * @param trainData    训练数据
     * @param dataInterval 数据间隔
     * @param size         大小,宽,高
     * @param boundryFile  四至
     * @param isclip       是否裁剪
     * @return
     */
    public String calEquiSurface(double[][] trainData,
                                 double[] dataInterval,
                                 int[] size,
                                 String boundryFile,
                                 boolean isclip) {
        String geojsonline = "";
        try {
            double _undefData = -9999.0;
            SimpleFeatureCollection polylineCollection = null;
            List<PolyLine> cPolylineList = new ArrayList<PolyLine>();
            List<Polygon> cPolygonList = new ArrayList<Polygon>();

            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];
            /**
             * double[][] S0,
             * double[] X,
             * double[] Y,
             * int[][] S1,
             * double undefData
             */
            List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y,
                    S1, _undefData);

            /**
             * double[][] S0,
             * double[] X,
             * double[] Y,
             * int nc,
             * double[] contour,
             * double undefData,
             * List<Border> borders,
             * int[][] S1
             */
            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 = clipFeatureCollection(fc, polylineCollection);
                geojsonline = getPolylineGeoJson(dc.getFeatures());
            }

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


    public String getPolylineGeoJson(FeatureCollection fc) {
        FeatureJSON fjson = new FeatureJSON();
        StringBuffer sb = new StringBuffer();
        try {
            sb.append("{\"type\": \"FeatureCollection\",\"features\": ");
            FeatureIterator itertor = fc.features();
            List<String> list = new ArrayList<String>();
            while (itertor.hasNext()) {
                SimpleFeature feature = (SimpleFeature) itertor.next();
                StringWriter writer = new StringWriter();
                fjson.writeFeature(feature, writer);
                list.add(writer.toString());
            }
            itertor.close();
            sb.append(list.toString());
            sb.append("}");
        } catch (Exception e) {
            e.printStackTrace();
        }
        return sb.toString();
    }

    public FeatureSource clipFeatureCollection(FeatureCollection fc,
                                               SimpleFeatureCollection gs) {
        FeatureSource 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();
                Geometry dataGeometry = (Geometry) dataFeature.getProperty(
                        "the_geom").getValue();
                while (contourFeatureIterator.hasNext()) {
                    Feature contourFeature = contourFeatureIterator.next();
                    Geometry contourGeometry = (Geometry) contourFeature
                            .getProperty("geometry").getValue();
                    double v = (Double) contourFeature.getProperty("value")
                            .getValue();
                    if (dataGeometry.intersects(contourGeometry)) {
                        Geometry geo = dataGeometry
                                .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();

            SimpleFeatureCollection sc = FeaureUtil
                    .creatSimpleFeatureByFeilds(
                            "polygons",
                            "crs:4326,the_geom:LineString,value:double",
                            values);
            cs = FeaureUtil.creatFeatureSourceByCollection(sc);

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

        return cs;
    }


    public String getPolylineGeoJson(List<PolyLine> cPolylineList) {
        String geo = null;
        String geometry = " { \"type\":\"Feature\",\"geometry\":";
        String properties = ",\"properties\":{ \"value\":";

        String head = "{\"type\": \"FeatureCollection\"," + "\"features\": [";
        String end = "  ] }";
        if (cPolylineList == null || cPolylineList.size() == 0) {
            return null;
        }
        try {
            for (PolyLine pPolyline : cPolylineList) {
                List<Object> ptsTotal = new ArrayList<Object>();

                for (PointD ptD : pPolyline.PointList) {
                    List<Double> pt = new ArrayList<Double>();
                    pt.add(ptD.X);
                    pt.add(ptD.Y);
                    ptsTotal.add(pt);
                }

                JSONObject js = new JSONObject();
                js.put("type", "LineString");
                js.put("coordinates", ptsTotal);

                geo = geometry + js.toString() + properties + pPolyline.Value + "} }" + "," + geo;
            }
            if (geo.contains(",")) {
                geo = geo.substring(0, geo.lastIndexOf(","));
            }

            geo = head + geo + end;
        } catch (Exception e) {
            e.printStackTrace();
            return geo;
        }
        return geo;
    }

    public static void main(String[] args) {
        long start = System.currentTimeMillis();

        EquiSurfaceLine equiSurface = new EquiSurfaceLine();
        CommonMethod cm = new CommonMethod();

        double[] bounds = {110.759, 31.23, 113.112, 32.6299};

        double[][] trainData = new double[3][100];

        for (int i = 0; i < 100; i++) {
            double x = bounds[0] + new Random().nextDouble() * (bounds[2] - bounds[0]),
                    y = bounds[1] + new Random().nextDouble() * (bounds[3] - bounds[1]),
                    v = 0 + new Random().nextDouble() * (45 - 0);
            trainData[0][i] = x;
            trainData[1][i] = y;
            trainData[2][i] = v;
        }

        double[] dataInterval = new double[]{20, 25, 30, 35, 40, 45};

        String boundryFile = rootPath + "/data/shp/XY_01_XZQH_PY.shp";

        int[] size = new int[]{100, 100};

        boolean isclip = true;

        try {
            String strJson = equiSurface.calEquiSurface(trainData, dataInterval, size, boundryFile, isclip);
            String strFile = rootPath + "/out/XY_01_XZQH_PY_line1.geojson";
            cm.append2File(strFile, strJson);
            System.out.println(strFile + "差值成功, 共耗时" + (System.currentTimeMillis() - start) + "ms");
        } catch (Exception e) {
            e.printStackTrace();
        }
    }
}

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏闵开慧

曾经做过的40道程序设计课后习题总结(一)

曾经做过的40道程序设计课后习题总结(一) 课后习题目录 1 斐波那契数列 2 判断素数 3 水仙花数 4 分解质因数 5 杨辉三角 6 学习成绩查询...

3408
来自专栏人工智能LeadAI

讨厌算法的程序员 | 第五章 合并算法

本篇介绍的“合并”算法,是为后面学习“归并排序”的一个准备。合并算法是归并排序中的一个子算法,请注意两者之间的关系和差异。 之所以把它独立成一篇,一方面是一旦了...

3645
来自专栏小詹同学

Leetcode打卡 | No.015 三数之和

欢迎和小詹一起定期刷leetcode,每周一和周五更新一题,每一题都吃透,欢迎一题多解,寻找最优解!这个记录帖哪怕只有一个读者,小詹也会坚持刷下去的!

1202
来自专栏海天一树

小朋友学C语言(4):单精度浮点数与双精度浮点数

上节课 简单介绍了浮点数。计算机程序中的浮点数分为单精度浮点数和双精度浮点数。 单精度和双精度精确的范围不一样。 计算机里的最基本的存储单位用位(bit)来表...

37812
来自专栏小文博客

十进制转换二进制(C语言)

2932
来自专栏数据结构与算法

BZOJ 1174: [Balkan2007]Toponyms

Description 给你一个字符集合,你从其中找出一些字符串出来. 希望你找出来的这些字符串的最长公共前缀*字符串的总个数最大化. Input 第一行给出...

2866
来自专栏来自地球男人的部落格

[LeetCode] 119. Pascal's Triangle II

【原题】 Given an index k, return the kth row of the Pascal’s triangle. For exampl...

2026
来自专栏ml

位运算的方法,小结

文章来源未知----再次声明为转载... 本文是针对使用位运算来实现一些方法,我们都知道位运算的代价比其他符号运算都低,所以当一个方法只使用位运算且运算次数与其...

37213
来自专栏数据结构与算法

P3808 【模版】AC自动机(简单版)

题目背景 这是一道简单的AC自动机模版题。 用于检测正确性以及算法常数。 为了防止卡OJ,在保证正确的基础上只有两组数据,请不要恶意提交。 题目描述 给定n个模...

2726
来自专栏苦逼的码农

从0打卡leetcode之day10--整数回文数

判断一个整数是否是回文数。回文数是指正序(从左向右)和倒序(从右向左)读都是一样的整数。

712

扫码关注云+社区