【GIS小案例】基于Java + Geotools + wContour的等值面图

1 引言

展示雨量数据的分布情况，等值面图是非常合适的方法。以锦州市某日的降雨量为例，雨量站点大概有100个，分布于整个锦州市。首先尝试的方法是利用kriging.js生成等值面，但是由于锦州市边界区域太大，浏览器客户端生成压力非常大，速度也很慢。所以尝试着在服务端生成，然后客户端调用显示。

2 等值面生成的两种方式

2.1 服务端生成

通用的方式是基于Java + Geotools + wContour在服务端处理数据，适用业务如下：

大范围高密度的空间数据插值
服务器端可以设置定时任务，流水线处理原始数据，生成目标数据（图片或者GeoJSON数据）

2.2 客户端生成

通用的方式是前端插值处理数据，生成网格数据，再进行等值线或等值面的生成，适用业务如下：

数据更新频率高，实时性强
按需服务，减轻服务器压力

可以使用如下方式：

使用turf.js内置函数通过离散点插值、等值面绘制、裁剪，生成矢量数据然后绘制数据
使用kriging.js内置的函数，离散点生成网格数据然后网格渲染
结合turf.js和kriging.js两者的优点，提升交互性能、插值效果和显示效果

3 基于Java + Geotools + wContour在服务端处理数据

生成结果图1（边界不裁剪）：

生成结果图2（边界裁剪）：

完整代码：

