大气污染物浓度风域空间分布图的绘制方法

杨洛祺 顾世杰 刘晓帅 胡引翠 武帅 田冰

摘要：风玫瑰图主要用以表达风向频率或各方向上的平均风速。若要可视化表达各风向上不同风速情境下的大气污染物浓度的分异特征，需要突破当前风玫瑰图仅展示风速风向信息的限制。为解决这一问题，文章在风玫瑰图基础上，采用晕渲法，用深浅不同的色调表示研究区域在统计时段内大气污染物浓度在各风向、不同风速下的特征。基于Python，使用contour函数以及numpy、matplotlib、pandas库，文章研发了大气污染物浓度风域空间分布图绘制工具，解决了绘图中由于数据缺失而造成的图像空白和计算报错的问题，并形成创新性的污染物质量浓度与风域的组合晕渲玫瑰图绘制方法。相较于传统的气象学制图软件，该工具绘图基础要求偏向底层绘图，更加灵活易用，为研究大气污染物区域特征以及跨区域传输特征提供了新的可视化方法和技术手段。

关键词：大气；污染物质量浓度；风玫瑰图；风域；空间

中图分类号：P412.36；P208      文献标识码：A

文章编号：1009-3044（2023）06-0113-05

开放科学（资源服务）标识码（OSID）：

0 引言

在大气污染物风域空间分布研究中，定量解析及污染物溯源研究是目前研究的重点方向之一[1-5]。大气污染物浓度分布可视化表达是辅助发现大气污染传输区域特征及研究成果共享交流的重要技术手段，专题地图和量化图表是常用的可视化表达方式。

专题地图在基于底图的基础上按照主题的要求，突出并完善地表示一种或几种与主题相关的要素信息，着重描述专题化的内容及专门化的用途[6]。通常采用地理信息系统软件分析制作大气污染专题地图[7-9]，如绘制大气污染物后向轨迹图，分析污染物潜在源区及轨迹分布图；采用办公语言或高级编程语言绘制量化图表[10-13]，如NCL绘制污染物质量浓度风玫瑰图。量化图表在量化信息的情况下突出一种或几种特定表达指标的价值和衡量标准。量化图表中，风玫瑰图是得到广泛应用的方法，其直观地描述了风速风向等气象要素的变化特征。其将地面风向分为8个或16个方位进行表示，用极坐标系绘制某一时间段内研究区域风速风向出现频率的综合统计图。其在建筑规划、风能资源评估、大气污染评价等领域具有广泛的应用[11， 14]。张顺尧[15]等在城市微气候对建筑物的影响研究中，通过小区域风玫瑰图的运用，揭示了建筑物开口方向与微气候中风场的相关关系；张克存[16]等利用风玫瑰图可视化了测点风向及起沙风频率，辅助性地展示了湖泊外围的输沙方向；张宇[17]利用风玫瑰图与Pearson系数表达了臭氧和VOCs与风速风向的相关性程度及拟合关系。

目前绘制风玫瑰图以C#、Python、NCL以及Extjs、Excel、Matlab为主要的高级编程语言和办公语言为基础[18-24]。风玫瑰图用以表达风向频率或各方向上的平均风速。本文选取量化图表的表示方法，分析研究区域各方向风速大小并加和污染物浓度信息要素，获得多级地理尺度的大气污染物时序变化数据集，构建研究区域污染物组合污染物质量浓度与风域的组合玫瑰图[25-28]。

目的是可視化表达各风向上不同风速情境下的大气污染物浓度的分异特征，突破当前风玫瑰图仅展示风速风向信息的限制。为解决这一问题，本文在风玫瑰图基础上，采用晕渲法，用深浅不同的色调表示研究区域在统计时段内大气污染物浓度在各风向、不同风速下的特征。

将依托风玫瑰图可视化展示的特点，在风向、风速所构成的风域空间中，晕渲呈现站点监测的大气污染物浓度在风域空间中的变化。为实现这一可视化技术思路，本文基于Python语言中anaconda3[29]配置环境并使用第三方库Simple Index库，统筹使用晕渲法及叠加法的原理，设计一种融合大气污染物（PM2.5、PM10、臭氧等）浓度、三维风向包括水平分量及垂直分量及风速信息的污染物质量浓度与风域的组合晕渲玫瑰图的绘制工具，绘制结果能更清晰地体现污染物的浓度变化及来源信息，为近地面大气污染物传输特征的研究提供一种更直观的可视化方法。

1 数据来源与基础数据预处理

1.1 数据来源

