最近需要在用java来分析数据生成等值面图片，在网上搜索各路资料后自己再修修改改搞定了。
在这里总结分享一下。
实现思路主要是 wContour分析等值数据，geotools用于转换分析结果和边界裁切。
分析部分参考了java实现NC数据等值线等值面可视化
用到的第三方包：wContour(点击下载)，geotools (从官方下载或通过maven安装一下)
没有积分的话可以到百度云下载：
wContour(百度云)——提取码：22wb
等值数据分析
private static GeometryFactory geometryFactory = new GeometryFactory();
//生成等值线数据，
/**
* 生成等值面图片
* @param trainData 等值点数据，3*n 的二维数组,按顺序是"x"组，"y"组,"值"组
* @param dataInterval 等值数组，会按照这个参数的值生成对应的等值线
* @param size 等值面插值点数量,越多越精细也越慢
* @param boundaryFile 等值面裁切shp文件路径
*/
public FeatureCollection equiSurface(double[][] trainData,
double[] dataInterval,
int[] size,
String boundryFile) throws IOException {
double _undefData = -9999.0;
List<PolyLine> cPolylineList;
List<Polygon> cPolygonList;
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.toURI().toURL());
//设置编码，根据shp文件设置，一般中文操作系统的arcgis导出的是GBK，也可能是utf-8或其它，自行测试
Charset charset = Charset.forName("GBK");
shpDataStore.setCharset(charset);
String typeName = shpDataStore.getTypeNames()[0];
SimpleFeatureSource featureSource = null;
featureSource = shpDataStore.getFeatureSource(typeName);
SimpleFeatureCollection fc = featureSource.getFeatures();
//获取shp边界，生成格网点
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 = dataInterval.length;
// 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);
// 等值面结果转换和裁切
FeatureCollection polygonCollection = getFeatureCollection(cPolygonList);
FeatureCollection source = clipFeatureCollection(fc, (SimpleFeatureCollection) polygonCollection);
// 结果生成geojson
// FeatureJSON fjson = new FeatureJSON();
// StringWriter writer = new StringWriter();
// fjson.writeFeatureCollection(source, writer);
// String json = writer.toString();
return source;结果转换
从 wContour 的 polygon 转换为 geotools 的 polygon ，以便后续裁切工作。
//结果转换
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 = FeatureUtil.creatFeatureCollection(
"polygons",
"the_geom:Polygon:srid=4326,value:double",
values);
} catch (Exception e) {
e.printStackTrace();
}
return cs;
}结果裁切
核心方法是 intersects 及 intersection，需要注意一下几何类型，如果是MultiPolygon的话不能转Geometry，因此做个判断。
//结果裁切
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 = FeatureUtil.creatFeatureCollection(
"MultiPolygons",
"the_geom:MultiPolygon:srid=4326,value:double",
values);
}catch (SchemaException e) {
e.printStackTrace();
}
return cs;
}用于创建FeatureCollection的工具方法
//geotools创建FeatureCollection
public static FeatureCollection creatFeatureCollection(String typeName, String typeSpec, List<Map<String, Object>> values) throws 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"));
featureBuilder.add(feat.get("value"));
SimpleFeature feature = featureBuilder.buildFeature(null);
collection.add(feature);
}
return collection;
}出图
出图的部分就是将等值面结果放入MapContent中并添加颜色样式，生成图片。
Shape2Image shp2img = new Shape2Image();
//添加图层
shp2img.addShapeLayer(fs,levelProps,OPACITY);
//输出图片
shp2img.getMapContent(params, outPath);private static MapContent map = new MapContent();
/**
* 添加featureCollection等值面图层
*
* @param featureCollection 等值面要素几何
* @param levelProps 色阶,结构如：{0.1:"#a5f38d"}
* @param opacity 透明度
*/
public void addShapeLayer(FeatureCollection featureCollection, Map<Double,String> levelProps,float opacity) {
try {
// 由坐标顺序引发坐标变换，这三行由于修正数据，不加的话会出现要素漏缺。
SimpleFeatureType simpleFeatureType = (SimpleFeatureType) featureCollection.getSchema();
String crs = CRS.lookupIdentifier(simpleFeatureType.getCoordinateReferenceSystem(), true);
featureCollection = new ForceCoordinateSystemFeatureResults(featureCollection,
CRS.decode(crs, true));
//创建样式
StyleFactory sf = new StyleFactoryImpl();
FilterFactory ff = new FilterFactoryImpl();
FeatureTypeStyle fts = sf.createFeatureTypeStyle();
for (Map.Entry entry:levelProps.entrySet()) {
double key = (Double) entry.getKey();
String value = (String) entry.getValue();
Fill fill = sf.createFill(ff.literal(value),ff.literal(opacity));
Stroke stroke = sf.createStroke(ff.literal("#ffffff"),ff.literal(0),ff.literal(0));
Symbolizer symbolizer = sf.createPolygonSymbolizer(stroke, fill, "the_geom");
Rule rule = sf.createRule();
rule.setName("dzm_"+key);
rule.symbolizers().add(symbolizer);
Filter filter = ECQL.toFilter("value="+key);
rule.setFilter(filter);
fts.rules().add(rule);
}
Style style = sf.createStyle();
style.setName("style_dzm");
style.featureTypeStyles().add(fts);
Layer layer = new FeatureLayer(featureCollection, style);
map.addLayer(layer);
} catch (Exception e) {
e.printStackTrace();
}
}/**
* 根据四至坐标、长、宽像素获取地图内容，并生成图片
*
* @param params
* @param imgPath
*/
public void getMapContent(Map params, String imgPath) {
try {
double[] bbox = (double[]) params.get("bbox");
double x1 = bbox[0], y1 = bbox[1],
x2 = bbox[2], y2 = bbox[3];
int width = Integer.parseInt(params.get("width").toString()) ,
height = Integer.parseInt(params.get("height").toString());
// 设置输出范围
CoordinateReferenceSystem crs = DefaultGeographicCRS.WGS84;
ReferencedEnvelope mapArea = new ReferencedEnvelope(x1, x2, y1, y2, crs);
// 初始化渲染器
StreamingRenderer sr = new StreamingRenderer();
sr.setMapContent(map);
// 初始化输出图像
BufferedImage bi = new BufferedImage(width, height,
BufferedImage.TYPE_INT_ARGB);
Graphics g = bi.getGraphics();
((Graphics2D) g).setRenderingHint(RenderingHints.KEY_ANTIALIASING,
RenderingHints.VALUE_ANTIALIAS_ON);
((Graphics2D) g).setRenderingHint(RenderingHints.KEY_TEXT_ANTIALIASING,
RenderingHints.VALUE_TEXT_ANTIALIAS_ON);
Rectangle rect = new Rectangle(0, 0, width, height);
// 绘制地图
sr.paint((Graphics2D) g, rect, mapArea);
//将BufferedImage变量写入文件中。
ImageIO.write(bi, "png", new File(imgPath));
} catch (Exception e) {
e.printStackTrace();
}
}结果展示
结果中我添加了标题和图例以及行政区划，行政区划是 geoserver 发布的 wms 服务，标题和图例使用 java 的 imageio 包加上去就行了，比较简单就不写了。
补充
有比较多的朋友问了重复的问题，在这里统一解答一下。
- 当结果分级较少时可能出现渲染颜色不正确的情况，原因是仅有两条等值线无法判断数值变化趋势，建议增加临界分级来提升渲染效果。
levPros.put(0.0,"#ffffff");
levPros.put(50.9999,"#ffffff");
levPros.put(51,"#ff0000");- 关于图例标题，主要使用的是Graphics2D、BufferedImage、ImageIO 这三个类完成的，网上相关文章很多。这里就写几个主要的接口：
// 核心java包
import javax.imageio.ImageIO;
import java.awt.*;
import java.awt.image.BufferedImage;
// 初始化画布
BufferedImage bi = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB);
Graphics g = bi.getGraphics();
g.fillRect // 填充矩形
g.drawRect // 画矩形框
g.setColor // 设置画笔颜色
g.drawString // 写字
ImageIO.read // 读图片文件
ImageIO.write // 写图片文件
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
