雷达数据可视化原理详解

Spatial Data 2021-11-04

1725

雷达广泛应用于气象、环境监测、军事方向，雷达数据通常非常的大，如何将雷达数据在Web端进行快速可视化一直是困扰很多技术开发人员的难题，本文选择一环境观测雷达测试原始数据，从0开始阐述如何数据处理和可视化表达。

一数据说明

雷达以一点O为圆心，以固定距离R为半径，以正北方向为0°，顺时针每隔若干角度的间隔发射一条射线，从圆心发出的射线，每隔固定距离d获取一个监测值，图示如下：

雷达数据示意图

以实际环境监测雷达数据的PM2.5数据为例，以雷达点为圆心，每隔2°发射一条射线，射线长度为4995米，每条射线上每隔7.5米一个测量值，可以轻易得到以下数据：

每条射线上监测值数量=4995/7.5=666。
射线总数=360/2=180。
雷达一次扫描获取监测数据总数=180*666=119880。
已知带坐标的数据只有雷达所在的圆形点，其他都是描述性数据。

雷达数据具有空间数据少，描述数据多，放射性数据距离圆心越远，数据密度越稀疏，总体数据量还特别大等特点，如果要在web上进行可视化，很多技术人员可能会一头雾水，因为我们常用于Web可视化的空间数据都是点线面等矢量数据，而雷达数据显然与常见矢量非常不同。

为了能在Web端的地图上对雷达数据可视化，起码需要技术人员解决如下几个问题：

如何将描述类数据转换成具体的空间数据表达？
离散性的射线监测数据如何转变成简单的规律数据？
数据量很大如何实现网络快速传输？
数据量大如何实现web端快速可视化并叠加到地图上？

本文将分章节阐述每个问题如何解决。

二数据空间表达

目前唯一已知一点是雷达站的经纬度点，且已知射线规则，可以暂时用一些代码将雷达扫描圆形区域，射线等描述数据转成实际属性的常用数据，本文作者使用PostGIS快速转换下数据，并可以通过QGIS或者ogr2ogr命令导出为GeoJSON格式备用。

以雷达点为圆心，以4995米为半径，获取雷达圆形区域数据：

select ST_Buffer(
ST_GeomFromText('Point(104.255386 30.8904285)', 4326)::geography,
 4995)::geometry as geom into radar_circle;

由于PostGIS比较强GIS，支持球面等坐标，所以本文严谨的话以球面来缓冲，结果如下：

雷达扫描区域

下一步制作雷达射线，雷达射线是以正北方向为0°，顺时针间隔2°夹角发射射线：

create table radar_line(
  gid serial primary key,
  angle int,
  geom geometry(LineString, 4326)
);
do language plpgsql $$
DECLARE 
  lon numeric;
  lat numeric;
  o geometry;
  _line geometry;
  _buffer geometry;
  _angle numeric;
BEGIN --缓冲一个圆
  o := ST_GeomFromText('Point(104.255386 30.8904285)', 4326);
  for i in 1..180 loop 
    _angle :=(i -1) * 2;
    --暂时以5公里为半径，之后裁剪下
    lon := 104.255386 + 5  111.0 * cos((90-_angle)  180.0 * pi());
    lat := 30.8904285 + 5  111.0 * sin((90-_angle)  180.0 * pi());
    _line := ST_MakeLine(
      o,ST_SetSrid(ST_MakePoint(lon, lat), 4326)
    );
    insert into radar_line(angle, geom) values (_angle, _line);
  end loop;
end;
$$;


--使用雷达圆形区域对射线简单裁剪下
update radar_line t1 set geom = 
ST_Intersection(t1.geom,t2.geom) 
from radar_circle t2;

射线制作结果如下：

雷达扫描射线

射线上根据射线方向，射线起点（圆心），间隔距离是可以算出每个监测点的空间位置的，但计算量实际有点大，所以工程上实际不作如此处理，具体处理可参考下一章节。目前在PostGIS中有了两张表，分别存储了雷达扫描区域和雷达扫描射线，使用QGIS导出为GeoJSON备用。

三制作格点数据

雷达的监测值是在每个射线上的，但是射线数据在图像表达上不够连续，例如由射线图可知，射线越长，射线间间隔越大，数据越稀疏，如果按照射线渲染，射线边缘间隔会有大量的分叉空白，但实际数据是空间连续的，因此，不能用射线去实际渲染，工程上实际是将射线数据转换成规律的格点数据，本章节简述一下转换原理。

3.1 制作格点

