Python——氣象數據分析

將對意大利北部沿海地區的氣象數據進行分析與可視化。我們在實驗過程中先會運用 Python 中 matplotlib 庫的對數據進行圖表化處理，然后調用 scikit-learn 庫當中的的 SVM 庫對數據進行回歸分析，最終在圖表分析的支持下得出我們的結論。

筆記來源

圖靈教育 的 《Python 數據分析實戰》 第 2 章

知識點

• matplotlib 庫畫出圖像
• scikit-learn 庫對數據進行回歸分析
• numpy 庫對數據進行切片

原理

氣象數據是在網上很容易找到的一類數據。很多網站都提供以往的氣壓、氣溫、濕度和降雨量等氣象數據。只需指定位置和日期，就能獲取一個氣象數據文件。這些測量數據是由氣象站收集的。氣象數據這類數據源涵蓋的信息范圍較廣。數據分析的目的是把原始數據轉化為信息，再把信息轉化為知識，因此拿氣象數據作為數據分析的對象來講解數據分析全過程再合適不過。

待檢驗的假設：靠海對氣候的影響

寫作本章時，雖正值夏初，卻已酷熱難耐，住在大城市的人感受更為強烈。於是周末很多人到山村或海濱城市去游玩，放松一下身心，遠離內陸城市的悶熱天氣。我常常想，靠海對氣候有什么影響？這個問題可以作為數據分析的一個不錯的出發點。我不想把本章寫成科學類讀物，只是想借助這樣一種方式，讓數據分析愛好者能夠把所學用於實踐，解決 “海洋對一個地區的氣候有何影響” 這個問題。

研究系統：亞得里亞海和波河流域

既然已定義好問題，就需要尋找適合研究數據的系統，提供適合回答這個問題的環境。首先，需要找到一片海域供你研究。我住在意大利，可選擇的海有很多，因為意大利是一個被海洋包圍的半島國家。為什么要把自己的選擇局限在意大利呢？因為我們所研究的問題剛好和意大利人的一種典行為相關，也就是夏天我們喜歡躲在海邊，以躲避內陸的酷熱。我不知道在其他國家這種行為是否也很普遍，因此我只把自己熟悉的意大利作為一個系統進行研究。但是你可能會考慮研究意大利的哪個地區呢？上面說過，意大利是半島國家，找到可研究的海域不是問題，但是如何衡量海洋對其遠近不同的地方的影響？這就引出了一個大問題。意大利其實多山地，離海差不多遠，可以彼此作為參照的內陸區域較少。為了衡量海洋對氣候的影響，我排除了山地，因為山地也許會引入其他很多因素，比如海拔。

意大利波河流域這塊區域就很適合研究海洋對氣候的影響。這一片平原東起亞得里亞海，向內陸延伸數百公里（見圖 9-1）。它周邊雖不乏群山環繞，但由於它很寬廣，削弱了群山的影響。此外，該區域城鎮密集，也便於選取一組離海遠近不同的城市。我們所選的幾個城市，兩個城市間的最大距離約為 400 公里。

第一步，選 10 個城市作為參照組。選擇城市時，注意它們要能代表整個平原地區（見圖 9-2）。

如圖 9-2 所示，我們選取了 10 個城市。隨后將分析它們的天氣數據，其中 5 個城市在距海 100 公里范圍內，其余 5 個距海 100～400 公里。

選作樣本的城市列表如下：

• Ferrara（費拉拉）
• Torino（都靈）
• Mantova（曼托瓦）
• Milano（米蘭）
• Ravenna（拉文納）
• Asti（阿斯蒂）
• Bologna（博洛尼亞）
• Piacenza（皮亞琴察）
• Cesena（切塞納）
• Faenza（法恩莎）

現在，我們需要弄清楚這些城市離海有多遠。方法有多種。這里使用 TheTimeNow 網站提供的服務，它支持多種語言（見圖 9-3）。

有了計算兩城市間距離這樣的服務，我們就可以計算每個城市與海之間的距離。

你可以選擇海濱城市 Comacchio 作為基點，計算其他城市與它之間的距離（見圖 9-2）。使用上述服務計算完所有距離后，得到的結果如表 9-1 所示。

