geotools中等值面的生成与OL3中的展示
概述：
本文讲述如何在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<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);
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<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 lv = (Double) contourFeature.getProperty("lvalue")
.getValue();
double hv = (Double) contourFeature.getProperty("hvalue")
.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("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<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();
}
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();
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();
}
}
}2、openlayers3
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<title>Ol3 wms</title>
<link rel="stylesheet" type="text/css" href="../../../plugin/ol3/css/ol.css"/>
<style type="text/css">
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;
}
</style>
<script type="text/javascript" src="../../../plugin/ol3/build/ol-debug.js"></script>
<script type="text/javascript" src="../../../plugin/jquery/jquery-1.8.3.js"></script>
<script type="text/javascript">
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){
color = "245,200,200,255";
}
else if(lvalue>=20&&lvalue<30){
color = "245,183,48,255";
}
else if(lvalue>=25&&lvalue<30){
color = "245,20,48,255";
}
else if(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());
});
}
</script>
</head>
<body onLoad="init()">
<div id="map">
<div id="infobox">
</div>
</div>
</body>
</html>-----------------------------------------------------------------------------------------------
技术博客
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)
转载自：https://blog.csdn.net/GISShiXiSheng/article/details/77687996