文章使用的数据包括：日期、风速、风向和污染物浓度。

数据来源主要为验证实例的污染物近地面监测点、气象站点的实时监测数据，欧洲哥白尼卫星数据服务中心提供的ERA-5数据和清华大学提供的TAP数据。

近地面监测点数据来源于河北师范大学地理科学学院监测点（37.99°N，114.52°E），该位置属于城市区域次交通枢纽和科研混合区。气象数据来源于石家庄市气象站（53698站）数据，包括实时风速、风向等信息。本研究中利用在线观测仪器测量了PM2.5逐时质量浓度及所使用的气象数据信息，测量时间为2022年8月1日。

ERA-5数据（https：//cds.climate.copernicus.eu/-！/home）来源于欧洲哥白尼卫星数据服务中心提供的数据集，选取石家庄市（38.05°N ，114.52°E）近地面10m的水平方向风速分量和垂直方向风速分量的数据；TAP（Tracking Air Pollution in China）数据（http：//tapdata.org.cn/）来源于清华大学提供的近地面1km PM2.5浓度数据，该数据综合了多种气溶胶数据反演及站点数据，自然源及人为源数据作为辅助数据，较高水平地反映了当时PM2.5浓度具体情况。

1.2 近地面监测点数据预处理

近地面监测点和气象站监测到的大气污染物浓度、风场等数据的储存均为文本格式，需将数据进行储存格式转换，并以污染物类型和具体时间为依据对基础数据进行统计，将Excel表格上传到PyCharm软件中加载数据。

首先，将污染物浓度的数据类型浮点型（float）转换为短整型（int），减少存储空间，提高计算效率。

其次，将数据由文本型转换为Excel格式，并以污染物类型和时段为依据对基础数据进行分类。

再次，将预处理后的数据（包括风向、风速和PM2.5浓度）加载至工具软件，将风向单位由度转换为弧度，将数据缺失的网格用0值进行填充。

最后，利用pivot_table函数绘制数据透视表，使用meshgrid函数建立相关矩阵。主要代码如下：

data['PM2.5'] = data['PM2.5'].astype（float）

dt=data.pivot_table（values='PM2.5'， index='风速（m/s）'， columns='风向（deg）'， aggfunc=np.mean）

dt.fillna（0， inplace=True）

dt = dt.reindex（index=speed， columns=deg， fill_value=0）

dt.head（31）

theta， r = np.meshgrid（deg， speed）

本程序案例应用日均数据绘制日均大气污染物浓度风域空间分布图，但该程序不仅限于绘制日均图，可根据自身研究需要，绘制月均或年均污染物风域空间分布图。在应用中需注意修改#header函数，使函数参数与Excel格式对齐，实现软件正常读取。例如小时数据需修改为#header=24；逐月数据为#header=12等。

1.3 ERA-5数据预处理

下载对应研究区域的ERA-5数据，对数据进行格式转换。数据中包括近地面10m的水平方向风速分量u-wind和垂直方向风速分量v-wind，可通过程序代码实现水平、垂直方向的风速风量和风向的转换。主要代码如下：

if （u > 0） and （v > 0）：

fx = 270 - math.atan（v / u） * 180 / math.pi

elif （u < 0） and （v > 0）：

fx = 90 - math.atan（v / u） * 180 / math.pi

elif （u < 0） and （v < 0）：

fx = 90 - math.atan（v / u） * 180 / math.pi

elif （u > 0） and （v < 0）：

fx = 270 - math.atan（v / u） * 180 / math.pi

elif （u == 0） and （v > 0）：

fx = 180

elif （u == 0） and （v < 0）：

fx = 0

elif （u > 0） and （v == 0）：

fx = 270

elif （u < 0） and （v == 0）：

fx = 90

elif （u == 0） and （v == 0）：

fx = 999.9

2 可视化实现

本文开发环境为PyCharm Community Edition2021.2。由于Python3不兼容Python2，所以在PyCharm Community Edition 2021.2中将基础环境设定为anaconda环境并安装Simple Index开源环境库（https：//www.lfd.uci.edu/～gohlke/pythonlibs/）由北京外国语大学提供。

利用传统的Windrose绘图工具绘制风玫瑰图，具有只包含风速风向的信息、操作复杂、易出现统计结果报错等局限性，为后续的专题化相关研究带来了困难。为解决这一问题，本文将渲染法和透视法创新性与风玫瑰图相结合，主要使用contour函数、numpy库、maplotlib库及pandas库进行绘制方法的创新，为风域空间专题呈现提供新思路。

