Python 氣象資料分析 -- 《Python 資料分析實戰》__Python

選取了10個城市。隨後將分析它們的天氣資料，其中5個城市在距海100公裡範圍內，其餘5個距海100～400公裡。

選作樣本的城市列表如下：
Ferrara（費拉拉）
Torino（都靈）
Mantova（曼托瓦）
Milano（米蘭）
Ravenna（拉文納）
Asti（阿斯蒂）
Bologna（博洛尼亞）
Piacenza（皮亞琴察）
Cesena（切塞納）
Faenza（法恩莎）

資料來源：http://openweathermap.org/

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

import matplotlib.pyplot as pltimport matplotlib.dates as mdatesfrom dateutil import parserimport pandas as pdimport numpy as npdf_ferrara = pd.read_csv('ferrara_270615.csv')df_milano = pd.read_csv('milano_270615.csv')df_mantova = pd.read_csv('mantova_270615.csv')df_ravenna = pd.read_csv('ravenna_270615.csv')df_torino = pd.read_csv('torino_270615.csv')df_asti = pd.read_csv('asti_270615.csv')df_bologna = pd.read_csv('bologna_270615.csv')df_piacenza = pd.read_csv('piacenza_270615.csv')df_cesena = pd.read_csv('cesena_270615.csv')df_faenza = pd.read_csv('faenza_270615.csv')# 讀取城市氣象資料# 取出要分析的溫度和日期資料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']# 把日期資料轉換成 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]dat_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]# 調用 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)#這裡需要畫出三根線，所以需要三組參數ax.plot(day_ravenna,y1,'r',day_faenza,y2,'r',day_cesena,y3,'r')ax.plot(dat_milano,y4,'g',day_asti,y5,'g',day_torino,y6,'g')#顯示映像fig

結果：

離海最近的三個城市的最高氣溫比離海最遠的三個城市低不少，而最低氣溫看起來差別較小。

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

#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_min:存放每個城市最低溫度的列表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 = [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')fig#scikit-learn庫的SVR方法from sklearn.svm import SVR# dist1是靠近海的城市集合，dist2是遠離海洋的城市集合dist1 = dist[0:5]dist2 = dist[5:10]# 改變列表的結構，dist1現在是5個列表的集合# 之後我們會看到 numpy 中 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)# 加入資料，進行擬合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 軸的取值範圍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')figprint svr_lin1.coef_  #斜率print svr_lin1.intercept_  # 截距print svr_lin2.coef_print svr_lin2.intercept_

結果：

離海60公裡以內，氣溫上升速度很快，從28度陡升至31度，隨後增速漸趨緩和（如果還繼續增長的話），更長的距離才會有小幅上升。這兩種趨勢可分別用兩條直線來表示，直線的運算式為：x = ax + b
其中a為斜率，b為截距。
考慮將這兩條直線的交點作為受海洋影響和不受海洋影響的地區的分界點，或者至少是海洋影響較弱的分界點。

#考慮將這兩條直線的交點作為受海洋影響和不受海洋影響的地區的分界點，或者至少是海洋影響較弱的分界點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))fig, ax = plt.subplots()x=np.linspace(0,300,31)ax.plot(x,line1(x),x,line2(x),result,line1(result),'ro')fig

結果：

執行上述代碼，將得到交點的座標[x,y] = [ 53, 30 ]
因此，可以說海洋對氣溫產生影響的平均距離（該天的情況）為53公裡。

現在，分析最低氣溫。

#最低溫fig, ax = plt.subplots()plt.axis((0,400,15,25))ax.plot(dist,temp_min,'bo')fig

結果：

很明顯夜間或早上6點左右的最低溫與海洋無關

2.濕度資料分析
可以考察當天三個近海城市和三個內陸城市的濕度趨勢。

import matplotlib.pyplot as pltimport matplotlib.dates as mdatesfrom dateutil import parserimport pandas as pdimport numpy as npdf_ferrara = pd.read_csv('ferrara_270615.csv')df_milano = pd.read_csv('milano_270615.csv')df_mantova = pd.read_csv('mantova_270615.csv')df_ravenna = pd.read_csv('ravenna_270615.csv')df_torino = pd.read_csv('torino_270615.csv')df_asti = pd.read_csv('asti_270615.csv')df_bologna = pd.read_csv('bologna_270615.csv')df_piacenza = pd.read_csv('piacenza_270615.csv')df_cesena = pd.read_csv('cesena_270615.csv')df_faenza = pd.read_csv('faenza_270615.csv')# 讀取城市濕度資料# 取出要分析的濕度和日期資料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']# 把日期資料轉換成 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]dat_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]# 調用 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)#這裡需要畫出三根線，所以需要三組參數ax.plot(day_ravenna,y1,'r',day_faenza,y2,'r',day_cesena,y3,'r')ax.plot(dat_milano,y4,'g',day_asti,y5,'g',day_torino,y6,'g')#顯示映像fig

