大氣污染物濃度風域空間分布圖的繪制方法

楊洛祺 顧世杰 劉曉帥 胡引翠 武帥 田冰

摘要：風玫瑰圖主要用以表達風向頻率或各方向上的平均風速。若要可視化表達各風向上不同風速情境下的大氣污染物濃度的分異特征，需要突破當前風玫瑰圖僅展示風速風向信息的限制。為解決這一問題，文章在風玫瑰圖基礎上，采用暈渲法，用深淺不同的色調表示研究區域在統計時段內大氣污染物濃度在各風向、不同風速下的特征。基于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.

【通聯編輯：謝媛媛】