前端-风场可视化研究1
实时风场可视化（如 earth.nullschool.net），非常漂亮，生动地展示了风向的变化
最近想研究一下风场可视化，风力图；然后看了以下三个工程
https://github.com/Esri/wind-js
https://github.com/mapbox/webgl-wind
https://github.com/RaymanNg/3D-Wind-Field
wind-js的数据处理与渲染是基于CPU 的，后两者是基于gpu；
下图是wind-js的运行效果图
wind-js 中的 evolve
function evolve() {
buckets.forEach(function(bucket) { bucket.length = 0; });
particles.forEach(function(particle) {
if (particle.age > MAX_PARTICLE_AGE) {
field.randomize(particle).age = 0;
}
var x = particle.x;
var y = particle.y;
var v = field(x, y); // vector at current position
var m = v[2];
if (m === null) {
particle.age = MAX_PARTICLE_AGE; // particle has escaped the grid, never to return...
}
else {
var xt = x + v[0];
var yt = y + v[1];
if (field(xt, yt)[2] !== null) {
// Path from (x,y) to (xt,yt) is visible, so add this particle to the appropriate draw bucket.
particle.xt = xt;
particle.yt = yt;
buckets[colorStyles.indexFor(m)].push(particle);
}
else {
// Particle isn't visible, but it still moves through the field.
particle.x = xt;
particle.y = yt;
}
}
particle.age += 1;
});
}
负责算出每个风粒子的下一帧的坐标，如果生命age 达到阈值，则用一个随机的位置赋予给该粒子；避免运行一段时间后，所有粒子挤在一起或所有粒子跑出屏幕。
该函数中的 field 对象有两次引用，
field.randomize(particle).age = 0;
用一个随机的位置赋予给 particle, 并将particle的age属性置为零
var v = field(x, y)
根据particle的屏幕像素坐标（x，y) 查找出该像素坐标下对应的风速 v，风速v是包含三个元素的数组 [速度的x轴分量，速度的y轴分量，速度的大小]
wind-js工程的gfs.json, 这个数据的原始分辨率为 360 x 181, 依据是
而我的 canvas的分辨率是 2560 x 1313。
field函数的作用是，把任意屏幕像素坐标 x,y 映射到 360 x 181 的网格中，映射后的x坐标记为 x1，映射后y坐标记为 y1; 通常情况下，x1 和 y1 都不是整数，需要 分别对x1 和 y1 使用双线性插件求出对应的 vx（速度的x分量） 和 vy （速度的x分量）
代码中的
var invert = function(x, y, windy)
将屏幕上的像素坐标 (x, y) 转换回地理坐标 (lon, lat)， 也就是负责 (x,y) 映射到 （x1, y1)，（x1, y1)是经纬度坐标，
使用了墨卡托投影公式的逆运算。
var invert = function(x, y, windy){
var mapLonDelta = windy.east - windy.west;
var worldMapRadius = windy.width / rad2deg(mapLonDelta) * 360/(2 * Math.PI);
var mapOffsetY = ( worldMapRadius / 2 * Math.log( (1 + Math.sin(windy.south) ) / (1 - Math.sin(windy.south)) ));
var equatorY = windy.height + mapOffsetY;
var a = (equatorY-y)/worldMapRadius;
var lat = 180/Math.PI * (2 * Math.atan(Math.exp(a)) - Math.PI/2);
var lon = rad2deg(windy.west) + x / windy.width * rad2deg(mapLonDelta);
return [lon, lat];
};
windy的 east, west, south, north 属性来自于 esri/map
extent 是一个经纬度地理范围，当用户缩放地图或平移地图时，extent都会变化
var extent = map.geographicExtent;
setTimeout(function(){
windy.start(
[[0,0],[map.width, map.height]],
map.width,
map.height,
[[extent.xmin, extent.ymin],[extent.xmax, extent.ymax]]
);
},500);
var project = function( lat, lon, windy)
var deg2rad = function( deg ){
return (deg / 180) * Math.PI;
};
// 墨卡托投影：纬度转 y 坐标
var mercY = function( lat ) {
return Math.log( Math.tan( lat / 2 + Math.PI / 4 ) );
};
// 将经纬度投影为屏幕坐标
var project = function( lat, lon, windy) { // both in radians, use deg2rad if neccessary
// 对纬度使用墨卡托投影，得到 y 值
var ymin = mercY(windy.south);
var ymax = mercY(windy.north);
// 每单位经度对应多少像素（x 方向）
var xFactor = windy.width / ( windy.east - windy.west );
// 每单位墨卡托 y 值对应多少像素（y 方向）
var yFactor = windy.height / ( ymax - ymin );
//反转 y 坐标方向（因为墨卡托 y 向上增长，屏幕坐标向下增长）
var y = mercY( deg2rad(lat) );
var x = (deg2rad(lon) - windy.west) * xFactor;
y = (ymax - y) * yFactor; // y points south
return [x, y];
};
project 是把地理坐标（经纬度）投影到 屏幕坐标，由于使用了墨卡托投影，会导致地图上的高纬度地区的像素拉伸，纬度越高，纵向拉伸地越厉害
投影逻辑：
将经纬度转换为弧度
对纬度使用墨卡托投影，得到 y 值
线性缩放 x 和 y 到屏幕坐标系
xFactor：每单位经度对应多少像素（x 方向）
yFactor：每单位墨卡托 y 值对应多少像素（y 方向）
反转 y 坐标方向（因为墨卡托 y 向上增长，屏幕坐标向下增长）
var interpolateField = function( grid, bounds, extent, callback )
var interpolateField = function( grid, bounds, extent, callback ) {
var velocityScale = VELOCITY_SCALE;
var columns = [];
var x = bounds.x;
function interpolateColumn(x) {
var column = [];
// 每次步进2个像素 降采样, 减少插值次数，提高性能
for (var y = bounds.y; y <= bounds.yMax; y += 2) {
var coord = invert( x, y, extent );
if (coord) {
var λ = coord[0], φ = coord[1];
if (isFinite(λ)) {
var wind = grid.interpolate(λ, φ);
if (wind) {
wind = distort(λ, φ, x, y, velocityScale, wind, extent);
column[y+1] = column[y] = wind;
}
}
}
}
columns[x+1] = columns[x] = column;
}
(function batchInterpolate() {
var start = Date.now();
while (x < bounds.width) {
interpolateColumn(x);
x += 2;
if ((Date.now() - start) > 1000) { //MAX_TASK_TIME) {
setTimeout(batchInterpolate, 25);
return;
}
}
createField(columns, bounds, callback);
})();
};
是对整个屏幕区域内所有点插值得到 field风速数据，field其实是一个二维数组，每个元素也是一个数组 [速度的x轴分量，速度的y轴分量，速度的大小]，
是一列 一列 存储的 [画布第一列像素，画布第二列像素，画布第三列像素......]
interpolateField的逻辑是：
使用 invert 得到每个像素点对应的地理坐标
使用 grid.interpolate 获取该点（地理坐标）的风速
使用 distort 修正project方法中使用的墨卡托投影 导致的向量变形
批量异步处理 batchInterpolate，防止页面卡顿
var distortion = function(λ, φ, x, y, windy)
/**
* Calculate distortion of the wind vector caused by the shape of the projection at point (x, y). The wind
* vector is modified in place and returned by this function.
*/
var distort = function(projection, λ, φ, x, y, scale, wind, windy) {
var u = wind[0] * scale;
var v = wind[1] * scale;
var d = distortion(projection, λ, φ, x, y, windy);
// Scale distortion vectors by u and v, then add.
wind[0] = d[0] * u + d[2] * v;
wind[1] = d[1] * u + d[3] * v;
return wind;
};
var distortion = function(λ, φ, x, y, windy) {
var τ = 2 * Math.PI;
var H = Math.pow(10, -5.2);
var hλ = λ < 0 ? H : -H;
var hφ = φ < 0 ? H : -H;
var pλ = project(φ, λ + hλ,windy);
var pφ = project(φ + hφ, λ, windy);
// 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 * τ);
return [
(pλ[0] - x) / hλ / k,
(pλ[1] - y) / hλ / k,
(pφ[0] - x) / hφ,
(pφ[1] - y) / hφ
];
};
计算投影引起的向量变形矩阵
使用微小偏移（hλ, hφ）来近似地图投影的局部变形
通过坐标差值来计算出局部的雅可比矩阵（Jacobian matrix）
这个矩阵描述了该点附近的地图投影是如何“拉伸”或“扭曲”的
distortion的返回值：[dx/dλ, dy/dλ, dx/dφ, dy/dφ]：一个 2x2 矩阵，用于修正风向向量
interpolateColumn 方法里为什么需要调用 distort方法？
因为 var wind = grid.interpolate(λ, φ); 这里双线性插值完得到的 wind 速度单位是 m/s, 而 esri/map 绘制的地图是墨卡托投影的，
上图为 esri/map 绘制的地图，观察这个地图，可以发现高纬度的纬线 和 赤道0°纬线 一样长，显然纬度越高，横向的拉伸越厉害。
双线性插值完得到的 wind代表的速度向量 也需要进行相对应的拉伸来 匹配地图的墨卡托投影，看起来才协调。
如果不调用distort 修改wind，风粒子的流动就和二维地图无法匹配。
打个比方，在北纬60°，恰好有一个风粒子，运动方向是正东方向，速率为1000m/s, cos 60° = 0.5
恰好在赤道上有另一个风粒子，运动方向也是正东方向，速率也为1000m/s
不调用distort时，是两个条平行线，并且随着时间的流逝，两条平行线轨迹始终一样长。
这样就会和地图没有匹配起来。因为当赤道上的风粒子绕地球一圈时，北纬60°的风粒子 应该绕地球两圈 才对，
不调用distort时，赤道上的风粒子 和 北纬60°的风粒子 同时绕地球一圈，这显然不对。应该对北纬60°的风粒子的速率进行两倍的拉伸
var wind = grid.interpolate(λ, φ); 得到的wind速度位于 均匀的 线性的空间中；需要调用 distort 把 wind 转换为墨卡托投影的 非均匀的 非线性的空间中，这样才和二维的墨卡托投影地图 匹配的上
充分理解wind-js示例，对理解后续两个webgl的 风力图示例很有帮助。
浙公网安备 33010602011771号