開發准備

定義好要研究的系統之后，我們就需要創建數據源，以獲取研究所需的數據。上網瀏覽一番，你就會發現很多網站都提供世界各地的氣象數據，其中就有 Open Weather Map，它的網址是 http://openweathermap.org/ （見圖 9-4）。

該網站提供以下功能：在請求的 URL 中指定城市，即可獲取該城市的氣象數據。我們已經准備好了數據，不需要大家再去調用該網站的 API。 下面，就先下載氣象數據。

網盤鏈接：https://pan.baidu.com/s/1AJRiv--ysy265_it_7s3Qg 提取碼：jyue

這時候，我們通過 tree 命令應該能夠再 WeatherData 中間看到 10 個城市的天氣數據文件（以 .csv 結尾）

!apt-get install tree

!tree WeatherData/

導入相關包開始實驗。

import numpy as np
import pandas as pd
import datetime

如果你想用本章的數據，需要加載寫作本章時保存的 10 個 CSV 文件。

df_ferrara = pd.read_csv('WeatherData/ferrara_270615.csv')
df_milano = pd.read_csv('WeatherData/milano_270615.csv')
df_mantova = pd.read_csv('WeatherData/mantova_270615.csv')
df_ravenna = pd.read_csv('WeatherData/ravenna_270615.csv')
df_torino = pd.read_csv('WeatherData/torino_270615.csv')
df_asti = pd.read_csv('WeatherData/asti_270615.csv')
df_bologna = pd.read_csv('WeatherData/bologna_270615.csv')
df_piacenza = pd.read_csv('WeatherData/piacenza_270615.csv')
df_cesena = pd.read_csv('WeatherData/cesena_270615.csv')
df_faenza = pd.read_csv('WeatherData/faenza_270615.csv')

我們把這些數據讀入內存，完成了實驗准備的部分。

實驗步驟

從數據可視化入手分析收集到的數據是常見的做法。前面講過，matplotlib 庫提供一系列圖表生成工具，能夠以可視化形式表示數據。數據可視化在數據分析階段非常有助於發現研究系統的一些特點。

導入以下必要的庫：

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from dateutil import parser

溫度數據分析

舉例來說，非常簡單的分析方法是先分析一天中氣溫的變化趨勢。我們以城市米蘭為例。

# 取出我們要分析的溫度和日期數據
y1 = df_milano['temp']
x1 = df_milano['day']

# 把日期數據轉換成 datetime 的格式
day_milano = [parser.parse(x) for x in x1]

# 調用 subplot 函數, fig 是圖像對象，ax 是坐標軸對象
fig, ax = plt.subplots()

# 調整x軸坐標刻度，使其旋轉70度，方便查看
plt.xticks(rotation=70)

# 設定時間的格式
hours = mdates.DateFormatter('%H:%M')

# 設定X軸顯示的格式
ax.xaxis.set_major_formatter(hours)

# 畫出圖像，day_milano是X軸數據，y1是Y軸數據，‘r’代表的是'red' 紅色
ax.plot(day_milano ,y1, 'r')

執行上述代碼，將得到如圖 9-8 所示的圖像。由圖可見，氣溫走勢接近正弦曲線，從早上開始氣溫逐漸升高，最高溫出現在下午兩點到六點之間，隨后氣溫逐漸下降，在第二天早上六點時達到最低值。

我們進行數據分析的目的是嘗試解釋是否能夠評估海洋是怎樣影響氣溫的，以及是否能夠影響氣溫趨勢，因此我們同時來看幾個不同城市的氣溫趨勢。這是檢驗分析方向是否正確的唯一方式。因此，我們選擇三個離海最近以及三個離海最遠的城市。

# 讀取溫度和日期數據
y1 = df_ravenna['temp']
x1 = df_ravenna['day']
y2 = df_faenza['temp']
x2 = df_faenza['day']
y3 = df_cesena['temp']
x3 = df_cesena['day']
y4 = df_milano['temp']
x4 = df_milano['day']
y5 = df_asti['temp']
x5 = df_asti['day']
y6 = df_torino['temp']
x6 = df_torino['day']

