geotools中等值面的生成与OL3中的展示

综合技术 2017-09-10

概述:

本文讲述如何在geotools中IDW插值生成等值面,并根据给定shp进行裁剪,并生成geojson数据,以及Openlayers3中展示。

效果:

插值数据

裁剪结果

裁剪区域数据

实现代码:

1、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 EquiSurface {
    /**
     * 生成等值面
     *
     * @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 geojsonpogylon = "";
        try {
            double _undefData = -9999.0;
            SimpleFeatureCollection polygonCollection = 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);// 平滑

            cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList,
                    _borders, dataInterval);

            geojsonpogylon = getPolygonGeoJson(cPolygonList);

            if (isclip) {
                polygonCollection = GeoJSONUtil.readGeoJsonByString(geojsonpogylon);
                FeatureSource dc = clipFeatureCollection(fc, polygonCollection);
                geojsonpogylon = getPolygonGeoJson(dc.getFeatures());
            }
        } catch (Exception e) {
            e.printStackTrace();
        }
        return geojsonpogylon;
    }

    private FeatureSource clipFeatureCollection(FeatureCollection fc,
                                                SimpleFeatureCollection gs) {
        FeatureSource cs = null;
        try {
            List<map> values = new ArrayList<map>();
            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 lv = (Double) contourFeature.getProperty("lvalue")
                            .getValue();
                    double hv = (Double) contourFeature.getProperty("hvalue")
                            .getValue();
                    if (dataGeometry.intersects(contourGeometry)) {
                        Geometry geo = dataGeometry
                                .intersection(contourGeometry);
                        Map map = new HashMap();
                        map.put("the_geom", geo);
                        map.put("lvalue", lv);
                        map.put("hvalue", hv);
                        values.add(map);
                    }

                }

            }

            contourFeatureIterator.close();
            dataFeatureIterator.close();

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

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

        return cs;
    }

    private String getPolygonGeoJson(FeatureCollection fc) {
        FeatureJSON fjson = new FeatureJSON();
        StringBuffer sb = new StringBuffer();
        try {
            sb.append("{"type": "FeatureCollection","features": ");
            FeatureIterator itertor = fc.features();
            List list = new ArrayList();
            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();
    }

    private 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 ptsTotal = new ArrayList();
                List pts = new ArrayList();

                PolyLine pline = pPolygon.OutLine;

                for (PointD ptD : pline.PointList) {
                    List pt = new ArrayList();
                    pt.add(ptD.X);
                    pt.add(ptD.Y);
                    pts.add(pt);
                }

                ptsTotal.add(pts);

                if (pPolygon.HasHoles()) {
                    for (PolyLine cptLine : pPolygon.HoleLines) {
                        List cpts = new ArrayList();
                        for (PointD ccptD : cptLine.PointList) {
                            List pt = new ArrayList();
                            pt.add(ccptD.X);
                            pt.add(ccptD.Y);
                            cpts.add(pt);
                        }
                        if (cpts.size() > 0) {
                            ptsTotal.add(cpts);
                        }
                    }
                }

                JSONObject js = new JSONObject();
                js.put("type", "Polygon");
                js.put("coordinates", ptsTotal);
                double hv = pPolygon.HighValue;
                double lv = pPolygon.LowValue;

                if (hv == lv) {
                    if (pPolygon.IsClockWise) {
                        if (!pPolygon.IsHighCenter) {
                            hv = hv - 0.1;
                            lv = lv - 0.1;
                        }

                    } else {
                        if (!pPolygon.IsHighCenter) {
                            hv = hv - 0.1;
                            lv = lv - 0.1;
                        }
                    }
                } else {
                    if (!pPolygon.IsClockWise) {
                        lv = lv + 0.1;
                    } else {
                        if (pPolygon.IsHighCenter) {
                            hv = hv - 0.1;
                        }
                    }

                }

                geo = geometry + js.toString() + properties + hv
                        + ", "lvalue":" + lv + "} }" + "," + 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();

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

        double[] bounds = {73.4510046356223, 18.1632471876417,
                134.976797646506, 53.5319431522236};

        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 = "D:\data\国家基础地理信息系统SHP文件\国界\国界\bou1_4p.shp";

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

        boolean isclip = true;

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

2、openlayers3


	
	Ol3 wms
	
	
		body, #map {
			border: 0px;
			margin: 0px;
			padding: 0px;
			width: 100%;
			height: 100%;
			font-size: 13px;
		}
		#infobox{
			position: absolute;
			top: 10px;
			right: 10px;
			z-index: 999;
			background: #ffffff;
			border: 1px solid #000;
			box-shadow: 1px 1px 1px #C0C0C0;
		}
	
	
	
	
		function init(){
			var bounds = [73.4510046356223, 18.1632471876417,
				134.976797646506, 53.5319431522236];
			var projection = new ol.proj.Projection({
				code: 'EPSG:4326',
				units: 'degrees'
			});
			$.get("data/china1.json",null,function(result){
				var features = (new ol.format.GeoJSON()).readFeatures(result);
				var vectorSource = new ol.source.Vector();
				vectorSource.addFeatures(features);
				var vector = new ol.layer.Vector({
					source: vectorSource,
					style:function(feature, resolution) {
						var lvalue = feature.get("lvalue"), color = "0,0,0,0";
						if(lvalue=20&&lvalue=25&&lvalue=30&&lvalue<35){
							color = "0,183,48,255";
						}
						else{
							color = "0,255,255,255";
						}
						return new ol.style.Style({
						    stroke: new ol.style.Stroke({
						        color: '#ffffff',
//						        lineDash: [10],
						        width: 1
						    }),
						    fill: new ol.style.Fill({
						        color: "rgba("+color+")",
						        opacity:0.6
						    })
						});
					}
				});
                var tiled = new ol.layer.Tile({
                    source: new ol.source.TileWMS({
                        url: 'http://192.168.10.185:8086/geoserver/lzugis/wms',
                        params: {'FORMAT': 'image/png',
                            'VERSION': '1.1.1',
                            tiled: true,
                            LAYERS: 'lzugis:capital',
                            STYLES: ''
                        },
                        serverType: 'geoserver'
                    })
                });
				var map = new ol.Map({
					controls: ol.control.defaults({
						attribution: false
					}),
					target: 'map',
					layers: [vector,tiled],
					view: new ol.View({
						projection: projection
					})
				});
				map.getView().fit(bounds, map.getSize());
			});
		}
	


-----------------------------------------------------------------------------------------------

技术博客

CSDN:http://blog.csdn.NET/gisshixisheng

博客园:http://www.cnblogs.com/lzugis/

在线教程

http://edu.csdn.Net/course/detail/799

Github

https://github.com/lzugis/

联系方式

q q:1004740957

e-mail:niujp08@qq.com

公众号:lzugis15

Q Q 群:452117357(webgis)

337469080(Android)

CSDN博客

责编内容by:CSDN博客 (源链)。感谢您的支持!