Java项目实战记录：雷达数据插值
业务背景
之前已经实现了雷达数据的解析和雷达数据后端渲染功能，现在又有一个新的需求。之前是将雷达数据点使用GeoTools渲染成PNG的图片，但这个数据返给前端后不能无极缩放，因为它是个栅格图片；并且如果前端直接渲染雷达点数据也不行，渲染后的数据放大会变成稀疏的点。因此要无极缩放并且不能稀疏，这样的话就需要使用插值方法了，使用插值方法将雷达点数据插值成等值面，这样前端只需要渲染等值面即可。
插值方法
Kriging插值法
克里金插值法是一种基于统计学的插值方法，又称空间局部插值法，原理是利用区域化变量为基础，以变异函数为基本工具，对未知样点进行线性无偏、最优化估计。主要包括计算样本变异函数、根据变异函数对待估计数据建模、利用所建模型进行克里金插值估计和估计方差四大部分。
理论变异函数的基本公式为： ( ,ℎ)=1/2 {[ ( )− ( +ℎ)]^2 }
( )和 ( +ℎ)分别是在点x、x+h处的观测值，h为步长；{[ ( )− ( +ℎ)]^2}为方差。
理论变异函数模型有：线性模型、指数模型、球面模型等。讲过大量实验，系统拟选取球状模型作为理论变异函数：
实验变异函数r(h)反映了区域化变量的空间自相关性：
其中，h为监测点之间的空间间隔距离；N(h)为距离等于h的点对数；Z( )为处于点 处变量的观测值；Z( +h)为与点偏离h处变量的实测值。
假设 0为未观测点， ( =1,2,…, )为其周围的观测点，Z表示计算的区域化变量，则对 0处某个区域化变量的估计值可以写为：
IDW插值法
反距离权重 (IDW) 插值：彼此距离较近的事物要比彼此距离较远的事物更相似。当为任何未测量的位置预测值时，反距离权重法会采用预测位置周围的测量值。与距离预测位置较远的测量值相比，距离预测位置最近的测量值对预测值的影响更大。反距离权重法假定每个测量点都有一种局部影响，而这种影响会随着距离的增大而减小。由于这种方法为距离预测位置最近的点分配的权重较大，而权重却作为距离的函数而减小，因此称之为反距离权重法。
代码实现
wContour
这里的插值方法，我使用的是wContour
wContour库专为地球科学数据分析而设计，支持等高线的生成和操作。等高线，也称作等值线，是一种表示三维空间中相同数值的连续线，在地图上显示为连接所有具有相同数据值的点，常用于气象、地质和其他地球科学领域来展示温度、降水、高度等数据的空间分布。以下是wContour库的一些关键特点和功能：
-
语言支持：wContour库提供了Java和C#版本，方便在不同的编程环境中使用，使得开发者可以在Java或.NET平台上集成等高线功能。
-
等高线生成：库提供了生成等高线的算法，这些算法可以从栅格数据或散点数据中计算出等高线。这意味着你可以从具有不同空间分辨率的数据集中生成等高线，无论数据是规则还是不规则分布的。
-
等高线填充：除了生成等高线外，wContour还支持等高线填充功能，允许开发者在等高线之间使用不同的颜色或图案进行填充，以便更直观地表示不同的数据区间。这在制作地图或地球科学数据可视化时非常有用。
-
算法优化：wContour的等高线算法经过优化，能够处理大规模数据集，保证了处理速度和效率，即使是复杂或大尺度的数据集也能快速生成等高线。
-
数据裁剪和转换：该库还提供了数据裁剪和转换工具，允许开发者裁剪出感兴趣的区域或将数据从一种格式转换为另一种格式，以适应不同的分析需求。
-
多种应用场景：wContour的应用场景广泛，包括气象学、海洋学、地质学、环境科学等多个领域，可以用于制作气候地图、地形图、环境监测图等。
wContour打包成maven
- 从官网下载源码
- 下载相关依赖，使用mvn命令打包成jar包，如果使用idea可以如下图所示来打包
- 使用mvn命令打包到本地maven仓库中，方便项目pom引入
mvn install:install-file -Dfile=D:\code\jar\wContour\wContour-1.7.1.jar -DgroupId=com.github -DartifactId=wContour -Dversion=1.7.1 -Dpackaging=jar
在maven项目中使用
<dependency>
<groupId>com.github</groupId>
<artifactId>wContour</artifactId>
<version>1.7.1</version>
</dependency>
核心代码
将计算所需要的点和插值范围传入，通过wContour使用IDW进行计算，具体步骤如下所示：
/**
* 生成等值面
*
* @param radarDataPojo 训练数据
* @param dataInterval 数据间隔
* @param size 大小，矩形分割的大小，越大越精准
* @param wktBoundary 插值范围 wkt
* @return 等值面JSON字符串
*/
public ArrayList<BisProductData> equivalentSurface(RadarDataPojo radarDataPojo,
double[] dataInterval,
int[] size,
String wktBoundary) {
long start = System.currentTimeMillis();
ArrayList<BisProductData> res = new ArrayList<>();
ArrayList<LonLatData> list = radarDataPojo.getLonLatDataList();
double[][] trainData = new double[list.size()][3];
for (int i = 0; i < list.size(); i++) {
trainData[i][0] = list.get(i).getLon();
trainData[i][1] = list.get(i).getLat();
trainData[i][2] = list.get(i).getVal();
}
try {
double _undefData = -9999.0;
SimpleFeatureCollection polygonCollection = null;
List<PolyLine> cPolylineList = new ArrayList<PolyLine>();
List<Polygon> cPolygonList = new ArrayList<>();
int width = size[0],
height = size[1];
double[] _X = new double[width];
double[] _Y = new double[height];
// 解析WKT字符串
WKTReader reader = new WKTReader();
Geometry boundaryGeom = reader.read(wktBoundary);
// 2. 获取边界的最小和最大坐标
double minX = boundaryGeom.getEnvelopeInternal().getMinX();
double minY = boundaryGeom.getEnvelopeInternal().getMinY();
double maxX = boundaryGeom.getEnvelopeInternal().getMaxX();
double maxY = boundaryGeom.getEnvelopeInternal().getMaxY();
// 3. 创建网格
Interpolate.createGridXY_Num(minX, minY, maxX, maxY, _X, _Y);
double[][] _gridData = new double[width][height];
int nc = dataInterval.length;
// 4. 进行插值计算（IDW 插值）
_gridData = Interpolate.interpolation_IDW_Neighbor(trainData,
_X, _Y, 12, _undefData);
int[][] S1 = new int[_gridData.length][_gridData[0].length];
// 5. 提取轮廓线
List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y,
S1, _undefData);
// 6. 生成等值线
cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc,
dataInterval, _undefData, _borders, S1);
// 7. 对等值线进行平滑
cPolylineList = Contour.smoothLines(cPolylineList);// 平滑
// 8. 生成等值面
cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList,
_borders, dataInterval);
// 9. 将生成的等值面转换为 GeoJSON 格式的字符串
res = getPolygonGeoJsonList(cPolygonList);
logger.info("差值成功, 共耗时" + (System.currentTimeMillis() - start) + "ms");
} catch (Exception e) {
logger.error("Error during calculation of equivalentSurface", e);
}
return res;
}
完整代码
VectorInterpolation插值工具类
import com.fasterxml.jackson.databind.ObjectMapper;
import com.zykj.radar_server.datahandle.pojo.LonLatData;
import com.zykj.radar_server.entity.BisProductData;
import com.zykj.radar_server.entity.pojo.RadarDataPojo;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.WKTReader;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
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.FileWriter;
import java.io.IOException;
import java.util.*;
/**
* 矢量插值方法
* <p>
* 使用IDW进行插值
*/
public class VectorInterpolation {
private static final Logger logger = LoggerFactory.getLogger(VectorInterpolation.class);
/**
* 生成等值面
*
* @param radarDataPojo 训练数据
* @param dataInterval 数据间隔
* @param size 大小，矩形分割的大小，越大越精准
* @param wktBoundary 插值范围 wkt
* @return 等值面JSON字符串
*/
public ArrayList<BisProductData> equivalentSurface(RadarDataPojo radarDataPojo,
double[] dataInterval,
int[] size,
String wktBoundary) {
long start = System.currentTimeMillis();
ArrayList<BisProductData> res = new ArrayList<>();
ArrayList<LonLatData> list = radarDataPojo.getLonLatDataList();
double[][] trainData = new double[list.size()][3];
for (int i = 0; i < list.size(); i++) {
trainData[i][0] = list.get(i).getLon();
trainData[i][1] = list.get(i).getLat();
trainData[i][2] = list.get(i).getVal();
}
try {
double _undefData = -9999.0;
SimpleFeatureCollection polygonCollection = null;
List<PolyLine> cPolylineList = new ArrayList<PolyLine>();
List<Polygon> cPolygonList = new ArrayList<>();
int width = size[0],
height = size[1];
double[] _X = new double[width];
double[] _Y = new double[height];
// 解析WKT字符串
WKTReader reader = new WKTReader();
Geometry boundaryGeom = reader.read(wktBoundary);
// 2. 获取边界的最小和最大坐标
double minX = boundaryGeom.getEnvelopeInternal().getMinX();
double minY = boundaryGeom.getEnvelopeInternal().getMinY();
double maxX = boundaryGeom.getEnvelopeInternal().getMaxX();
double maxY = boundaryGeom.getEnvelopeInternal().getMaxY();
// 3. 创建网格
Interpolate.createGridXY_Num(minX, minY, maxX, maxY, _X, _Y);
double[][] _gridData = new double[width][height];
int nc = dataInterval.length;
// 4. 进行插值计算（IDW 插值）
_gridData = Interpolate.interpolation_IDW_Neighbor(trainData,
_X, _Y, 12, _undefData);
int[][] S1 = new int[_gridData.length][_gridData[0].length];
// 5. 提取轮廓线
List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y,
S1, _undefData);
// 6. 生成等值线
cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc,
dataInterval, _undefData, _borders, S1);
// 7. 对等值线进行平滑
cPolylineList = Contour.smoothLines(cPolylineList);// 平滑
// 8. 生成等值面
cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList,
_borders, dataInterval);
// 9. 将生成的等值面转换为 GeoJSON 格式的字符串
res = getPolygonGeoJsonList(cPolygonList);
logger.info("差值成功, 共耗时" + (System.currentTimeMillis() - start) + "ms");
} catch (Exception e) {
logger.error("Error during calculation of equivalentSurface", e);
}
return res;
}
private ArrayList<BisProductData> getPolygonGeoJsonList(List<Polygon> cPolygonList) {
ArrayList<BisProductData> res = new ArrayList<>();
// List<String> geoJsonList = new ArrayList<>();
// 如果多边形列表为空，返回空列表
if (cPolygonList == null || cPolygonList.isEmpty()) {
return res;
}
try {
for (Polygon pPolygon : cPolygonList) {
BisProductData bisProductData = new BisProductData();
List<Object> ptsTotal = new ArrayList<>();
List<Object> pts = new ArrayList<>();
PolyLine pline = pPolygon.OutLine;
for (PointD ptD : pline.PointList) {
List<Double> 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<Object> cpts = new ArrayList<>();
for (PointD ccptD : cptLine.PointList) {
List<Double> pt = new ArrayList<>();
pt.add(ccptD.X);
pt.add(ccptD.Y);
cpts.add(pt);
}
if (!cpts.isEmpty()) {
ptsTotal.add(cpts);
}
}
}
Map<String, Object> geometry = new HashMap<>();
geometry.put("type", "Polygon");
geometry.put("coordinates", ptsTotal);
String geometryJson = new ObjectMapper().writeValueAsString(geometry);
double hv = pPolygon.HighValue;
double lv = pPolygon.LowValue;
bisProductData.setGeom(geometryJson);
bisProductData.setHvalue(hv);
bisProductData.setLvalue(lv);
res.add(bisProductData);
}
} catch (Exception e) {
e.printStackTrace();
}
System.out.println(res.get(0));
return res;
}
/**
* 将指定的 JSON 字符串追加到文件末尾。
*
* @param strFile 要追加内容的文件路径
* @param strJson 要追加的 JSON 字符串
*/
public static void append2File(String strFile, String strJson) {
// 创建文件对象
File f = new File(strFile);
// 文件写入流
FileWriter fw;
try {
// 初始化文件写入流
fw = new FileWriter(f);
// 将 JSON 字符串写入文件
fw.write(strJson);
// 关闭文件写入流
fw.close();
} catch (IOException e) {
e.printStackTrace();
}
}
}
测试方法
这里使用的是之前的雷达数据
@Test
public void testGeom() {
long start = System.currentTimeMillis();
VectorInterpolation equiSurface = new VectorInterpolation();
ArrayList<LonLatData> list = RadarDataParser.handlePolarProduct("C:\\data\\雷达测试数据\\数据\\20190305_110801_ZHBJ_Z_VOL_2.50.dat");
double[][] trainData = new double[list.size()][3];
for (int i = 0; i < list.size(); i++) {
trainData[i][0] = list.get(i).getLon();
trainData[i][1] = list.get(i).getLat();
trainData[i][2] = list.get(i).getVal();
}
double[] dataInterval = new double[]{-10, -5, 0, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70};
String boundryFile = "D:\\temp\\layers\\POLYGON.shp";
int[] size = new int[]{200, 200};
boolean isclip = false;
try {
String strJson = equiSurface.equivalentSurface(trainData, dataInterval, size, boundryFile, isclip);
// System.out.println(strJson);
String strFile = "D:\\temp\\layers\\20200618_073326_ZHBJ_Z_VOL_250_100.json";
VectorInterpolation.append2File(strFile, strJson);
System.out.println(strFile + "差值成功, 共耗时" + (System.currentTimeMillis() - start) + "ms");
} catch (Exception e) {
e.printStackTrace();
}
}
效果
使用IDW插值生成的效果如下所示：
根据等值面的值进行渲染，效果如下所示：