# 把日期從 string 類型轉化為標准的 datetime 類型
day_ravenna = [parser.parse(x) for x in x1]
day_faenza = [parser.parse(x) for x in x2]
day_cesena = [parser.parse(x) for x in x3]
day_milano = [parser.parse(x) for x in x4]
day_asti = [parser.parse(x) for x in x5]
day_torino = [parser.parse(x) for x in x6]

# 調用 subplots() 函數，重新定義 fig, ax 變量
fig, ax = plt.subplots()
plt.xticks(rotation=70)

hours = mdates.DateFormatter('%H:%M')
ax.xaxis.set_major_formatter(hours)

#這里需要畫出三根線，所以需要三組參數， 'g'代表'green'
ax.plot(day_ravenna,y1,'r',day_faenza,y2,'r',day_cesena,y3,'r')
ax.plot(day_milano,y4,'g',day_asti,y5,'g',day_torino,y6,'g')

上述代碼將生成如圖 9-9 所示的圖表。離海最近的三個城市的氣溫曲線使用紅色，而離海最遠的三個城市的曲線使用綠色。

如圖 9-9 所示，結果看起來不錯。離海最近的三個城市的最高氣溫比離海最遠的三個城市低不少，而最低氣溫看起來差別較小。

我們可以沿着這個方向做深入研究，收集 10 個城市的最高溫和最低溫，用線性圖表示氣溫最值點和離海遠近之間的關系。

# dist 是一個裝城市距離海邊距離的列表
dist = [df_ravenna['dist'][0],
    df_cesena['dist'][0],
    df_faenza['dist'][0],
    df_ferrara['dist'][0],
    df_bologna['dist'][0],
    df_mantova['dist'][0],
    df_piacenza['dist'][0],
    df_milano['dist'][0],
    df_asti['dist'][0],
    df_torino['dist'][0]
]

# temp_max 是一個存放每個城市最高溫度的列表
temp_max = [df_ravenna['temp'].max(),
    df_cesena['temp'].max(),
    df_faenza['temp'].max(),
    df_ferrara['temp'].max(),
    df_bologna['temp'].max(),
    df_mantova['temp'].max(),
    df_piacenza['temp'].max(),
    df_milano['temp'].max(),
    df_asti['temp'].max(),
    df_torino['temp'].max()
]

# temp_min 是一個存放每個城市最低溫度的列表
temp_min = [df_ravenna['temp'].min(),
    df_cesena['temp'].min(),
    df_faenza['temp'].min(),
    df_ferrara['temp'].min(),
    df_bologna['temp'].min(),
    df_mantova['temp'].min(),
    df_piacenza['temp'].min(),
    df_milano['temp'].min(),
    df_asti['temp'].min(),
    df_torino['temp'].min()
]

先把最高溫畫出來。

fig, ax = plt.subplots()
ax.plot(dist,temp_max,'ro')

結果如圖 9-10 所示。

如圖 9-10 所示，現在你可以證實，海洋對氣象數據具有一定程度的影響這個假設是正確的（至少這一天如此）。進一步觀察上圖，你會發現海洋的影響衰減得很快，離海 60～70 公里開外，氣溫就已攀升到高位。

用線性回歸算法得到兩條直線，分別表示兩種不同的氣溫趨勢，這樣做很有趣。我們可以使用 scikit-learn 庫的 SVR 方法。

！注意：這段代碼會跑比較久的時間，請耐心等待

from sklearn.svm import SVR

# dist1是靠近海的城市集合，dist2是遠離海洋的城市集合
dist1 = dist[0:5]
dist2 = dist[5:10]

# 改變列表的結構，dist1現在是5個列表的集合
# 之后我們會看到 nbumpy 中 reshape() 函數也有同樣的作用
dist1 = [[x] for x in dist1]
dist2 = [[x] for x in dist2]

# temp_max1 是 dist1 中城市的對應最高溫度
temp_max1 = temp_max[0:5]
# temp_max2 是 dist2 中城市的對應最高溫度
temp_max2 = temp_max[5:10]

