1. 程式人生 > >web前端利用leaflet生成粒子風場,類似windy

web前端利用leaflet生成粒子風場,類似windy

wind.js如下:

$(function() {
    var dixing = L.tileLayer.chinaProvider('Google.Satellite.Map', {
        maxZoom: 18,
        minZoom: 2
    });
    var map = L.map("map", {
        center: [33.5, 114.6],
        zoom: 10,
        maxZoom: 20,
        minZoom: 3,
        layers: dixing,
        zoomControl: 
false, attributionControl: false, // crs:L.CRS.EPSG4326 }); map.setView([-19, 150], 5); //各種定義和引數 var grid = []; //定義經緯度網格陣列 var particles = []; //存放粒子陣列 var PARTICLE_MULTIPLIER = 0.4 / 300; //粒子數量引數,預設1/300,可以根據實際調 var max_age = 100; //定義粒子的最大運動次數 //沒有風的情況 var
NULL_WIND_VECTOR = [NaN, NaN, null]; var min_color = 0; // 風速為0使用的顏色 var max_color = 10; // 風速最大使用的顏色 var linewidth = 1.5; //定義粒子粗細 var FRAME_RATE = 35, //定義每秒執行的次數 FRAME_TIME = 1000 / FRAME_RATE; // desired frames per second var galpha = 0.92; //定義透明度,透明度越大,尾巴越長 //存放顏色的陣列 var
colorScale = ["rgb(36,104, 180)", "rgb(60,157, 194)", "rgb(128,205,193 )", "rgb(151,218,168 )", "rgb(198,231,181)", "rgb(238,247,217)", "rgb(255,238,159)", "rgb(252,217,125)", "rgb(255,182,100)", "rgb(252,150,75)", "rgb(250,112,52)", "rgb(245,64,32)", "rgb(237,45,28)", "rgb(220,24,32)", "rgb(180,0,35)"]; //粒子動畫的速度,為引數乘以螢幕的畫素密度開三次方 //如果不支援畫素密度,就乘以1, //預設0.0005 // var vscale = (0.01) * (Math.pow(window.devicePixelRatio, 1 / 3) || 1); // scale for wind velocity (completely arbitrary--this value looks nice) var vscale = 0.01; var animationLoop; //定義動畫 //第一步生成經緯度網格/////////////////////////////////////////////////// $.getJSON("leaflet-velocity-master/demo/wind-global.json", function(data) { setTimeout(function() { buildGrid(data); map.fire("moveend"); }, 100) // bulid_grid(); // create_new_lizi(); // var then = Date.now(); // (function frame() { // animationLoop = requestAnimationFrame(frame); // var now = Date.now(); // var delta = now - then; // if(delta > FRAME_TIME) { // then = now - delta % FRAME_TIME; // evolve(); // draw(); // } // })(); // map.fire("moveend"); }) var lon_min, lat_min, x_kuadu, y_kuadu, ni, nj; function createWindBuilder(uComp, vComp) { var uData = uComp.data, vData = vComp.data; return { header: uComp.header, //recipe: recipeFor("wind-" + uComp.header.surface1Value), data: function data(i) { return [uData[i], vData[i]]; } }; }; function createBuilder(data) { var uComp = null, vComp = null, scalar = null; //從這裡判斷出,資料中的這兩個引數,必須是固定的,否則會不執行 data.forEach(function(record) { switch(record.header.parameterCategory + "," + record.header.parameterNumber) { case "1,2": case "2,2": uComp = record; break; case "1,3": case "2,3": vComp = record; break; default: scalar = record; } }); return createWindBuilder(uComp, vComp); }; function buildGrid(data, callback) { builder = createBuilder(data); var header = builder.header; lon_min = header.lo1; //小的經度 lat_min = header.la1; // the grid's origin (e.g., 0.0E, 90.0N) //小緯度 x_kuadu = header.dx; //x方向的跨度 y_kuadu = header.dy; // distance between grid points (e.g., 2.5 deg lon, 2.5 deg lat) //y方向的跨度 ni = header.nx; //x方向的總數 nj = header.ny; // number of grid points W-E and N-S (e.g., 144 x 73) // //y方向的總數 // date = new Date(header.refTime); // date.setHours(date.getHours() + header.forecastTime); // // // Scan mode 0 assumed. Longitude increases from lon_min, and latitude decreases from lat_min. // // http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-4.shtml var p = 0; var isContinuous = Math.floor(ni * x_kuadu) >= 360; //x方向的跨度乘以x方向的數量是否大於360 for(var j = 0; j < nj; j++) { var row = []; for(var i = 0; i < ni; i++, p++) { row[i] = builder.data(p); } if(isContinuous) { // For wrapped grids, duplicate first column as last column to simplify interpolation logic row.push(row[0]); } grid[j] = row; } //grid是一個三維陣列 //第一緯表示行數 //第二緯表示列數 //第三緯表示每一個網格點的uv }; //////////////// //第二步生成插值網格////////////////////////////////////// var columns = []; //存放插值網格陣列 var start = Date.now(); var bwidth = map.getSize().x, bheight = map.getSize().y; var bx = 0, by = 0; var projection = {}; var mbounds = map.getBounds(); var now_bounds = { east: deg2rad(mbounds._northEast.lng), west: deg2rad(mbounds._southWest.lng), north: deg2rad(mbounds._northEast.lat), south: deg2rad(mbounds._southWest.lat), height: bheight, width: bwidth } var mapArea = (now_bounds.south - now_bounds.north) * (now_bounds.west - now_bounds.east); var velocityScale = vscale * Math.pow(mapArea, 0.4); function bulid_grid() { columns = []; //存放插值網格陣列 start = Date.now(); bx = 0, by = 0; projection = {}; mbounds = map.getBounds(); now_bounds = { east: deg2rad(mbounds._northEast.lng), west: deg2rad(mbounds._southWest.lng), north: deg2rad(mbounds._northEast.lat), south: deg2rad(mbounds._southWest.lat), height: bheight, width: bwidth } mapArea = (now_bounds.south - now_bounds.north) * (now_bounds.west - now_bounds.east); velocityScale = vscale * Math.pow(mapArea, 0.4); batchInterpolate() } function batchInterpolate() { while(bx < bwidth) { interpolateColumn(bx); bx += 2; if(Date.now() - start > 1000) { //MAX_TASK_TIME) { setTimeout(batchInterpolate, 25); //如果粒子生成時間大於1秒了,就不再生成 return; } } } function interpolateColumn(x) { var column = []; //畫布上的每兩個畫素是一個格點 for(var y = 0; y <= bheight; y += 2) { var coord = invert(x, y, now_bounds); //求出每一個格點所對應的經緯度 if(coord) { var λ = coord[0], //經度 φ = coord[1]; //緯度 //監測經度是不是數字 if(isFinite(λ)) { var wind = interpolate(λ, φ); //每一個格點的uv和風速大小 if(wind) { wind = distort(projection, λ, φ, x, y, velocityScale, wind, now_bounds); //wind 表示粒子x方向的畫素速度,y方向上的畫素速度和風速 column[y + 1] = column[y] = wind; } } } } //相鄰兩個格點用相同的速度 columns[x + 1] = columns[x] = column; } function invert(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]; }; function deg2rad(deg) { return deg / 180 * Math.PI; }; function rad2deg(ang) { return ang / (Math.PI / 180.0); }; function floorMod(a, n) { return a - n * Math.floor(a / n); }; var isValue = function isValue(x) { return x !== null && x !== undefined; }; var mercY = function mercY(lat) { return Math.log(Math.tan(lat / 2 + Math.PI / 4)); }; function interpolate(λ, φ) { //如果風場網格資料不存在,就不執行 if(!grid) return null; //λ, φ分別是每個網格點的經度和緯度 var i = floorMod(λ - lon_min, 360) / x_kuadu; // calculate longitude index in wrapped range [0, 360) //lat_min應該是上面的緯度,大的一個 var j = (lat_min - φ) / y_kuadu; // calculate latitude index in direction +90 to -90 //計算網格點在從上到下,從左到右,以最小刻度為0的第幾個經緯度格點上 var fi = Math.floor(i), //格點的上一行 ci = fi + 1; //格點的下一行 var fj = Math.floor(j), //格點的前一列 cj = fj + 1; //格點的後一列 var row; if(row = grid[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]; //計算出格點周圍4個點 if(isValue(g01) && isValue(g11)) { // All four points found, so interpolate the value. return bilinearInterpolateVector(i - fi, j - fj, g00, g10, g01, g11); } } } return null; }; function bilinearInterpolateVector(x, y, g00, g10, g01, g11) { //y表示格點距離上一個經緯度格點的距離 //x表示格點距離前一個經緯度給點的距離 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)]; //返回格點插值出來u和v還有速度 }; function distort(projection, λ, φ, x, y, scale, wind, windy) { //projection是一個空的物件 // λ, φ格點的經緯度 //x, y格點所在的畫素點 //格點的風向風速 //windy //scale 一個引數,每次粒子運動的距離 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; }; function distortion(projection, λ, φ, 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φ]; }; var project = function project(lat, lon, windy) { // both in radians, use deg2rad if neccessary var ymin = mercY(windy.south); var ymax = mercY(windy.north); var xFactor = windy.width / (windy.east - windy.west); var yFactor = windy.height / (ymax - ymin); var y = mercY(deg2rad(lat)); var x = (deg2rad(lon) - windy.west) * xFactor; var y = (ymax - y) * yFactor; // y points south return [x, y]; }; //第三步初始化生成粒子/////////////////////////////////// randomize = function(o) { // UNDONE: this method is terrible var x, y; var safetyNet = 0; do { x = Math.round(Math.floor(Math.random() * bwidth)); y = Math.round(Math.floor(Math.random() * bheight)); } while (field(x, y)[2] === null && safetyNet++ < 30); o.x = x; o.y = y; return o; }; function field(x, y) { var column = columns[Math.round(x)]; return column && column[Math.round(y)] || NULL_WIND_VECTOR; } function create_new_lizi() { particles.length = 0; var particleCount = Math.round(bwidth * bheight * PARTICLE_MULTIPLIER); //計算粒子的數量 for(var i = 0; i < particleCount; i++) { particles.push(randomize({ age: Math.floor(Math.random() * max_age) })); } } ////////////////////// //第四步,定義每次粒子的變化 var colorStyles = windIntensityColorScale(min_color, max_color); var buckets = colorStyles.map(function() { return []; }); function windIntensityColorScale(min, max) { colorScale.indexFor = function(m) { // map velocity speed to a style return Math.max(0, Math.min(colorScale.length - 1, Math.round((m - min) / (max - min) * (colorScale.length - 1)))); }; return colorScale; } function evolve() { //清除風點陣列 buckets.forEach(function(bucket) { bucket.length = 0; }); //particles是存放所有點的陣列 particles.forEach(function(particle) { //如果粒子執行的次數大於設定的最大次數 //重新隨機粒子的位置,並把次數定義成0 if(particle.age > max_age) { 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) { //如果沒有速度,就讓age是最大值,重新隨機點 particle.age = max_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; }); } ////////////////////// //第五步,初始化畫布,並畫粒子 var canvas = L.DomUtil.create("canvas", "can"); canvas.height = bheight; canvas.width = bwidth; // var wind_pane = map.getPane("markerPane"); // // canvas.style.background = "red"; // wind_pane.appendChild(canvas); var animated = map.options.zoomAnimation && L.Browser.any3d; //新增下面的class之後,圖層可以隨著地圖縮放變化 L.DomUtil.addClass(canvas, 'leaflet-zoom-' + (animated ? 'animated' : 'hide')); //把畫布新增到overlayPane圖層中 map._panes.overlayPane.appendChild(canvas); var fadeFillStyle = "rgba(0, 0, 0, 0.97)"; var g = canvas.getContext("2d"); g.lineWidth = linewidth; g.fillStyle = fadeFillStyle; g.globalAlpha = 0.6; map.on("zoomanim",function(e){ var scale =map.getZoomScale(e.zoom); var bound = map.getBounds(); var yun_lat1 = bound._northEast.lat; var yun_lon1 = bound._northEast.lng; var yun_lat2 = bound._southWest.lat; var yun_lon2 = bound._southWest.lng; var new_position = map.latLngToLayerPoint([yun_lat1, yun_lon2]); // -- different calc of offset in leaflet 1.0.0 and 0.0.7 thanks for 1.0.0-rc2 calc @jduggan1 var offset = L.Layer ? map._latLngToNewLayerPoint(map.getBounds().getNorthWest(), e.zoom, e.center) : map.getCenterOffset(e.center).multiplyBy(-scale).subtract(map.getMapPanePos()); L.DomUtil.setTransform(canvas, offset, scale); }) map.on("moveend", function() { window.cancelAnimationFrame(animationLoop); bwidth = map.getSize().x, bheight = map.getSize().y; g.clearRect(0, 0, bwidth, bheight); canvas.height = bheight; canvas.width = bwidth; var bound = map.getBounds(); var yun_lat1 = bound._northEast.lat; var yun_lon1 = bound._northEast.lng; var yun_lat2 = bound._southWest.lat; var yun_lon2 = bound._southWest.lng; var new_position = map.latLngToLayerPoint([yun_lat1, yun_lon2]); L.DomUtil.setPosition(canvas, new_position); bulid_grid(); create_new_lizi(); setTimeout(function() { var then = Date.now(); (function frame() { animationLoop = requestAnimationFrame(frame); var now = Date.now(); var delta = now - then; if(delta > FRAME_TIME) { then = now - delta % FRAME_TIME; evolve(); draw(); } })(); },200) }) function draw() { // Fade existing particle trails. var prev = "lighter"; g.globalCompositeOperation = "destination-in"; g.fillRect(0, 0, bwidth, bheight); g.globalCompositeOperation = prev; g.globalAlpha = galpha; // Draw new particle trails. //buckets是把風點按照不同顏色分級,分成的陣列 //陣列的每一項是一個物件, buckets.forEach(function(bucket, i) { if(bucket.length > 0) { g.beginPath(); g.strokeStyle = colorStyles[i]; bucket.forEach(function(particle) { g.moveTo(particle.x, particle.y); g.lineTo(particle.xt, particle.yt); particle.x = particle.xt; particle.y = particle.yt; }); g.stroke(); } }); } })