用spmap看中国空气质量

本文作者：赵冰洁, 中南财经政法大学金融学院

本文编辑：王玉洁

技术总编：戴   雯

Stata&Python云端课程来啦！

为了感谢大家长久以来的支持和信任，爬虫俱乐部为大家送福利啦！！！Stata&Python特惠课程双双上线腾讯课堂～原价2400元的Python编程培训课程，现在仅需100元，详情请查看推文《Python云端课程福利大放送！0基础也能学~》。关于Stata，爬虫俱乐部推出了系列课程，内容包括字符串函数、正则表达式、爬虫专题和文本分析，可以随心搭配，价格美丽，物超所值，更多信息可查看Stata系列推文《与春天有个约会，爬虫俱乐部重磅推出正则表达式网络课程！》、《与春天有个约会，爬虫俱乐部重磅推出基本字符串函数网络课程》等。变的是价格，不变的是课程质量和答疑服务。对报名有任何疑问欢迎在公众号后台和腾讯课堂留言哦

导读

前段时间，小编正津津有味地在网上冲浪，突然发现北京的同学在朋友圈里，讨论起了沙尘暴附体事件，有的说误以为来到了碟中谍4片场，有的说出门一趟2层口罩全黄。

的确，空气质量问题一直以来都是大家最关注列表中的一员，它与我们的生活息息相关。空气质量优时，我们是“自由的”，可以不用担心是否会吸入有害颗粒，是否会对正常出行产生影响；空气质量差时，我们是“被禁锢的”，为了避免有害颗粒的伤害，我们会减少出行，整座城市都被蒙上了“神秘”的袈裟。

今天，敬业的小编将从数据的角度，为大家展示3月份中国的空气质量状况。来吧，上才艺~

01空气质量指数（AQI）爬取