# 我們調用SVR函數，在參數中規定了使用線性的擬合函數
# 並且把 C 設為1000來盡量擬合數據（因為不需要精確預測不用擔心過擬合）
svr_lin1 = SVR(kernel='linear', C=1e3)
svr_lin2 = SVR(kernel='linear', C=1e3)

# 加入數據，進行擬合（這一步可能會跑很久，大概10多分鍾，休息一下:) ）
svr_lin1.fit(dist1, temp_max1)
svr_lin2.fit(dist2, temp_max2)

# 關於 reshape 函數請看代碼后面的詳細討論
xp1 = np.arange(10,100,10).reshape((9,1))
xp2 = np.arange(50,400,50).reshape((7,1))
yp1 = svr_lin1.predict(xp1)
yp2 = svr_lin2.predict(xp2)

然后繪圖

# 限制了 x 軸的取值范圍
fig, ax = plt.subplots()
ax.set_xlim(0,400)

# 畫出圖像
ax.plot(xp1, yp1, c='b', label='Strong sea effect')
ax.plot(xp2, yp2, c='g', label='Light sea effect')
ax.plot(dist,temp_max,'ro')

這里 np.arange(10,100,10) 會返回 [10, 20, 30,..., 90]，如果把列表看成是一個矩陣，那么這個矩陣是 1 9 的。這里 reshape((9,1)) 函數就會把該列表變為 9 1 的， [[10], [20], ..., [90]]。這么做的原因是因為 predict() 函數的只能接受一個 N 1 的列表，返回一個 1 N 的列表。

上述代碼將生成如圖 9-11 所示的圖像。

如上所見，離海 60 公里以內，氣溫上升速度很快，從 28 度陡升至 31 度，隨后增速漸趨緩和（如果還繼續增長的話），更長的距離才會有小幅上升。這兩種趨勢可分別用兩條直線來表示，直線的表達式為：y=ax+by=ax+b 其中 a 為斜率，b 為截距。

print(svr_lin1.coef_)  #斜率
print(svr_lin1.intercept_)  # 截距
print(svr_lin2.coef_)
print(svr_lin2.intercept_)

你可能會考慮將這兩條直線的交點作為受海洋影響和不受海洋影響的區域的分界點，或者至少是海洋影響較弱的分界點。

from scipy.optimize import fsolve

# 定義了第一條擬合直線
def line1(x):
    a1 = svr_lin1.coef_[0][0]
    b1 = svr_lin1.intercept_[0]
    return a1*x + b1

# 定義了第二條擬合直線
def line2(x):
    a2 = svr_lin2.coef_[0][0]
    b2 = svr_lin2.intercept_[0]
    return a2*x + b2

# 定義了找到兩條直線的交點的 x 坐標的函數
def findIntersection(fun1,fun2,x0):
    return fsolve(lambda x : fun1(x) - fun2(x),x0)

result = findIntersection(line1,line2,0.0)
print("[x,y] = [ %d , %d ]" % (result,line1(result)))

# x = [0,10,20, ..., 300]
x = np.linspace(0,300,31)
plt.plot(x,line1(x),x,line2(x),result,line1(result),'ro')

執行上述代碼，將得到交點的坐標 [x,y]=[53,30][x,y]=[53,30] 並得到如圖 9-12 所示的圖表。

因此，你可以說海洋對氣溫產生影響的平均距離（該天的情況）為 53 公里。現在，我們可以轉而分析最低氣溫。

# axis 函數規定了 x 軸和 y 軸的取值范圍
plt.axis((0,400,15,25))
plt.plot(dist,temp_min,'bo')

在這個例子中，很明顯夜間或早上 6 點左右的最低溫與海洋無關。如果沒記錯的話，小時候老師教給大家的是海洋能夠緩和低溫，或者說夜間海洋釋放白天吸收的熱量。但是從我們得到情況來看並非如此。我們剛使用的是意大利夏天的氣溫數據，而驗證該假設在冬天或其他地方是否也成立，將會非常有趣。

濕度數據分析

10 個 DataFrame 對象中還包含濕度這個氣象數據。因此，你也可以考察當天三個近海城市和三個內陸城市的濕度趨勢。

