气象风场技术（年代久远，归档）
气象风场技术
一、NCEP提供的GRIB2格式的GFS数据介绍
1、名词介绍：
NCEP：National Weather Service
GRIB2: 是一种二进制文件格式
GFS：Global Forecast System
2、风场数据：
GRIB2格式经过转换成易读的JSON格式，风场JSON是包含2个JSON对象的JSON数组，一个表示风场的U组件，一个表示V组件， U/V可理解为二维坐标系的x/y轴。
[ { "header": { "discipline": 0, "disciplineName": "Meteorological products", "gribEdition": 2, "gribLength": 5800, "center": 7, "centerName": "US National Weather Service - NCEP(WMC)", ... "parameterNumber": 2, "parameterNumberName": "U-component_of_wind", "parameterUnit": "m.s-1", ... "gridDefinitionTemplateName": "Latitude_Longitude", "numberPoints": 4088, //格网中共计4088个格点 "shape": 6, "shapeName": "Earth spherical with radius of 6,371,229.0 m", "gridUnits": "degrees", "resolution": 48, "winds": "true", "scanMode": 64, "nx": 73, //格网横轴73个点（73 = 139 - 67 + 1） "ny": 56, //格网纵轴56个点（56 = 55 - 0 + 1） "basicAngle": 0, "lo1": 67, //经度起始位置 "la1": 0, //纬度起始位置 "lo2": 139, //经度结束位置 "la2": 55, //纬度结束位置 "dx": 1, // 间隔经纬度为1度 "dy": 1 //间隔经纬度为1度 }, "data": [-6.284585, -5.914585, -5.754585, -5.914585, -5.634585,...] //data表示风场在相应经纬度位置处的风速(m/s) }, { "header": {...}, "data": [...] }
二、根据风场数据构造二维网格，网格中的每个格点对应经纬度和相应的风场参数(风速m/s的u和v分量）
function buildGrid(json){ var header = json.header; var λ0 = header.lo1, φ0 = header.la1 > header.la2 ? header.la1 : header.la2; // 网格原点,取投影平面左上角 (e.g., 0.0E, 90.0N) var Δλ = header.dx, Δφ = header.dy; // 格点之间的距离 (e.g., 1 deg lon, 1 deg lat) var ni = header.nx, nj = header.ny; // 东西向/南北向格点的数量 (e.g., 73 x 56) // 经度从λ0开始递增, 纬度从φ0开始递减. var grid = [], p = 0; var isContinuous = Math.floor(ni * Δλ) >= 360;//是否连续，若经度范围超过360了，就是饶了地球一圈，首尾肯定会连起来，即表示连续 if(header.la1 > header.la2) {//起始纬度大于结束纬度，即方向是从北半球指向南半球的，GFS数据相应是按二维数组的行/列规则进行存储的 for (let j = 0; j < nj; j++) { var row = []; for (let i = 0; i < ni; i++, p++) { row[i] = json.data(p); } if (isContinuous) { row.push(row[0]); //对网格数据是收尾连续的情形，复制第一列作为最后一列，以简化插值逻辑 } grid[j] = row; } } else {//起始纬度小于结束纬度，即方向是从南半球指向北半球的，GFS数据存储的顺序不完全按二维数组的行列规则来的，行的索引是反的，由下向上递增 for (let k = nj -1; k >= 0; k--) { var row = []; for (let i = 0; i < ni; i++, p++) { row[i] = json.data(p); } if (isContinuous) { row.push(row[0]); //对网格数据是收尾连续的情形，复制第一列作为最后一列，以简化插值逻辑 } grid[k] = row; } } }
三、插值计算
1、我们要将给定经纬度范围的风场数据绘制到平面像素空间，首先是把给定的经纬度范围投影到平面坐标系，通常对WGS84地理坐标系运用墨卡托投影得到像素坐标，从而获得像素坐标范围bounds；
2、bounds有x轴的起始值x和最大值xMax，y轴的起始值y和最大值yMax， 单位是像素，从而可得到空间分辨率为(xMax - x) * (yMax - y)；
3、由于风场数据的地理间隔是1度一个数值，经过投影成像素坐标后，一度之间还是有很多点没数据，若忽略不计，则风场最终的渲染效果可能不好，所以还得通过插值补全它；
4、考虑到插值计算量过大会产生严重的性能问题，对于给定分辨率的像素空间，可以设置合适的步长，比如x和y轴每隔2个像素插值一个，因此就有下面这样的代码：
for (let x = bounds.x; x < bounds.xMax;) { interpolateColumn(x); x += 2; //x轴方向递增步长为2 } function interpolateColumn() { for (let y = bounds.y; y <= bounds.yMax; y += 2) { //y轴方向递增步长为2 //TODO interpoate } }
5、对风场数据应用了双线性插值算法，风场数据包含U和V两个分量，正好对应双线性算法插值计算需要2个变量。
1）线性插值（Linear Interpolate)
在数学中，线性插值是一种曲线拟合方法，利用线性多项式在已知数据点的离散集合范围内构造新的数据点。
两个已知点之间的线性插值：
已知两点由坐标(x0，y0)和(x1，y1)给出，线性插值就是两点之间的直线。对于区间(x0，x1)中的x值，由方程给出沿直线的y值
2）双线性插值（Bilinear Interpolate)
在数学上，双线性插值是线性插值的拓展，是针对2维直线网格上的2个自变量函数z = f(x, y)的插值。
主要思想是先在一个方向上执行线性插值，再在另一个方向上操作一次。虽然样本值和样本位置中的每一步插值是线性的，但是整体上是非线性的。
注：红色点是原有的数据点，绿色点是我们想要插值的结果。
算法思路：
令我们要找的未知值的函数关系 f 和相应坐标（x, y)， 假设已知的4个点是Q11 = (x1, y1)，Q12=(x2, y2)，Q21=(x2, y1)， Q22=(x2, y2)。
首先在x轴上做线性插值，得到如下式子：
我们继续在y轴上插值得到估算式子：
注：先在x轴还是y轴上插值的最终结果是一致的，顺序不影响插值结果。
另一种解决插值问题的方法是解下面的方程：
系数a0,a1,a2,a3可以通过解线性系统得到：
若一个解用f(Q)表示，我们可以写成：
系数可以通过解下面的线性系统获得：
上式推导过程如下：
若我们选择一个坐标系统来表示这4个已知的点，则插值公式可以简化为：
注：风场那块插值就用了这个公式
或者等价的矩阵操作：
3）对风场的二维网格应用双线性插值，生成预测数据
function interpolate(λ, φ) { //λ表示longitude参数, φ表示latitude参数 var i = floorMod(λ - λ0, 360) / Δλ; //计算longitude参数在[0, 360)范围内的索引 var j = (φ0 - φ) / Δφ; //计算纬度在+90到-90范围中的索引 // 1 2 在把λ和φ转换成网格分数索引，我们找到围绕着点(i,j)的4个拐点G，这4个拐点是通过对i和j进行 // fi i ci floor和ceil取整操作获得的。例如， i = 1.4， j = 8.3时，四个拐点为(1, 8), (2, 8), // | =1.4 | (1, 9) , (2, 9). // ---G---|---G--- fj 8 // j ___|_ _. | 注意：对于包裹着的网格，第一列被复制为最后一列，故可以使用索引ci而无须取模 // =8.3 | | // ---G-------G--- cj 9 // | | var fi = Math.floor(i), ci = fi + 1; var fj = Math.floor(j), cj = fj + 1; var row; if ((row = grid[fj])) {//网格中能取到行索引为fj的行 var g00 = row[fi]; var g10 = row[ci]; if (isValue(g00) && isValue(g10) && (row = grid[cj])) { var g01 = row[fi]; var g11 = row[ci]; if (isValue(g01) && isValue(g11)) { // 网格中的四个点找到了，开始用这4个拐点插值. return bilinearInterpolateVector(i - fi, // 以(fi, fj)为原点时，此时(fi, fj) = (0, 0), i的位置索引 j - fj, // 以(fi, fj)为原点时，此时(fi, fj) = (0, 0), j的位置索引 g00, g10, g01, g11); //四个拐点 } } } return null; } // 参考前面介绍的双线性插值算法简化公式：f(x,y) = f(0,0) * (1-x) * (1-y) + f(1,0) * x * (1-y) + f(0,1) * (1-x) * y + f(1,1) * x * y // g00、g10、g01、g11分别是f(0,0)、f(1,0)、f(0,1)、f(1,1)，因此代入多项式很容易算得f(x,y)的值，即(x,y)索引位置处的数值。 function bilinearInterpolateScalar(x, y, g00, g10, g01, g11) { // 双线性插值，针对标量，一个坐标 var rx = (1 - x); var ry = (1 - y); return g00 * rx * ry + g10 * x * ry + g01 * rx * y + g11 * x * y; // } // 对于矢量我们要分别对2个分量进行求解 // 双线性插值，针对矢量，U/V共2个分量的计算，风场用这个 function bilinearInterpolateVector(x, y, g00, g10, g01, g11) { var rx = (1 - x); var ry = (1 - y); var a = rx * ry, b = x * ry, c = rx * y, d = x * y; var u = g00[0] * a + g10[0] * b + g01[0] * c + g11[0] * d; var v = g00[1] * a + g10[1] * b + g01[1] * c + g11[1] * d; return [u, v, Math.sqrt(u * u + v * v)]; //根据风场插值得到位置(x,y)处的U/V两个分量，并用勾股定理算得U/V的合向量的模。 }
四、风场粒子空间计算
1、有限差分估计(Finite Difference Estimate)
通过插值计算得到风场粒子的预测数据wind = interpolate(longitude, latitude)，然而，由于地理坐标系在投影的过程中发生了失真(distortion)，因此直接用wind绘制风场图是错误的！
我们需要对失真程度进行估算以减小误差，得到近似正确的风场粒子向量。
这里大致介绍下用微积分学上的有限差分来估算投影失真程度的代码实现。
/** * 解决平面像素坐标（x,y）投影地理坐标系失真造成的风场向量的畸变 * crsUtils是自定义投影坐标系工具类，可实现地理坐标和平面像素坐标互转，由于没有加密问题，性能比高德地图API要好 */ function distort(crsUtils, λ, φ, x, y, scale, wind, options) { let u = wind[0] * scale; //放大U分量 let v = wind[1] * scale; //放大v分量 let d = __distortion(crsUtils, λ, φ, x, y, options); //根据估算的结果重新构造U/V wind[0] = d[0] * u + d[2] * v; wind[1] = d[1] * u + d[3] * v; return wind; } /** * 返回在给定点应用某个特定投影产生的失真 * * 该方法使用了有限差分估计思想，通过在经度和纬度上增加非常小的量h来创建两条线段以计算扭曲程度。这些线段随后被投影到像素空间， * 在那里它们变成三角形的对角线以表示在那个位置投影扭曲经度和纬度的程度。 * * <pre> * (λ, φ+h) (xλ, yλ) * . . * | ==> \ * | \ __. (xφ, yφ) * (λ, φ) .____. (λ+h, φ) (x, y) .-- * </pre> * * See: * Map Projections: A Working Manual, Snyder, John P: pubs.er.usgs.gov/publication/pp1395 * gis.stackexchange.com/questions/5068/how-to-create-an-accurate-tissot-indicatrix * www.jasondavies.com/maps/tissot * * @returns {Array} array of scaled derivatives [dx/dλ, dy/dλ, dx/dφ, dy/dφ] */ function __distortion(crsUtils, λ, φ, x, y) {//λ表示经度, φ表示纬度 var hλ = λ < 0 ? H : -H; var hφ = φ < 0 ? H : -H; var pλ = __project([λ + hλ, φ], crsUtils); var pφ = __project([λ, φ + hφ], crsUtils); // Meridian scale factor (see Snyder, equation 4-3), where R = 1. This handles issue where length of 1° λ // changes depending on φ. Without this, there is a pinching effect at the poles. var k = Math.cos(φ / 360 * 2*Math.PI); //计算 return [ (pλ.x - x) / hλ / k, (pλ.y - y) / hλ / k, (pφ.x - x) / hφ, (pφ.y - y) / hφ ]; } function __project(lngLatArr, crsUtils) { var p = crsUtils.project(Array.from(lngLatArr).reverse()); return { x: p.x - crsUtils.pixelOrigin.x, y: p.y - crsUtils.pixelOrigin.y }; }
浙公网安备 33010602011771号