首先根据圆形区域，计算其bbox，由射线上666个监测值，可得到南北，东西方向格点数量皆为666*2=1332，则可以很容易得到每个格点值的坐标：

// 雷达中心点
const circle_pt = [104.255386, 30.8904285];
// 雷达圆形区域bbox，可以通过postgis或者turf都能快速计算出来
const bbox = [
    104.20988464355469, 30.84493637084961, 104.30088806152344, 30.93592071533203
];
// 格网分辨率
const res_x = (bbox[2] - bbox[0])  1332;
const res_y = (bbox[3] - bbox[1]) / 1332;


 for (let i = 0; i < 1332; i++) {
    for (let j = 0; j < 1332; j++) {
         格点坐标
        const gridGeom = [
            bbox[0] + res_x * j,
            bbox[1] + res_y * i
        ];
        // 格点序号
        const gridIndex = j + i * 1332;
    }
}

格点生成示意图

3.2 格点赋值

格点坐标获取很容易，格点赋值就稍微麻烦点，这里首先简述下格点赋值原理：

格点赋值说明图

首先，对每一个格点，计算其是否落在圆形区域，不在圆形区域则不计算，减少计算量。

其次，落在圆形区域，则计算圆形到格点的相对正北方向角，每条射线都有其正北方向角，简单计算对比可得到格点与哪条射线最接近，则该格点值从该射线上获取。

最后，获取格点与射线的垂线交点，计算垂足与圆心位置距离l1，射线总长度l2，则可得到垂足在射线上的位置百分比per=l1/l2，由于射线上有666个监测值，共同组成一个数组datas，而监测点分辨率res=1/666，那么可以得到监测值序号再去datas中捞取格点值：

const index = Math.round(per/ res);
// 射线监测值数据为数组datas
const gridvalue = datas[index];

逻辑完整代码如下基于NodeJS+Turf做的计算：

// qgis导出的雷达扫射圆形区域json数据
const circle = require('./circle.json');
// qgis导出的雷达射线json数据
const lines = require('./line.json');
const point_in_polygon = require('@turf/boolean-point-in-polygon');
const turfBearing = require('@turf/bearing');
const nearest_point_on_line = require('@turf/nearest-point-on-line');
const turfDistance = require('@turf/distance');
const xlsx = require('node-xlsx');
const fs = require('fs');


// 雷达扫射区域bbox
const circle_pt = [104.255386, 30.8904285];
const bbox = [
    104.20988464355469, 30.84493637084961, 104.30088806152344, 30.93592071533203
];


let lineAngles = new Map();
lines.features.forEach(feature => {
    const length = turfDistance.default(feature.geometry.coordinates[0], feature.geometry.coordinates[1]);
    feature.properties.length = length;
    lineAngles.set(feature.properties.angle, feature);
});


// 解析雷达射线数据
let sheetList = xlsx.parse('./radar_county_analysis_data.xlsx');
let lineDatas = new Map();
const datas = sheetList[0].data;
for (let i = 1; i < datas.length; i++) {
    const angle = parseInt(datas[i][15]);
    let values = datas[i][17];
    values = values.split(',');
    lineDatas.set(angle, values);
}




async function main() {
    const res_x = (bbox[2] - bbox[0]) / 1332;
    const res_y = (bbox[3] - bbox[1]) / 1332;
    // 测点分辨率
    const res = 1.0 / 666;
    // 格点值数组
    let datas = new Array(1332 * 1332);
    for (let i = 0; i < 1332; i++) {
        for (let j = 0; j < 1332; j++) {
            const gridGeom = [
                bbox[0] + res_x * j,
                bbox[1] + res_y * i
            ];
            const gridIndex = j + i * 1332;
            // 只有点在圆形区域内才计算
            if (point_in_polygon.default(gridGeom, circle)) {
                // 计算夹角
                let angle = turfBearing.default(circle_pt, gridGeom);
                //console.log(circle_pt,gridGeom,angle);
                let mod = angle % 2;
                if (mod >= 1)
                    angle = angle - mod + 2;
                else
                    angle = angle - mod;
                if (angle < 0)
                    angle = 360 + angle;


                // 计算线上最近的点
                const nearest_feature = lineAngles.get(angle);
                const nearest_pt = nearest_point_on_line.default(nearest_feature, gridGeom);
                // 最近点距离起点距离，归并到最近点位上
                const _length = turfDistance.default(circle_pt, nearest_pt);
                const per = _length / nearest_feature.properties.length;
                const index = Math.round(per / res);
                let value = Number(lineDatas.get(angle)[index]);
                if (!value || value > 20000)
                    continue;
                datas[gridIndex] = value;
            }
        }
    }