# 讀取濕度數據
y1 = df_ravenna['humidity']
x1 = df_ravenna['day']
y2 = df_faenza['humidity']
x2 = df_faenza['day']
y3 = df_cesena['humidity']
x3 = df_cesena['day']
y4 = df_milano['humidity']
x4 = df_milano['day']
y5 = df_asti['humidity']
x5 = df_asti['day']
y6 = df_torino['humidity']
x6 = df_torino['day']

# 重新定義 fig 和 ax 變量
fig, ax = plt.subplots()
plt.xticks(rotation=70)

# 把時間從 string 類型轉化為標准的 datetime 類型
day_ravenna = [parser.parse(x) for x in x1]
day_faenza = [parser.parse(x) for x in x2]
day_cesena = [parser.parse(x) for x in x3]
day_milano = [parser.parse(x) for x in x4]
day_asti = [parser.parse(x) for x in x5]
day_torino = [parser.parse(x) for x in x6]

# 規定時間的表示方式
hours = mdates.DateFormatter('%H:%M')
ax.xaxis.set_major_formatter(hours)

#表示在圖上
ax.plot(day_ravenna,y1,'r',day_faenza,y2,'r',day_cesena,y3,'r')
ax.plot(day_milano,y4,'g',day_asti,y5,'g',day_torino,y6,'g')

上述代碼將生成如圖 9-14 所示的圖表。

乍看上去好像近海城市的濕度要大於內陸城市，全天濕度差距在 20% 左右。我們再來看一下濕度的極值和離海遠近之間的關系，是否跟我們的第一印象相符。

# 獲取最大濕度數據
hum_max = [df_ravenna['humidity'].max(),
df_cesena['humidity'].max(),
df_faenza['humidity'].max(),
df_ferrara['humidity'].max(),
df_bologna['humidity'].max(),
df_mantova['humidity'].max(),
df_piacenza['humidity'].max(),
df_milano['humidity'].max(),
df_asti['humidity'].max(),
df_torino['humidity'].max()
]

plt.plot(dist,hum_max,'bo')

我們把 10 個城市的最大濕度與離海遠近之間的關系做成圖表，請見圖 9-15。

# 獲取最小濕度
hum_min = [
df_ravenna['humidity'].min(),
df_cesena['humidity'].min(),
df_faenza['humidity'].min(),
df_ferrara['humidity'].min(),
df_bologna['humidity'].min(),
df_mantova['humidity'].min(),
df_piacenza['humidity'].min(),
df_milano['humidity'].min(),
df_asti['humidity'].min(),
df_torino['humidity'].min()
]
plt.plot(dist,hum_min,'bo')

再來把 10 個城市的最小濕度與離海遠近之間的關系做成圖表，請見圖 9-16。

由圖 9-15 和圖 9-16 可以確定，近海城市無論是最大還是最小濕度都要高於內陸城市。然而,在我看來，我們還不能說濕度和距離之間存在線性關系或者其他能用曲線表示的關系。我們采集的數據點數量（10）太少，不足以描述這類趨勢。

風向頻率玫瑰圖

在我們采集的每個城市的氣象數據中，下面兩個與風有關：

• 風力（風向）
• 風速

分析存放每個城市氣象數據的 DataFrame 就會發現，風速不僅跟一天的時間段相關聯，還與一個介於 0~360 度的方向有關。

例如，每一條測量數據也包含風吹來的方向（圖 9-17）。

為了更好地分析這類數據，有必要將其做成可視化形式，但是對於風力數據，將其制作成使用笛卡兒坐標系的線性圖不再是最佳選擇。 要是把一個 DataFrame 中的數據點做成散點圖

plt.plot(df_ravenna['wind_deg'],df_ravenna['wind_speed'],'ro')

就會得到圖 9-18 這樣的圖表，很顯然該圖的表現力也有不足。

要表示呈 360 度分布的數據點，最好使用另一種可視化方法：極區圖。

首先，創建一個直方圖，也就是將 360 度分為八個面元，每個面元為 45 度，把所有的數據點分到這八個面元中。

hist, bins = np.histogram(df_ravenna['wind_deg'],8,[0,360])
print(hist)
print(bins)

histogram() 函數返回結果中的數組 hist 為落在每個面元的數據點數量。[0 5 11 1 0 1 0 0]

