本實(shí)驗(yàn)將對意大利北部沿海地區(qū)的氣象數(shù)據(jù)進(jìn)行分析與可視化鬓催。我們在實(shí)驗(yàn)過程中先會運(yùn)用 Python 中 matplotlib 庫的對數(shù)據(jù)進(jìn)行圖表化處理改淑，然后調(diào)用 scikit-learn 庫當(dāng)中的的 SVM 庫對數(shù)據(jù)進(jìn)行回歸分析突勇，最終在圖表分析的支持下得出我們的結(jié)論年栓。
下載氣象數(shù)據(jù)
!wget -nc http://labfile.oss.aliyuncs.com/courses/780/WeatherData.zip
裝 unzip 解壓縮
!apt-get install unzip
解壓縮
!unzip -o WeatherData.zip
這時候惹盼，我們通過 tree 命令應(yīng)該能夠再 WeatherData 中間看到 10 個城市的天氣數(shù)據(jù)文件（以 .csv 結(jié)尾）
!apt-get install tree
!tree WeatherData/
導(dǎo)入相關(guān)包開始實(shí)驗(yàn)
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')
我們把這些數(shù)據(jù)讀入內(nèi)存春瞬，完成了實(shí)驗(yàn)準(zhǔn)備的部分
實(shí)驗(yàn)步驟
導(dǎo)入以下必要的庫：
matplotlib 庫提供一系列圖表生成工具，能夠以可視化形式表示數(shù)據(jù)
%matplotlib inline
import matplotlib.pylot as plt
import matplotlib.dates as mdates
from dateutil import parser
#取出milano的溫度和日期數(shù)據(jù)
y1 = df_milano['temp']
x1 = df_milano['day']
#將日期數(shù)據(jù)改為datetime的格式
day_milano = [parser.parse(x) for x in x1]
#調(diào)用subplot函數(shù)鸭津，fig是圖像對象彤侍，ax是坐標(biāo)軸對象
fig,ax = plt.subplots()
#調(diào)整x軸坐標(biāo)刻度，使其旋轉(zhuǎn)70度逆趋，方便查看
plt.xticks(rotation=70)
#設(shè)定時間的格式
hours=mdates.DateFormatter('%H:%M')
#設(shè)定X軸顯示的格式
ax.xaxis.set_major_formatter(hours)
#畫出圖像盏阶，day_milano是X軸數(shù)據(jù)，temp是Y軸數(shù)據(jù)闻书，‘r’代表紅色
ax.plot(day_milano,y1,'r')
由圖可見名斟，氣溫走勢接近正弦曲線，從早上開始?xì)鉁刂饾u升高魄眉，最高溫出現(xiàn)在下午兩點(diǎn)到六點(diǎn)之間砰盐，隨后氣溫逐漸下降，在第二天早上六點(diǎn)時達(dá)到最低值坑律。
我們進(jìn)行數(shù)據(jù)分析的目的是嘗試解釋是否能夠評估海洋是怎樣影響氣溫的岩梳，以及是否能夠影響氣溫趨勢，因此我們同時來看幾個不同城市的氣溫趨勢。這是檢驗(yàn)分析方向是否正確的唯一方式冀值。因此也物，我們選擇三個離海最近以及三個離海最遠(yuǎn)的城市。
#選擇三個離海最近以及三個離海最遠(yuǎn)的城市列疗。
#用循環(huán)會不會更方便……滑蚯？
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]
#和上面例子一樣
fig,ax = plt.subplots()
plt.xticks(rotation=70)
hours=mdates.DateFormatter('%H:%M')
ax.xaxis.set_major_formatter(hours)
#近海的3條，遠(yuǎn)海的3條
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')
離海最近的三個城市的氣溫曲線使用紅色抵栈，而離海最遠(yuǎn)的三個城市的曲線使用綠色告材。
離海最近的三個城市的最高氣溫比離海最遠(yuǎn)的三個城市低不少，而最低氣溫看起來差別較小古劲。
我們可以沿著這個方向做深入研究创葡，收集 10 個城市的最高溫和最低溫，用線性圖表示氣溫最值點(diǎn)和離海遠(yuǎn)近之間的關(guān)系绢慢。
# 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(),]
#把最高溫畫出來,x軸為距離。y軸為最高溫度
fig,ax=plt.subplots()
ax.plot(dist,temp_max,'ro')
進(jìn)一步觀察上圖洛波，你會發(fā)現(xiàn)海洋的影響衰減得很快胰舆，離海 60～70 公里開外，氣溫就已攀升到高位蹬挤。
用線性回歸算法得到兩條直線缚窿，分別表示兩種不同的氣溫趨勢，這樣做很有趣焰扳。我們可以使用 scikit-learn 庫的 SVR 方法倦零。
from sklearn.svm import SVR
# dist1是靠近海的城市集合，dist2是遠(yuǎn)離海洋的城市集合
dist1=dist[0:5]
dist2=dist[5:10]
#改變列表的結(jié)構(gòu)吨悍，dist1現(xiàn)在是5個列表的集合
# 之后我們會看到 nbumpy 中 reshape() 函數(shù)也有同樣的作用
dist1=[[x]for x in dist1]
dist2=[[x]for x in dist2]
#temp_max1 是 dist1 中城市的對應(yīng)最高溫度,temp_max2同理
temp_max1=temp_max[0:5]
temp_max2=temp_max[5:10]
#調(diào)用SVR函數(shù)扫茅，規(guī)定為線性擬合
#并且把 C 設(shè)為1000來盡量擬合數(shù)據(jù)（因?yàn)椴恍枰_預(yù)測不用擔(dān)心過擬合）
#C是什么？我并不清楚
svr_lin1=SVR(kernel='linear',C=1e3)
svr_lin2=SVR(kernel='linear',C=1e3)
#帶入數(shù)據(jù)進(jìn)行擬合
svr_lin1.fit(dist1,temp_max1)
svr_lin2.fit(dist2,temp_max2)
#關(guān)于 reshape 函數(shù)請看代碼后面的詳細(xì)討論
xp1=np.arange(10,100,10).reshape((9,1))
xp1=np.arange(50,400,50).reshape((7,1))
yp1=svr_lin1.predict(xp1)
yp2=svr_lin2.predict(xp2)
然后畫圖
#限制x軸的取值范圍
fig,ax=plot.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)) 函數(shù)就會把該列表變?yōu)?9 1 的躏仇， [[10], [20], ..., [90]]恋脚。這么做的原因是因?yàn)?predict() 函數(shù)的只能接受一個 N 1 的列表，返回一個 1 N 的列表焰手。
如上所見糟描，離海 60 公里以內(nèi)，氣溫上升速度很快书妻，從 28 度陡升至 31 度船响，隨后增速漸趨緩和（如果還繼續(xù)增長的話），更長的距離才會有小幅上升。這兩種趨勢可分別用兩條直線來表示灿意，直線的表達(dá)式為：
查看斜率和截距：
print(svr_lin1.coef_)#斜率
print(svr_lin1.intercept_)#截距
print(svr_lin2.coef_)
print(svr_lin2.intercept_)
你可能會考慮將這兩條直線的交點(diǎn)作為受海洋影響和不受海洋影響的區(qū)域的分界點(diǎn)估灿，或者至少是海洋影響較弱的分界點(diǎn)。
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
#定義找到兩條直線交點(diǎn)坐標(biāo)x坐標(biāo)的函數(shù)
def findIntersection(fun1,fun2,x0):
return fsolve(lambda x :fun1(x)-fun2(x),x0)
result = findIntersection(line1,line2,0.0)
x=np.linspace(0,300,31)
plt.plot(x,line1(x),x,line2(x),result,line1(result),'ro')
#result,line1(result),為交點(diǎn)坐標(biāo)
print(result)
print(line1(result))
可以說海洋對氣溫產(chǎn)生影響的平均距離（該天的情況）為 53 公里
現(xiàn)在缤剧，我們可以轉(zhuǎn)而分析最低氣溫馅袁。
#規(guī)定了x軸的范圍，y軸的范圍
plt.axis((0,400,15,25))
#繪出最低溫度的圖表
plt.plot(dist,temp_min,'bo')
在這個例子中荒辕，很明顯夜間或早上 6 點(diǎn)左右的最低溫與海洋無關(guān)
濕度數(shù)據(jù)分析
10 個 DataFrame 對象中還包含濕度這個氣象數(shù)據(jù)汗销。因此，你也可以考察當(dāng)天三個近海城市和三個內(nèi)陸城市的濕度趨勢抵窒。
#讀取濕度和時間數(shù)據(jù)
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']
# 調(diào)用 subplots() 函數(shù)弛针，重新定義 fig, ax 變量
fig, ax = plt.subplots()
plt.xticks(rotation=70)
# 把日期從 string 類型轉(zhuǎn)化為標(biāo)準(zhǔn)的 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]
#規(guī)定時間表達(dá)的方式
hours = mdates.DateFormatter('%H:%M')
ax.xaxis.set_major_formatter(hours)
#這里需要畫出三根線，所以需要三組參數(shù)李皇， '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')
紅色是近海城市削茁，綠色為遠(yuǎn)海城市
乍看上去好像近海城市的濕度要大于內(nèi)陸城市，全天濕度差距在 20% 左右掉房。我們再來看一下濕度的極值和離海遠(yuǎn)近之間的關(guān)系茧跋，是否跟我們的第一印象相符。
#獲取最大濕度數(shù)據(jù)
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')
#探尋最低濕度
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')
可以確定卓囚，近海城市無論是最大還是最小濕度都要高于內(nèi)陸城市.
然而,我們還不能說濕度和距離之間存在線性關(guān)系或者其他能用曲線表示的關(guān)系瘾杭。我們采集的數(shù)據(jù)點(diǎn)數(shù)量（10）太少，不足以描述這類趨勢哪亿。
風(fēng)向頻率玫瑰圖
在我們采集的每個城市的氣象數(shù)據(jù)中粥烁，下面兩個與風(fēng)有關(guān)：
- 風(fēng)力（風(fēng)向）
- 風(fēng)速
分析存放每個城市氣象數(shù)據(jù)的 DataFrame 就會發(fā)現(xiàn)，風(fēng)速不僅跟一天的時間段相關(guān)聯(lián)蝇棉，還與一個介于 0~360 度的方向有關(guān)讨阻。
例如，每一條測量數(shù)據(jù)也包含風(fēng)吹來的方向银萍。
df_ravenna[['wind_deg','wind_speed','day']]
為了更好地分析這類數(shù)據(jù)变勇，有必要將其做成可視化形式，但是對于風(fēng)力數(shù)據(jù)贴唇，將其制作成使用笛卡兒坐標(biāo)系的線性圖不再是最佳選擇搀绣。
要是把一個 DataFrame 中的數(shù)據(jù)點(diǎn)做成散點(diǎn)圖
plt.plot(df_ravenna['wind_deg'],df_ravenna['wind_speed'],'ro')
很顯然該圖的表現(xiàn)力也有不足。
要表示呈 360 度分布的數(shù)據(jù)點(diǎn)戳气，最好使用另一種可視化方法：極區(qū)圖链患。
首先，創(chuàng)建一個直方圖瓶您，也就是將 360 度分為八個面元麻捻，每個面元為 45 度纲仍，把所有的數(shù)據(jù)點(diǎn)分到這八個面元中。
hist,bins=np.histogram(df_ravenna['wind_deg'],8,[0,360])
print(hist)
print(bins)
histogram() 函數(shù)返回結(jié)果中的數(shù)組 hist 為落在每個面元的數(shù)據(jù)點(diǎn)數(shù)量贸毕。
返回結(jié)果中的數(shù)組 bins 定義了 360 度范圍內(nèi)各面元的邊界郑叠。
#繪制玫瑰風(fēng)圖
def showRoseWind(values,city_name,max_value):
N=8
#從8/Π到2Π
theta=np.arange(2*np.pi/16,2*np.pi,2*np.pi/8)
radii=np.array(values)
# 繪制極區(qū)圖的坐標(biāo)系
plt.axes([0.025,0.025,0.95,0.95],polar=True)
# 列表中包含的是每一個扇區(qū)的 rgb 值，x越大明棍，對應(yīng)的color越接近藍(lán)色
colors = [(1-x/max_value,1-x/max_value,0.75)for x in radii]
#畫出每個扇區(qū)
plt.bar(theta,radii,width=(2*np.pi/N),bottom=0.0,color=colors)
#設(shè)置極區(qū)圖的標(biāo)題
plt.title(city_name,x=0.2,fontsize=20)
畫圖
showRoseWind(hist,'Ravenna',max(hist))
這里乡革，扇形的顏色越接近藍(lán)色，值越大.
整個 360 度的范圍被分成八個區(qū)域（面元）摊腋，每個區(qū)域弧長為 45 度沸版，此外每個區(qū)域還有一列呈放射狀排列的刻度值。
在每個區(qū)域中兴蒸，用半徑長度可以改變的扇形表示一個數(shù)值视粮，半徑越長，扇形所表示的數(shù)值就越大橙凳。
為了增強(qiáng)圖表的可讀性蕾殴，我們使用與扇形半徑相對應(yīng)的顏色表。半徑越長岛啸，扇形跨度越大区宇，顏色越接近于深藍(lán)色。
#定義好 showRoseWind() 函數(shù)之后值戳，查看其他城市的風(fēng)向情況也非常簡單。
hist,bin=np.histogram(df_ferrara['wind_deg'],8,[0,360])
print(hist)
showRoseWind(hist,'Ferrara',max(hist))
計(jì)算風(fēng)速均值的分布情況
即使是跟風(fēng)速相關(guān)的其他數(shù)據(jù)炉爆，也可以用極區(qū)圖來表示堕虹。
定義 RoseWind_Speed 函數(shù)，計(jì)算將 360 度范圍劃分成的八個面元中每個面元的平均風(fēng)速芬首。
def RoseWind_Speed(df_city):
#……赴捞？
# degs = [45, 90, ..., 360]
degs=np.arange(45,361,45)
tmp=[]
for deg in degs:
# 獲取 wind_deg 在指定范圍的風(fēng)速平均值數(shù)據(jù)
tmp.append(df_city[(df_city['wind_deg']>(deg-46))&(df_city['wind_deg']<deg)]
['wind_speed'].mean())
#這里 df_city[(df_city['wind_deg']>(deg-46)) & (df_city['wind_deg']<deg)]
#獲取的是風(fēng)向大于 deg-46 度和風(fēng)向小于 deg 的數(shù)據(jù)。
return np.array(tmp)
RoseWind_Speed() 函數(shù)返回一個包含八個平均風(fēng)速值的 NumPy 數(shù)組郁稍。該數(shù)組將作為先前定義的 showRoseWind() 函數(shù)的第一個參數(shù)赦政，這個函數(shù)是用來繪制極區(qū)圖的。
showRoseWind(RoseWind_Speed(df_ravenna),'Ravenna',max(hist))
所示的風(fēng)向頻率玫瑰圖表示風(fēng)速在 360 度范圍內(nèi)的分布情況耀怜。