绘制过程中，对风速风向数据进行统计，利用linspace函数将数据进行分类，其中风速按16等分分类，风向按32等分分类。风向信息以圆形来表示，作为基础底图。绘制出的0°方向为正北方向，风向度数对应正北方向与风向间的夹角并按顺时针方向测量，范围为0°～360°，单位为角度数。风速通过不同大小的圆进行表示，根据研究区域的风速统计值域范围设定圆的半径R（程序案例默认R=10，即所获数据中最大风速不超过10m/s），将数据进行网格化并把网格内的数据按较小值划分，利用def函數创建for循环，形成嵌套圆图形。主要代码如下：

data['deg'] = np.radians（data['deg']）

v = data['speed']

R=10

speed = np.linspace（0，R， endpoint=True， num=16）

deg = np.linspace（0， 2*np.pi， endpoint=True， num=32）

def maker（s， sequence）：

'''

Divide the data in the grid into smaller values

'''

fori， val in enumerate（sequence）：

if s <= sequence[i]：

returnval

d = data['deg']

data['speed'] = v.apply（maker， sequence=speed）

data['deg'] = d.apply（maker， sequence=deg）

data.head（）

ax.set_theta_zero_location（"N"）

ax.set_theta_direction（'clockwise'）

根据实际需要在subplot函数中设置画布，调整figsize参数；自定义设置底图颜色，具体颜色参照基础色变换，本案例选用16进制颜色。该步骤核心为使用matplotlib库的pyplot模块中的contourf函数绘制轮廓。添加右侧图例色带，根据所获得的污染物浓度值域设置色带范围，调整为闭合状态并利用cmap参数自定义色带颜色类型，该参数的候选值有：Accent、Accent_r、Blues等（程序案例默认色带为'CMRmap_r'）。主要代码如下：

fig = plt.figure（figsize=（8，6））

ax = plt.subplot（projection='polar'， facecolor="white"）

ax.set（ylim=（0，R），yticklabels=（[]））

pos=ax.contourf（theta，r，dt.to_numpy（），levels=np.linspace（0，150，endpoint=True，num=10） ，cmap='CMRmap_r'，）

cb=plt.colorbar（pos， ax=ax，pad=0.1）

为更清晰地展现风速，利用set_ylabel函数在图面左侧添加坐标轴用来标识风速并使用set_yticklabels来标注刻度，为避免程序报错，刻度值数量需与风速值数量相同。如需进一步美化图形，可利用plt.title函数设置标题，设置条带字体大小，坐标轴字体大小，字体样式等。Python里默认的编码格式是ASXCII格式，而中文是GBK格式，因此在PyCharm中输入带有中文的路径时会出现系统报错，现提出了可利用font manager函数来解决使用中文时出现乱码的问题。主要代码如下：

from matplotlib import font_manager as fm

plt.rcParams['font.sans-serif'] = ['Time New Roman']

plt.rcParams['axes.unicode_minus'] = False

font1 = fm.FontProperties（fname='C：＼Windows＼Fonts＼msyh.ttc'）

se=ax.secondary_yaxis（-0.1）

se.set_ylabel（'Speeed'+r'$＼mathrm{/（m·s^{-1}）}$'，font=font1）

se.set_yticks（[0，1，2，3，4，5，6，7，8，9，10]）

se.set_yticklabels（[-10，-8，-6，-4，-2，0，2，4，6，8，10]）

frommatplotlib.ticker import MultipleLocator

xminor = MultipleLocator（1）

se.yaxis.set_minor_locator（xminor）

se.tick_params（axis="y"， direction="in"， which="both"）

plt.title（'Combination Rose Chart of Average Wind and $＼mathrm{PM_{2.5}（μg·m ^{-3}）}$ Mass Concentration'）

cb.ax.set_title（r'$＼mathrm{PM_{2.5}（μg·m ^{-3}）}$'，pad=10，size=10）

調整后的效果显示为正确，对图形进行输出，案例输出结果如图1所示。

3 应用个例

文章案例选取2022年8月1日河北师范大学地理科学学院监测点及石家庄气象站点（53698站）的实时监测数据，监测距地1km的逐时PM2.5浓度和石家庄市近地面平均风速风向信息，并基于ERA-5数据的2021年1月石家庄市的局部风场信息及TAP数据按照编写的程序对个例数据进行绘制测试。

3.1 基于地面监测数据的PM2.5浓度风域空间分布图