返回結果中的數組 bins 定義了 360 度范圍內各面元的邊界。[0. 45. 90. 135. 180. 225. 270. 315. 360.]

要想正確定義極區圖，離不開這兩個數組。我們將創建一個函數來繪制極區圖，其中部分代碼在第 7 章已講過。我們把這個函數定義為 showRoseWind()，它有三個參數：values 數組，指的是想為其作圖的數據，也就是這里的 hist 數組；第二個參數 city_name 為字符串類型，指定圖表標題所用的城市名稱；最后一個參數 max_value 為整型，指定最大的藍色值。

定義這樣一個函數很有用，它既能避免多次重復編寫相同的代碼，還能增強代碼的模塊化程度，便於你把精力放到與函數內部操作相關的概念上。

def showRoseWind(values,city_name,max_value):
    N = 8

    # theta = [pi*1/4, pi*2/4, pi*3/4, ..., pi*2]
    theta = np.arange(2 * np.pi / 16, 2 * np.pi, 2 * np.pi / 8)
    radii = np.array(values)
    # 繪制極區圖的坐標系
    plt.axes([0.025, 0.025, 0.95, 0.95], polar=True)

    # 列表中包含的是每一個扇區的 rgb 值，x越大，對應的color越接近藍色
    colors = [(1-x/max_value, 1-x/max_value, 0.75) for x in radii]

    # 畫出每個扇區
    plt.bar(theta, radii, width=(2*np.pi/N), bottom=0.0, color=colors)

    # 設置極區圖的標題
    plt.title(city_name, x=0.2, fontsize=20)

你需要修改變量 colors 存儲的顏色表。這里，扇形的顏色越接近藍色，值越大。定義好函數之后，調用它即可：

showRoseWind(hist,'Ravenna',max(hist))

運行上述函數，將得到如圖 9-19 所示的極區圖。

由圖 9-19 可見，整個 360 度的范圍被分成八個區域（面元），每個區域弧長為 45 度，此外每個區域還有一列呈放射狀排列的刻度值。在每個區域中，用半徑長度可以改變的扇形表示一個數值，半徑越長，扇形所表示的數值就越大。為了增強圖表的可讀性，我們使用與扇形半徑相對應的顏色表。半徑越長，扇形跨度越大，顏色越接近於深藍色。

從剛得到的極區圖可以得知風向在極坐標系中的分布方式。該圖表示這一天大部分時間風都吹向西南和正西方向。定義好 showRoseWind() 函數之后，查看其他城市的風向情況也非常簡單。

hist, bin = np.histogram(df_ferrara['wind_deg'],8,[0,360])
print(hist)
showRoseWind(hist,'Ferrara', max(hist))

計算風速均值的分布情況

即使是跟風速相關的其他數據，也可以用極區圖來表示。

定義 RoseWind_Speed 函數，計算將 360 度范圍划分成的八個面元中每個面元的平均風速。

def RoseWind_Speed(df_city):
    # degs = [45, 90, ..., 360]
    degs = np.arange(45,361,45)
    tmp = []
    for deg in degs:
        # 獲取 wind_deg 在指定范圍的風速平均值數據
        tmp.append(df_city[(df_city['wind_deg']>(deg-46)) & (df_city['wind_deg']<deg)]
        ['wind_speed'].mean())
    return np.array(tmp)

這里 df_city[(df_city['wind_deg']>(deg-46)) & (df_city['wind_deg'] 獲取的是風向大於 deg-46度和風向小於deg` 的數據。

RoseWind_Speed() 函數返回一個包含八個平均風速值的 NumPy 數組。該數組將作為先前定義的 showRoseWind() 函數的第一個參數，這個函數是用來繪制極區圖的。

showRoseWind(RoseWind_Speed(df_ravenna),'Ravenna',max(hist))

圖 9-21 所示的風向頻率玫瑰圖表示風速在 360 度范圍內的分布情況。

實驗總結

本章主要目的是演示如何從原始數據獲取信息。其中有些信息無法給出重要結論，而有些信息能夠驗證假設，增加我們對系統狀態的認識，而找出這種信息也就意味着數據分析取得了成功。