結果：

乍看上去好像近海城市的濕度要大於內陸城市，全天濕度差距在20%左右。再來看一下濕度的極值和離海遠近之間的關係

#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]        ]# 擷取最大濕度資料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()]fig, ax = plt.subplots()plt.plot(dist,hum_max,'bo')# 擷取最小濕度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()]#fig, ax = plt.subplots()plt.plot(dist,hum_min,'ro')

結果：

近海城市無論是最大還是最小濕度都要高於內陸城市。然而，還不能說濕度和距離之間存線上性關係或者其他能用曲線表示的關係。採集的資料點數量（10）太少，不足以描述這類趨勢。

3.風向頻率玫瑰圖
在採集的每個城市的氣象資料中，下面兩個與風有關：
風力（風向） 風速

分析存放每個城市氣象資料的DataFrame就會發現，風速不僅跟一天的時間段相關聯，還與一個介於0~360度的方向有關。例如，每一條測量資料也包含風吹來的方向
為了更好地分析這類資料，有必要將其做成可視化形式，但是對於風力資料，將其製作成使用笛卡兒座標系的線性圖不再是最佳選擇。

要是把一個DataFrame中的資料點做成散佈圖

#散佈圖不直觀，給個例子fig, ax = plt.subplots()plt.plot(df_ravenna['wind_deg'],df_ravenna['wind_speed'],'ro')#figplt.show()

結果：

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

首先，建立一個長條圖，也就是將360度分為八個面元，每個面元為45度，把所有的資料點分到這八個面元中。

#表示呈360度分布的資料點，使用另一種可視化方法：極區圖#360°八等份，分為八個面元，每份45°，把所有的資料點分到這八個面元中#histogram()函數返回結果中的數組hist為落在每個面元的資料點數量。[ 0 5 11 1 0 1 0 0]#返回結果中的數組bins定義了360度範圍內各面元的邊界。[ 0. 45. 90. 135. 180. 225. 270. 315. 360.]hist,bins=np.histogram(df_ravenna['wind_deg'],8,[0,360])print histprint bins

histogram()函數返回結果中的數組hist為落在每個面元的資料點數量。

[ 0 5 11 1 0 1 0 0]

返回結果中的數組bins定義了360度範圍內各面元的邊界。

[ 0. 45. 90. 135. 180. 225. 270. 315. 360.]

要想正確定義極區圖，離不開這兩個數組。建立一個函數來繪製極區圖，把這個函數定義為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(0.,2 * np.pi, 2 * np.pi / N)    radii = np.array(values)    fig,ax=plt.subplots()    # 繪製極區圖的座標系    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)    fig

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

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

結果：

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

從剛得到的極區圖可以得知風向在極座標系中的分布方式。該圖表示這一天大部分時間風都

吹向西南和正西方向。

定義好showRoseWind()函數之後，查看其他城市的風向情況也非常簡單。

#其他城市hist, bin = np.histogram(df_ferrara['wind_deg'],8,[0,360])print histshowRoseWind(hist,'Ferrara', max(hist))hist, bin = np.histogram(df_milano['wind_deg'],8,[0,360])print histshowRoseWind(hist,'milano', max(hist))

結果：


計算風速均值的分布情況

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

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

#計算風速均值的分布情況def RoseWind_Speed(df_city):    degs=np.arange(45,361,45)    tmp=[]    for deg in degs:        #擷取wind_deg在指定範圍的平均風速        #擷取的是風向大於 'deg-46' 度和風向小於 'deg' 的資料。        tmp.append(df_city[(df_city['wind_deg']>(deg-46)) & (df_city['wind_deg']<deg)]        ['wind_speed'].mean())    return np.array[tmp]

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

hist,bins=np.histogram(df_ravenna['wind_deg'],8,[0,360])print histprint binsshowRoseWind(RoseWind_Speed(df_ravenna),'Ravenna',max(hist))