本次爬虫所需数据均来自于空气质量数据库(http://pm.fuhaodaquan.org/)。该数据库由首师大李教授提议，并由首师大李教授空气质量实验和中科院联合赞助设立，其中的空气质量数据主要来自环保部实时监测平台。不同等级的空气质量对应不同的配色，绿色代表一级空气质量（优）、黄色代表二级空气质量（良）、橙色代表三级空气质量（轻度污染）、红色代表四级空气质量（中度污染）、紫色代表五级空气质量（重度污染）、褐红色代表六级空气质量（严重污染）。

首先，我们进入空气质量数据库中的全国空气质量指数（AQI）月均值排行榜，打开的页面是这样的：

那么，该如何获取排行榜中所有城市的空气质量指数（AQI）呢？

打开网页源代码，找到我们所要提取的有效信息的源代码位置，如下（红色标框部分）：

接下来，使用copy命令获取网页源代码，转码，并提取源代码信息。

clear all
cap mkdir "D:\Stata"
cd "D:\Stata"
copy "http://pm.fuhaodaquan.org/paihangbang/paihangbang_mon.php?by=aqi&order=desc&area=%C8%AB%B9%FA" temp.txt, replace
infix strL v 1-100000 using temp.txt, clear
replace v = ustrfrom(v, "gb18030", 1)
keep if index(v, "<tr><td>")
drop in 1
split v, p("<td><a ")
sxpose, clear
drop in 1/2
split _var1, p("</td>")
drop _var1 _var13

运行结果如下：

然后，利用正则表达式提取各省省份名称和空气质量指数（AQI）两个有效数据，再进行相应的数据整理以便后续绘图。

** 删除所有 <> 中的内容
foreach v of varlist _all {
replace `v' = ustrregexra(`v', "<.*?>", "")
replace `v' = ustrregexra(`v', ".*?>", "")
}
rename _var11 city
rename _var12 AQI

** 提取各省的省份名称以及相应的城市名称
gen prov = ustrregexs(1) if ustrregexm(city, "\((.*?)\)")
replace prov = city if prov == ""
drop city
order prov, before(AQI)

** 在空气质量数据库中，“安徽省”被命名为了“安微省”，为了后续数据合并，这里用字符串函数进行省份名称的调整
replace prov = subinstr(prov, "安微", "安徽", .)
destring AQI, replace

** 生成各城市的空气质量指数AQI均值，并取整数
sort prov AQI
by prov: egen v1 = mean(AQI)
replace AQI = v1
drop v1
duplicates drop
replace AQI = int(AQI) 

** 保存数据
compress
save 全国空气质量指数（AQI）3月均值.dta, replace

运行结果如下：

02空气质量指数（AQI）数据可视化

为了更加直观的观测到中国空气质量状况，我们引入了这个空间数据可视化命令—— spmap。

2.1 spmap原理介绍

在绘制地图之前，先进行命令的安装：

ssc install spmap
ssc install shp2dta
ssc install mif2dta

绘制地图一般需要一份带标签的数据和一份地理坐标数据。首先，我们先下载一份中国的省级地图并解压压缩文件：

clear all
copy "http://fmwww.bc.edu/repec/bocode/c/china_map.zip" china_map.zip, replace 
unzipfile china_map, replace  

在文件解压缩之后，我们的文件夹里会有三份文件：标签数据china_label.dta 、地理坐标数据 china_map.dta以及一份绘图示例的do文件 map_example.do。

在打开标签数据china_label.dta之后，我们会发现部分中文字符会出现乱码的情况，这时就需要进行相应的转码：

use china_label.dta, clear
replace name = ustrfrom(name,"gb18030",1)
save china_label.dta, replace

在完成命令安装以及数据处理之后，我们来简单了解一下绘制地图的原理。分别打开标签数据和地理坐标数据：

我们可以发现，标签数据和地理坐标数据是通过变量 id （或者 _ID ）连接的。

标签数据中，id=1时，对应黑龙江省；地理坐标数据中，ID=1对应的X、Y即为黑龙江省的地理坐标，由这两个坐标所围成的图形即为黑龙江省的地理图。

因此，Stata绘制地图的原理就是通过标签数据中的变量id投射到地理坐标数据中其相对应的地理坐标所围成的地理图形上。简单来说，就是先通过标签数据中的变量id确定省份，再在地理坐标数据中找到与变量id相对应的坐标数据，然后根据坐标数据绘制地图。

在绘制地图之前，需要将前面爬取到的数据全国空气质量指数（AQI）3月均值.dta 与 标签数据china_label.dta 按照省份进行合并。在合并之前，我们会发现两份数据文件中省份不一致，需要进行调整。

** 用正则表达式调整标签数据中的省份名称以使两份数据文件中的省份名称一致
use china_label.dta, clear
replace name = ustrfrom(name,"gb18030",1)
gen prov = ustrregexra(name, "省|市|自治区|回族|维吾尔|壮族|特别行政区","")
order prov, after(name)
save china_prov.dta, replace

** 按照省份prov进行合并
use 全国空气质量指数（AQI）3月均值.dta, clear
merge m:m prov using china_prov.dta
list if _merge == 2
drop if _merge != 3
drop _merge

** 保存数据
save 1.dta, replace

运行结果如下：

完成数据合并之后，便可绘制地图啦。

2.2 绘制地图

spmap 命令是通过叠加三个图层来实现绘图的：

• 第一层：底层地图（basemap）。其中，选项 id(var) 为必选项，其用来识别地图中不同区域所划分的_ID。其他的选项用来控制地图的另外四个方面：Cartogram （底图的背景色等）、Choropleth map（将作图数据分为几组、分组方式等）、Format（地图的填充颜色等）、Legend（是否显示图例、图例主题、图例标签等）。

• 第二层：在底图的基础上，添加其他对象。polygon（添加其他多边形）、line（添加线条）、diagram（添加条形/扇形统计图）、arrow（添加箭头）、label（添加标签）等。

• 第三层：在第一层、第二层的基础上添加元素进一步美化地图。比如 title（标题）、note（注释）、scalebar（比例尺）等。

  spmap AQI using "china_map.dta", ///
  id(id)  clmethod(custom)  ///
  clbreaks(0 50 100 150 200 300 500) fcolor(green yellow orange red purple maroon) ///
  legenda(on)  legorder(lohi)  ///
  graphr(margin(medium))  ///
  legend(title("AQI", size(*0.5) pos(11) color(black))  ///
         order(2 "优" 3 "良" 4 "轻度污染" 5 "中度污染" 6 "重度污染" 7 "严重污染") color(black))  ///
  label(label(prov) xcoord(x_coord) ycoord(y_coord))  ///
  title(中国空气质量图--3月均值, size(*0.8))  /// 
  subtitle("2021年03月01日-2021年03月31日", size(*0.8))  ///
  scalebar(units(500) scale(1.0) xpos(100) label("全国空气质量指数（AQI）"))

整体来看，三月份全国的空气质量中等，北方地区（除黑、吉、辽外）较南方地区而言，空气质量较差，其中甘肃、宁夏两个省份最为严重，海南省空气质量最佳。

2.3 绘制地图进阶版

带气泡的地图（Bubble Map）：可以通过添加 point()  选项完成带气泡的地图的绘制。

spmap using "china_map.dta", id(id) fcolor(white)  ///
  legenda(off) ///** 第二层：在底图的基础上，加上气泡图形，同时添加标签 **/
  point(data(1.dta) xcoord(x_coord) ycoord(y_coord) proportional(AQI) fcolor(green) ocolor(white) size(*2))  ///
  label(data(1.dta) label(AQI) xcoord(x_coord) ycoord(y_coord) color(white) size(*0.7))  ///** 第三层：进一步美化，添加标题、副标题、比例尺 **/
  title(中国空气质量图--3月均值, size(*0.8))  /// 
  subtitle("2021年03月01日-2021年03月31日", size(*0.8))

附加条形图的地图：可以通过diagram()选项完成附加条形图的地图的绘制。

** 第一层：底层地图 **
spmap using "china_map.dta", id(id) fcolor(stone) ///
  legenda(off)  ///** 第二层：在底图的基础上，添加条形图
  diagram(data(1.dta) variable(AQI) range(0 151) xcoord(x_coord) ycoord(y_coord) fcolor(red) size(1.2))  ///** 第三层：进一步美化，添加标题、副标题、比例尺 **/
  title(中国空气质量图--3月均值, size(*0.8))  /// 
  subtitle("2021年03月01日-2021年03月31日", size(*0.8))

最后为大家呈现一个3月华南地区的空气质量图：

从图中可以看出，华南地区普遍空气质量较好，大部分地区空气质量在[50,100]的范围内，甚至还有一些地区空气质量指数在[0,50]，即空气质量十分优良。在一个北方被蒙古黄沙席卷的天气里，华南地区的小伙伴们依然可以愉快的享受美好的天气，不得不说空气质量受地理位置的影响也是很大的呀

Tips：由于本期推文所涉及到的数据库中未包含台湾的空气质量指数，所以台湾的AQI是通过其他方式获取到的；其次，由于中国地图九段线部分的数据缺失，所以为保证中国地图的完整性，小编通过其他方式将其补充完整。

以上就是本期推文为大家呈现的有关全国3月份空气质量地图的全部内容啦，后台回复spmap可以将相关数据发给大家哦~

希望大家都可以在阳光明媚的天气里，遥望蓝天和云朵。保持一整天的好心情~

欢迎大家点赞、转发、评论哦！

对我们的推文累计打赏超过1000元，我们即可给您开具发票，发票类别为“咨询费”。用心做事，不负您的支持！

往期推文推荐         OpenCV库——轻松更换证件照背景色
这些年，爬虫俱乐部走出的博士们！
         看这里，近五年各校高被引论文上榜名单！

关于我们 

微信公众号“Stata and Python数据分析”分享实用的Stata、Python等软件的数据处理知识，欢迎转载、打赏。我们是由李春涛教授领导下的研究生及本科生组成的大数据处理和分析团队。

此外，欢迎大家踊跃投稿，介绍一些关于Stata和Python的数据处理和分析技巧。

投稿邮箱：[e-mail]投稿要求：
1）必须原创，禁止抄袭；
2）必须准确，详细，有例子，有截图；
注意事项：
1）所有投稿都会经过本公众号运营团队成员的审核，审核通过才可录用，一经录用，会在该推文里为作者署名，并有赏金分成。
2）邮件请注明投稿，邮件名称为“投稿+推文名称”。
3）应广大读者要求，现开通有偿问答服务，如果大家遇到有关数据处理、分析等问题，可以在公众号中提出，只需支付少量赏金，我们会在后期的推文里给予解答。