package com.example.demo.meteoInfoTest;import cn.hutool.json.JSONArray;
import com.tt.CommonMethod;
import com.tt.GeoJSONUtil;
import org.geotools.data.DataUtilities;
import org.geotools.data.collection.ListFeatureCollection;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.geotools.feature.SchemaException;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.geojson.feature.FeatureJSON;
import org.json.simple.JSONObject;
import org.locationtech.jts.geom.Geometry;
import org.opengis.feature.simple.SimpleFeature;import org.opengis.feature.simple.SimpleFeatureType;
import wcontour.Contour;
import wcontour.Interpolate;
import wcontour.global.Border;
import wcontour.global.PointD;
import wcontour.global.PolyLine;
import wcontour.global.Polygon;import java.io.File;
import java.io.IOException;
import java.io.StringWriter;
import java.nio.charset.Charset;
import java.util.*;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<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];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);geojsonpogylon = getPolygonGeoJson(cPolygonList);if (isclip) {polygonCollection = GeoJSONUtil.readGeoJsonByString(geojsonpogylon);SimpleFeatureCollection sm = clipPolygonFeatureCollection(fc, polygonCollection);FeatureCollection featureCollection = sm;FeatureJSON featureJSON = new FeatureJSON();StringWriter writer = new StringWriter();featureJSON.writeFeatureCollection(featureCollection, writer);geojsonpogylon = writer.toString();}} catch (Exception e) {e.printStackTrace();}return geojsonpogylon;}/*** 裁剪图形* @param fc* @return*/private SimpleFeatureCollection clipPolygonFeatureCollection(FeatureCollection fc,SimpleFeatureCollection gs) throws SchemaException, IOException {SimpleFeatureCollection simpleFeatureCollection = null;SimpleFeatureType TYPE = DataUtilities.createType("polygons","the_geom:MultiPolygon,lvalue:double,hvalue:double");List<SimpleFeature> list = new ArrayList<>();SimpleFeatureIterator contourFeatureIterator = gs.features();FeatureIterator dataFeatureIterator = fc.features();while (dataFeatureIterator.hasNext()){SimpleFeature dataFeature = (SimpleFeature) dataFeatureIterator.next();Geometry dataGeometry = (Geometry) dataFeature.getDefaultGeometry();contourFeatureIterator = gs.features();while (contourFeatureIterator.hasNext()){SimpleFeature contourFeature = contourFeatureIterator.next();Geometry contourGeometry = (Geometry) contourFeature.getDefaultGeometry();Double lv = (Double) contourFeature.getProperty("lvalue").getValue();Double hv = (Double) contourFeature.getProperty("hvalue").getValue();if (dataGeometry.getGeometryType() == "MultiPolygon"){for (int i = 0; i <dataGeometry.getNumGeometries(); i++){Geometry geom = dataGeometry.getGeometryN(i);if (geom.intersects(contourGeometry)){Geometry geo = geom.intersection(contourGeometry);SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);featureBuilder.add(geo);featureBuilder.add(lv);featureBuilder.add(hv);SimpleFeature feature = featureBuilder.buildFeature(null);list.add(feature);}}} else {if (dataGeometry.intersects(contourGeometry)){Geometry geo = dataGeometry.intersection(contourGeometry);SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(TYPE);featureBuilder.add(geo);featureBuilder.add(lv);featureBuilder.add(hv);SimpleFeature feature = featureBuilder.buildFeature(null);list.add(feature);}}}}contourFeatureIterator.close();dataFeatureIterator.close();simpleFeatureCollection = new ListFeatureCollection(TYPE, list);return simpleFeatureCollection;}private String getPolygonGeoJson(List<Polygon> 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<Object> ptsTotal = new ArrayList<Object>();List<Object> pts = new ArrayList<Object>();PolyLine pline = pPolygon.OutLine;for (PointD ptD : pline.PointList) {List<Double> pt = new ArrayList<Double>();pt.add(ptD.X);pt.add(ptD.Y);pts.add(pt);}ptsTotal.add(pts);if (pPolygon.HasHoles()) {for (PolyLine cptLine : pPolygon.HoleLines) {List<Object> cpts = new ArrayList<Object>();for (PointD ccptD : cptLine.PointList) {List<Double> pt = new ArrayList<Double>();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();// 雨量站点经纬度及雨量值JSONArray jsonArray = GeoJSONUtil.yuliang();double[][] trainData = new double[jsonArray.size()][3];for (int i = 0; i < jsonArray.size(); i++){trainData[i][0] = jsonArray.getJSONObject(i).getDouble("LGTD");trainData[i][1] = jsonArray.getJSONObject(i).getDouble("LTTD");trainData[i][2] = jsonArray.getJSONObject(i).getDouble("DYP");}double[] dataInterval = new double[]{0.0001, 2, 4, 6, 9, 10, 20, 30};String boundryFile = "F:\\我的文档\\等值面图研究\\jz-data\\jz-shape\\锦州_R.shp";int[] size = new int[]{100, 100};boolean isclip = false;try {String strJson = equiSurface.calEquiSurface(trainData, dataInterval, size, boundryFile, isclip);String strFile = "C:\\Users\\admin\\Desktop\\临时\\jinzhou.json";cm.append2File(strFile, strJson);System.out.println(strFile + "差值成功, 共耗时" + (System.currentTimeMillis() - start) + "ms");} catch (Exception e) {e.printStackTrace();}}
}

    public static SimpleFeatureCollection readGeoJsonByString(String geojsonpogylon) throws IOException {FeatureJSON fjson = new FeatureJSON();SimpleFeatureCollection featureCollection = (SimpleFeatureCollection) fjson.readFeatureCollection(geojsonpogylon);return featureCollection;}public void append2File(String strFile, String strJson){File f=new File(strFile);//新建一个文件对象，如果不存在则创建一个该文件FileWriter fw;try {fw = new FileWriter(f);fw.write(strJson);//将字符串写入到指定的路径下的文件中fw.close();} catch (IOException e) { e.printStackTrace(); }}

<script type="text/javascript">var bounds = [120.711715945734340, 40.783993710141660,122.595342945920280, 42.126609285989790];var projection = new ol.proj.Projection({code: 'EPSG:4326',units: 'degrees'});$.get("../data/jinzhou/jinzhou-cj2.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<4){color = "166,242,143,255";}else if(lvalue>=4&&lvalue<6){color = "87,238,40,255";}else if(lvalue>=6&&lvalue<10){color = "22,210,22,255";}else if(lvalue>=10&&lvalue<20){color = "97,184,255,255";}else{color = "28,125,204,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 map = new ol.Map({controls: ol.control.defaults({attribution: false}),target: 'map',layers: [vector],view: new ol.View({projection: projection})});map.getView().fit(bounds, map.getSize());});
</script>

参考文献:

[1] 等值线图以及色斑图实现方式综述（https://www.jianshu.com/p/04362ff87641）

[2] geotools中等值面的生成与OL3中的展示_LZUGIS-CSDN博客_geoserver 等值线

[3] java实现反距离权重插值算法生成geojson矢量数据_兴诚的博客-CSDN博客_geotools

[4] 2021-10-28 实现等值面的一些问题 - 简书

【GIS小案例】基于Java + Geotools + wContour的等值面图相关推荐

【GIS小案例】台风烟花的轨迹动图
效果如图所示: 代码: <!DOCTYPE html> <html lang="en"> <head><meta charset=&quo ...
微信小程序基于java实现v2支付，提现，退款
v2支付 v2微信官方文档封装支付请求实体 import io.swagger.annotations.ApiModelProperty; import lombok.Getter; import ...
英雄联盟的小案例理解Java中如何实现OCP原则
案例: 英雄联盟的英雄.道具.地图,每年都会进行频繁变更如果没有使用软件工程的开发思想,随便改其中一个道具的属性,就可能会导致非常严重的错误要实现变更/增加英雄时,可选英雄数量和玩家开始一局游戏时 ...
【GIS小案例】点聚合效果的实现
1,效果图 2,实现代码 <script>var viewer = new Cesium.Viewer('cesiumContainer');viewer.scene.open(" ...
【GIS小案例】CesiumHeatmap热力图
1,CesiumHeatmap: GitHub - danwild/CesiumHeatmap: A library to add heatmaps (using heatmap.js) to the ...
Opencv4学习-2、小案例之绿布抠图-视频背景图替换
一.绿布抠图-背景图替换主要是基于图像色彩空间,转换为HSV色彩空间实现mask层提取.然后通过一些简单的图像处理比如形态学开闭操作.高斯模糊等到完整的mask区域,运用mask区域生成权重系数,对 ...
基于java的圆通快递单号自动识别api接口代码实例
一.产品介绍快递单号识别,输入运单号自动识别物流公司,实时返回对应物流公司编码.查询单号时,返回的结果可能存在一个或多个物流公司编码,快递鸟大数据平台通过智能分析,实时更新单号库,保障物流公司编码准 ...
基于Java小案例家庭收入支出记录
基于Java小案例家庭收入支出记录跟着视频来写的,不喜勿喷,谢谢 package JavaText;public class FamilyAccount {public static void ma ...
Java Web应用小案例：查询城市天气信息
Java Web应用小案例:查询城市天气信息本期上大数据1班动态网站设计与开发课,经过半期的学习,学生已经可以利用所学的JSP知识开发简单的基于后台数据库操作的动态网站,但是这是远远不够的,课程教学 ...
猜物品游戏java编程_小猿圈Java初学者练习小案例：猜数字游戏
对于Java初学者,如果没有好的引导,可能会觉得自己学什么都不好,学什么都不会,这个时候就要给他们一下小的案例,让他们去实践一下,让他们知道自己学的东西是可以用到的,小猿圈java讲师为你准备了Jav ...

【GIS小案例】基于Java + Geotools + wContour的等值面图

【GIS小案例】基于Java + Geotools + wContour的等值面图相关推荐

最新文章

热门文章