十行代碼看到空氣質量指數

作者：統計之都 專業、人本、正直的中國統計學門戶網站

個人公眾號：統計之都

出處：十行代碼看到空氣質量指數

故事部分

我錯了, 我承認我是標題黨, 怎麼可能用十行代碼完成全國三百個多個城市空氣質量指數的抓取, 清洗與可視化呢！

我仔細數了數, 去掉注釋, 一共是9行, 湊個整才是10行 耶~

空氣質量指數（英文全稱Air Quality Index，簡稱AQI）是定量描述空氣質量狀況的無量綱指數.

關於空氣質量的段子已經層出不窮， 連呆在上海的我都已經開始關注北京的天氣了:

受朋友委託，大家幫個忙：北京人，女，26歲，未婚，1.68米，體重50公斤，英國海歸。貌美，愛好健身。目前在一家世界500強做產品經理，工作穩定，年薪近90萬。三環內有四套學區房，一套按揭，三套全款。 名下有一輛寶馬7系，上班時開。父母均是國家領導幹部。朋友和家人現在非常著急，想讓介紹一個效果比較好的防霾口罩。

正常不過三秒, [手動債見]

咳咳, 回歸正題, 代碼之前, 先介紹下出場的一個新包,leafletCN.

leafletCN【1】 是一個基於leaflet【2】包做的大中華擴展, 優勢在於有細分到縣級市級別的區劃數據, 地理區劃數據來源於github的geojson-map-china【3】項目。雖然沒那麼准, 但是也基本能用了~

AQI的數據來源於pm25.in, 網頁上是一個html的表格, 可以簡單的用XML的readHTMLtable來完成讀取.

好了, 以下是代碼, 分成三個部分：

1. 載入包, 並讀取網頁的數據

2. 整理數據並進行命名, 包含了獲取其中的城市, AQI以及將其轉化為數值

3. 使用geojsonMap繪製細分到城市的污染情況

這部分代碼運行之後會出現一個以高德地圖為底圖, 可縮放可點擊的地圖~

library(XML)nlibrary(leafletCN)n# 讀取網頁的表格n# Sorry for 爬了你家網站ntable = readHTMLTable("http://www.pm25.in/rank", n encoding = "UTF-8", stringsAsFactors = F)[[1]]nn# 整理數據並命名ndat = table[ , 2:3]nnames(dat) = c("city","AQI")ndat$AQI = as.numeric(dat$AQI)nn# 調用geojsonMap進行繪製ngeojsonMap(dat, "city",n popup = paste0(dat$city, ":", dat$AQI),n palette = "Reds", legendTitle = "AQI")n

猛擊圖片查看這個可交互的地圖吧！

安裝

## 穩定版n install.packages("leafletCN") n## 開發版 ndevtools::install_github("lchiffon/leafletCN")n

基本使用

常用的函數

• regionNames 返回某個地圖的區域名

• demomap 傳入地圖名繪製示例地圖

• geojsonMap 將一個分層顏色圖繪製在一個實時地圖上

• amap 在leaflet地圖上疊加高德地圖

• read.geoShape 讀取一個geojson的對象,保存成spdataframe,以方便leaflet調用

• leafletGeo 用地圖名以及一個數據框創建一個sp的對象

regionNames

傳入需要查看的城市名, 顯示這個城市支持的區域信息, 比如查看成都:

regionNames("成都") n[1] "成華區" "崇州市" "大邑縣" "都江堰市" "金牛區" n[6] "金堂縣" "錦江區" "龍泉驛區" "彭州市" "蒲江縣" n[11] "青白江區" "青羊區" "雙流縣" "溫江區" "武侯區" n[16] "新都區" "新津縣" "邛崍市" "郫縣"n

如果不傳入對象, 會自動返回300多個支持的名字列表,包括各個城市,省,以及三個特殊的名字:

1. world 世界地圖

2. china 中國分省份地圖

3. city 中國分城市地圖

demomap

傳入城市名(省名),顯示這個城市(省)的示例地圖

demomap("台灣")n

geojsonmap

將一個數據框顯示在需要展示的地圖上. 在函數中做了一些有趣的設置, leafletCN會自動匹配傳入的前兩個字元來尋找合適的位置進行繪製, 所以基本不需要糾結是寫』上海市』還是』上海』了。PS：不用擔心張家口和張家界這對名義雙胞胎，因為遇到這種情況默認會重查~~：）

圖做出來可以在上面點點點…

dat = data.frame(name = regionNames("china"),nvalue = runif(34)) ngeojsonMap(dat,"china")n

amap

amap為leaflet疊加一個高德地圖, 使用:

leaflet() %>% namap() %>%naddMarkers(lng = 116.3125774825, lat = 39.9707249401, n popup = "The birthplace of COS")n

read.geoShape

read.geoShape這個函數可以把一個geojson格式的數據讀取為一個SpatialPolygonsDataFrame對象, 方便sp或者leaflet包中的調用.

if(require(sp)){n filePath = system.file("geojson/china.json",package = "leafletCN")nmap = read.geoShape(filePath) n plot(map)n}n

leafletGeo

leafletGeo這個函數可以把一個數據框和一個地圖組合在一起, 方便用leaflet調用, 其中名字的 變數為name, 數值的變數為value~

if(require(leaflet)){n dat = data.frame(regionNames("china"),n runif(34))n map = leafletGeo("china", dat)nn pal <- colorNumeric(n palette = "Blues",n domain = map$value)nn leaflet(map) %>% addTiles() %>%n addPolygons(stroke = TRUE,n smoothFactor = 1,n fillOpacity = 0.7,n weight = 1,n color = ~pal(value),n popup = ~htmltools::htmlEscape(popup)n ) %>%n addLegend("bottomright", pal = pal, values = ~value,n title = "legendTitle",n labFormat = leaflet::labelFormat(prefix = ""),n opacity = 1)n}n