程序支持日均、月均、年均等多种格式的Excel表格数据（包括风速m/s、风向°及污染物浓度值μg/m3），案例程序选择日均格式的Excel表格数据，所用数据为2022年8月1日00：00至8月31日23：00河北师范大学近地面监测点监测到的PM2.5污染物浓度值及石家庄市国家基础气象站（53894站）监测到的同天近地面风速风向数据，如表1所示。

通过绘制所选监测点的大气污染物浓度风域空间分布图（图3），对风向与污染物质量浓度进行相关性分析：2022年8月1日河北师范大学周围1km内总体AQI值较低，全天各时段PM2.5质量浓度均不足50μg/m3。其中风场是影响大气污染物累积和消散的关键性因素，不同的风速风向会影响大气中污染物的浓度和传输方向。当日监测点的主要风频为西南风，根据PM2.5污染物的分布位置可知，当日PM2.5主要聚集在水平风为偏南风和东南风的位置；PM2.5浓度较高的地方风速较小，最高风速不足4m/s，最低风速约为1m/s。风速不足使污染物缺少跨区域传输的动力，因此所选区域当天PM2.5的主要来源为当地自身排放，有少许污染物来自于河北师范大学南部方向。

3.2 基于ERA-5和TAP数据的PM2.5浓度风域空间分布图

实例应用选取2021年1月ERA-5的石家庄市月平均近地面10m的水平、垂直方向的风速风量数据，利用转化程序得到风向度数值，后导入到Excel表格中，格式要求与监测站相同。利用ERA-5数据中包含的风速数据、转换后的风向数据，对应的地点PM2.5浓度的TAP数值共同绘制大气污染物浓度风域空间分布图（图4）。由图4可知，2021年1月石家庄市PM2.5质量浓度较高，空气质量状况属于轻度污染，PM2.5浓度值最高达到了170μg/m3，均值在100μg/m3左右。污染物在该时间段主要聚集在风向为正南、西南以及正西方向，并在风速小于4m/s的位置出现堆积现象。当风速小于等于2m/s时，PM2.5质量浓度的数值最高，说明当月石家庄市PM2.5污染物主要来自本地排放；少许污染物来源于跨区域传输，且主要为偏南方向及西南方向传输来的污染物，其中受偏南方向的地区污染物的跨区域传输的影响较大。因此石家庄市该时段PM2.5浓度值升高主要受到邢台市和邯郸市的污染物跨区域传输的影响。

4 结束语

文章提出了一种大气污染物浓度风域空间分布图绘制方法，并研发了制图工具：1）制图中将原数据中的空值部分填补为0值，使制图结果更加平滑、美观。2）创新性的将傳统的风玫瑰图进行优化，结合渲染法和叠加法，将风速风向数据与污染物浓度数据进行融合，使风玫瑰图专题化，更加直观地体现了污染物的聚集位置和传输方向，并结合了具体的案例进行应用分析。但是目前制图工具尚存在一些不足，比如导入数据形式较为单一，暂时只支持Excel表格，需今后加以改进，提高程序加载数据的包容性。

参考文献：

[1] 李婷苑，陈靖扬，龚宇，等.2022年广东省冬季一次臭氧污染过程的气象成因及潜在源区分析[J/OL].环境科学.[2022-11-10].https：//kns.cnki.net/kcms2/article/abstract？v=aoPRB0Fc1 faWAu-EXeRzyMxUuj5GhO0pt8nbA6uDJI6fvWKAnxvMUkK Xlr_OzYQUlayrCHhpecA4CyHKtyU8Y3AbZwdNoPqmutJB6 3LNWZfibjA1EJ1L3g==&uniplatform=NZKPT.

[2] 王海宁，于大江，刘硕，等.龙凤山大气本底站CO浓度特征及传输路径分析[J].环境科学学报，2022，42（9）：392-400.

[3] 武威，单铁良.基于后向轨迹的秋冬季漯河重污染输送及典型个例分析[J].气象与环境学报，2022，38（3）：65-74.

[4] 谢放尖，郑新梅，窦焘焘，等.南京地区细颗粒物污染输送影响及潜在源区[J/OL].环境科学. [2022-11-10].https：//kns.cnki.net/kcms2/article/abstract？v=aoPRB0Fc1fZYNQy3XUdgTVvss 7uswCns2_qOPmgG0jy4UkMqNSG33Q7qRoj8t7zcAQHeptjxX zGZl12pMzHwcBZeTNe_JQkT3uJT2Jz0qsMSCKSa8pmv_A==&uniplatform=NZKPT.

