用有限差分估计(Finite Difference Estimate)解决地理坐标与平面像素坐标转换过程的误差造成风场粒子向量失真问题

下载NCEP的气象场grib2数据，然后用双线性插值可以绘制相关的气象场图谱https://www.cnblogs.com/davidxu/p/13764587.html，但是风场是二维的向量，包含u和v两个分量。这个用经纬度投影到像素坐标会产生误差，直接绘制效果不太对:(

 

通过插值计算得到风场粒子的预测数据wind = interpolate(longitude, latitude)，然而，由于地理坐标系在投影的过程中发生了失真(distortion)，因此直接用wind绘制风场图是错误的！

我们需要对失真程度进行估算以减小误差，得到近似正确的风场粒子向量。

这里大致记录下用微积分学上的有限差分来估算投影失真程度的代码实现。

/**
 * 解决平面像素坐标（x,y）投影地理坐标系失真造成的风场向量的畸变
 * crsUtils是自定义投影坐标系工具类，可实现地理坐标和平面像素坐标互转
 */
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
    };
}

 改善后的渲染效果：