技術 | 使用Mapbox做AQI空氣質量分布圖

4 人贊了文章

空氣污染指數（AQI），就是根據環境空氣質量標準和各項污染物對人體健康、生態、環境的影響，將常規監測的幾種空氣污染物濃度簡化成為單一的概念性指數值形式。

使用工具：

1. Node.js；

2. Node.js插件：Mapnik、Turfjs、proj4等；

3. QGIS；

IDW插值

IDW(Inverse Distance Weighted)是一種常用而簡便的空間插值方法，它以插值點與樣本點間的距離為權重進行加權平均，離插值點越近的樣本點賦予的權重越大。一般使用IDW插值方法來分析空氣質量指數的空間變化趨勢。

Turfjs自帶標準的IDW插值演算法turf.IDW()，但主要針對矢量GeoJSON數據，同時為了方便後面參數設置，我們在這裡自定義了一個IDW方法。主要代碼如下：

數據裁剪

我們使用一個全國的多邊形圖層chinaregion.geojson，直接在IDW插值過程中對數據進行裁剪，忽略非多邊形區域外像素的計算，從而減少計算量。

此處，裁剪的計算方法使用turf.inside()來判斷判斷像素是否在多邊形內部：

var isInside = turf.inside(thecell, chinaareajson.features[0]);//true or false

註：由於後面每一個時間步長的數據，在計算時都是固定的空間範圍和像素分比例，因此我們可以使用turf.inside()預處理生成一個固定範圍圖片，來保存所要計算的像素。這樣每一次計算僅使用圖片中合法像素來進行計算。

var img = new mapnik.Image(grid.columns, grid.rows, { type: mapnik.imageType.gray8});......if (isInside) img.setPixel(column, grid.rows - 1 - row, 100);......img.save(theimage.tif, tif, function(err) { if (err) throw err;});

數據配色&導出

一般來說，空氣污染指數的取值範圍定為0～500，其中0～50、51～100、101～200、201～300和大於300，分別對應國家空氣質量標準中日均值的 I級、II級、III級、IV級和V級標準的污染物濃度限定數值。我們對計算後的全國空氣質量指數灰度圖按上述等級進行配色，並導出為PNG圖片。

優化IDW插值

由上圖可知，使用標準的IDW進行插值，對於沒有監測點覆蓋的大範圍區域（如上圖中新疆、西藏交接區域），受最近監測站數據的影響，插值結果比較大，不符合實際情況。

為了剔除測站數據對較遠區域的影響，我們將插值範圍限制在一定的範圍內，比如100KM、150KM：

var radus = 150;//150 km......if (d <= radus) { var w = 1.0 / Math.pow(d, b); sw += w; zw += w * controlPoints.features[j].properties[valueField];}

計算過程中，可以對測站點數據創建空間索引，加快計算速度。使用geokdbush來創建空間索引，並計算過程中，僅考慮一定半徑Radus內的測站點：

var index = kdbush(points, (p) => p.lon, (p) => p.lat);......var nearest = geokdbush.around(index, thecell[0], thecell[1], radus);for (var l = 0; l < nearest.length; l++) { var d = nearest[l].dist; if (d < mindis) mindis = d; var w = 1.0 / Math.pow(d, b); sw += w; zw += w * points[nearest[l].left].val;}

由上圖可知，雖然剔除了測站數據對較遠區域的影響，但由此造成結果數據「突變」嚴重，空氣質量指數數據擴散不連續，不符合實際情況。為了優化空氣質量指數擴散的連續性，我們在IDW插值過程中，對外側一定範圍內（100km- 150km）的數據進行再處理：

var secondradus = 100;//100 km, secondradus < radus......if (d < mindis) { mindis = d;}...... thevalue = zw / sw;var alpha = 255;if (mindis > secondradus){ thevalue = thevalue * (1 - (mindis - secondradus) / (radus - secondradus)); alpha = parseInt(255 - (mindis - secondradus) / (radus - secondradus) * 255);}

通過Mapbox GL JS載入渲染

Mapbox GL JS 是一個 JavaScript 庫，使用 WebGL 對矢量切片（vector tiles）和 Mapbox樣式（styles），構成的互動式地圖進行渲染。在這裡我們使用Mapbox GL JS將上述生成的空氣質量指標圖載入並渲染到網頁地圖中。為了有一個比較美觀的效果，你可以通過Studio生成並發布一張地圖作為底圖。

載入底圖

mapboxgl.accessToken = pk.xxxx.xxxx;//access tokenvar map = new mapboxgl.Map({ container: map, style: mapbox://styles/mapbox/streets-v9, center: [106, 36], zoom: 3.6, touchZoomRotatez: false, attributionControl: false});

載入數據：

將數據作為Raster類型載入到地圖中。

map.addLayer({ id: AQI + i, source: { type: image, url: data/low/a + (i + 1) + .png, coordinates: [ [73.50, 53.55], [134.80, 53.55], [134.80, 18.1775], [73.50, 18.1775] ] }, type: raster, paint: { raster-opacity: 0.7, raster-opacity-transition: { duration: 0 } }}, place-village);

插值方法除了在空氣質量指數AQI方面的應用以外，在其他領域也應用廣泛，例如區域氣溫分布、金屬元素地理統計分析等，歡迎大家在各自的領域使用Mapbox的技術。

除了使用插值方法來分析空氣質量指數AQI的分布以外，目前越來越多的人使用Heatmap熱點圖的方式來顯示空氣質量分布，Mapbox最近使用了Mapbox GL技術直接在客戶端使用Heatmap的方式渲染數據點值和密度的分布。詳情請戳這裡。