[5] 杨芳园，潘娅婷，康道俊，等.2019年昆明市一次臭氧污染过程特征及成因分析[J].环境科学学报，2022，42（10）：71-79.

[6] 张军海，李仁杰，傅学庆.地理信息系统原理与实践[M].2版.北京：科学出版社，2015.

[7] 陈旭，邹滨，李沈鑫，等.个体大气污染暴露测量与可视化系统[J].地理与地理信息科学，2019，35（4）：21-26，56.

[8] 杜诗薇，陈庆涛，向竟德.基于ArcGIS Engine的桂林市七星区大气污染扩散模拟[J].城市地理，2017（10）：188-189.

[9] 王维，程灏旻，叶秀云，等.基于ArcGIS的2020年内江市城市环境空气质量特征分析[J].环境保护与循环经济，2021，41（10）：58-60.

[10] 徐伟嘉，尹萌，李红霞.基于GIS的珠三角区域空气质量实况发布平台介绍[J].中国环境监测，2015，31（4）：146-152.

[11] 张鑫，曹蕾，韩基良.基于Python气象数据处理与可视化分析[J].气象灾害防御，2020，27（1）：29-33.

[12] 赵奎锋.GrADS网络交互绘图技术及应用[J].陕西气象，2019（3）：48-50.

[13] Naghibi S A，Hashemi H，Pradhan B.APG：a novel python-based ArcGIS toolbox to generate absence-datasets for geospatial studies[J].Geoscience Frontiers，2021，12（6）：101232.

[14] 谢志祥，秦耀辰，郑智成，等.京津冀大气污染传输通道城市PM2.5污染的死亡效应评估[J].环境科学学报，2019，39（3）：843-852.

[15] 张顺尧，陈易.基于城市微气候测析的建筑外部空间围合度研究——以上海市大连路总部研发集聚区国歌广场为例[J].华东师范大学学报（自然科学版），2016（6）：1-26.

[16] 张克存，奥迎焕，屈建军，等.巴丹吉林沙漠湖泊-沙山近地表风沙动力环境[J].干旱区地理，2013，36（5）：790-794.

[17] 张宇.赤壁市臭氧污染特征及其前体物VOCs臭氧生成潜势分析[D].武汉：华中科技大学，2021.

[18] 阿不都克依木·热合曼.试探阿克陶县风玫瑰图的绘制[J].城市建设理论研究（电子版），2019（29）：53.

[19] 曹丽娟，周玲，陈艳丽，等.风玫瑰图的研究与程序自动成图设计[J].计算机与现代化，2013（2）：134-137，146.

[20] 陈宇罡，汪永青，李琳，等.Java、Python和Matlab混合編程及其在气象中的应用[J].电子世界，2017（16）：15-16.

[21] 褚红健，李佑文，俞铭.一种基于ExtJs的风玫瑰图绘制方法[J].江苏科技信息，2020，37（17）：46-48.

[22] 段炼，岳炼，张展硕，等.基于python的WRF模式后处理研究[J].信息技术与信息化，2020（9）：50-53.

[23] 申伟强，马欣，李剑.C#与Matlab混合编程及其在气象数据可视化中的应用[J].科技创新导报，2013，10（3）：56-57.

[24] 周积强，魏巍，姚肖萌，等.基于Python的三维风场风玫瑰图绘制方法研究与应用[J].气象水文海洋仪器，2022，39（3）：83-86.

[25] Hu Weiyang.Importance of regional PM2.5 transport and precipitation washout in heavy air pollution in the Twain-Hu Basin over Central China：Observational analysis and WRF-Chem simulation[J].Science of the Total Environment，2021，758：143710.

[26] Mazzeo A，Zhong J，Hood C，et al.Modelling the impact of national vs.local emission reduction on PM2.5 in the west Midlands，UK using WRF-CMAQ[J].Atmosphere，2022，13（3）：377.

[27] Wang Q，Luo K，Fan J R，et al.Spatial distribution and multiscale transport characteristics of PM2.5 in China[J].Aerosol and Air Quality Research，2019，19（9）：1993-2007.

[28] Zhang Y L，Zhu B，Gao J H，et al.The source apportionment of primary PM2.5 in an aerosol pollution event over Beijing-Tianjin-Hebei region using WRF-chem，China[J].Aerosol and Air Quality Research，2017，17（12）：2966-2980.

[29] 阙金煌.基于Anaconda环境下的Python数据分析及可视化[J].信息技术与信息化，2021（4）：215-218.

【通联编辑：谢媛媛】