本实验将对意大利北部沿海地区的气象数据进行分析与可视化。我们在实验过程中先会运用 Python 中 matplotlib 库的对数据进行图表化处理，然后调用 scikit-learn 库当中的的 SVM 库对数据进行回归分析，最终在图表分析的支持下得出我们的结论。
下载气象数据
!wget -nc http://labfile.oss.aliyuncs.com/courses/780/WeatherData.zip
image.png
装 unzip 解压缩
!apt-get install unzip
image.png
解压缩
!unzip -o WeatherData.zip
image.png
这时候，我们通过 tree 命令应该能够再 WeatherData 中间看到 10 个城市的天气数据文件（以 .csv 结尾）
!apt-get install tree
!tree WeatherData/
image.png
导入相关包开始实验
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.pylot as plt
import matplotlib.dates as mdates
from dateutil import parser
#取出milano的温度和日期数据
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轴数据，temp是Y轴数据，‘r’代表红色
ax.plot(day_milano,y1,'r')
image.png
由图可见，气温走势接近正弦曲线，从早上开始气温逐渐升高，最高温出现在下午两点到六点之间，随后气温逐渐下降，在第二天早上六点时达到最低值。
我们进行数据分析的目的是尝试解释是否能够评估海洋是怎样影响气温的，以及是否能够影响气温趋势，因此我们同时来看几个不同城市的气温趋势。这是检验分析方向是否正确的唯一方式。因此，我们选择三个离海最近以及三个离海最远的城市。
#选择三个离海最近以及三个离海最远的城市。
#用循环会不会更方便……？
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条，远海的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')
image.png
离海最近的三个城市的气温曲线使用红色，而离海最远的三个城市的曲线使用绿色。
离海最近的三个城市的最高气温比离海最远的三个城市低不少，而最低气温看起来差别较小。
我们可以沿着这个方向做深入研究，收集 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(),]
#把最高温画出来,x轴为距离。y轴为最高温度
fig,ax=plt.subplots()
ax.plot(dist,temp_max,'ro')
image.png
进一步观察上图，你会发现海洋的影响衰减得很快，离海 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_max2同理
temp_max1=temp_max[0:5]
temp_max2=temp_max[5:10]
#调用SVR函数，规定为线性拟合
#并且把 C 设为1000来尽量拟合数据（因为不需要精确预测不用担心过拟合）
#C是什么？我并不清楚
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))
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')
image.png
这里
np.arange(10,100,10) 会返回 [10, 20, 30,..., 90]，如果把列表看成是一个矩阵，那么这个矩阵是 1 9 的。这里 reshape((9,1)) 函数就会把该列表变为 9 1 的， [[10], [20], ..., [90]]。这么做的原因是因为 predict() 函数的只能接受一个 N 1 的列表，返回一个 1 N 的列表。
如上所见，离海 60 公里以内，气温上升速度很快，从 28 度陡升至 31 度，随后增速渐趋缓和（如果还继续增长的话），更长的距离才会有小幅上升。这两种趋势可分别用两条直线来表示，直线的表达式为：
查看斜率和截距：
print(svr_lin1.coef_)#斜率
print(svr_lin1.intercept_)#截距
print(svr_lin2.coef_)
print(svr_lin2.intercept_)
image.png
你可能会考虑将这两条直线的交点作为受海洋影响和不受海洋影响的区域的分界点，或者至少是海洋影响较弱的分界点。
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)
x=np.linspace(0,300,31)
plt.plot(x,line1(x),x,line2(x),result,line1(result),'ro')
#result,line1(result),为交点坐标
print(result)
print(line1(result))
image.png
可以说海洋对气温产生影响的平均距离（该天的情况）为 53 公里
现在，我们可以转而分析最低气温。
#规定了x轴的范围，y轴的范围
plt.axis((0,400,15,25))
#绘出最低温度的图表
plt.plot(dist,temp_min,'bo')
image.png
在这个例子中，很明显夜间或早上 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']
# 调用 subplots() 函数，重新定义 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)
#这里需要画出三根线，所以需要三组参数， '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')
image.png
红色是近海城市，绿色为远海城市
乍看上去好像近海城市的湿度要大于内陆城市，全天湿度差距在 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')
image.png
#探寻最低湿度
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')
image.png
可以确定，近海城市无论是最大还是最小湿度都要高于内陆城市.
然而,我们还不能说湿度和距离之间存在线性关系或者其他能用曲线表示的关系。我们采集的数据点数量（10）太少，不足以描述这类趋势。
风向频率玫瑰图
在我们采集的每个城市的气象数据中，下面两个与风有关：
- 风力（风向）
- 风速
分析存放每个城市气象数据的 DataFrame 就会发现，风速不仅跟一天的时间段相关联，还与一个介于 0~360 度的方向有关。
例如，每一条测量数据也包含风吹来的方向。
df_ravenna[['wind_deg','wind_speed','day']]
image.png
为了更好地分析这类数据，有必要将其做成可视化形式，但是对于风力数据，将其制作成使用笛卡儿坐标系的线性图不再是最佳选择。
要是把一个 DataFrame 中的数据点做成散点图
plt.plot(df_ravenna['wind_deg'],df_ravenna['wind_speed'],'ro')
image.png
很显然该图的表现力也有不足。
要表示呈 360 度分布的数据点，最好使用另一种可视化方法：极区图。
首先，创建一个直方图，也就是将 360 度分为八个面元，每个面元为 45 度，把所有的数据点分到这八个面元中。
hist,bins=np.histogram(df_ravenna['wind_deg'],8,[0,360])
print(hist)
print(bins)
image.png
histogram() 函数返回结果中的数组 hist 为落在每个面元的数据点数量。
返回结果中的数组 bins 定义了 360 度范围内各面元的边界。
#绘制玫瑰风图
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)
# 绘制极区图的坐标系
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)
画图
showRoseWind(hist,'Ravenna',max(hist))
image.png
这里，扇形的颜色越接近蓝色，值越大.
整个 360 度的范围被分成八个区域（面元），每个区域弧长为 45 度，此外每个区域还有一列呈放射状排列的刻度值。
在每个区域中，用半径长度可以改变的扇形表示一个数值，半径越长，扇形所表示的数值就越大。
为了增强图表的可读性，我们使用与扇形半径相对应的颜色表。半径越长，扇形跨度越大，颜色越接近于深蓝色。
#定义好 showRoseWind() 函数之后，查看其他城市的风向情况也非常简单。
hist,bin=np.histogram(df_ferrara['wind_deg'],8,[0,360])
print(hist)
showRoseWind(hist,'Ferrara',max(hist))
image.png
计算风速均值的分布情况
即使是跟风速相关的其他数据，也可以用极区图来表示。
定义 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())
#这里 df_city[(df_city['wind_deg']>(deg-46)) & (df_city['wind_deg']<deg)]
#获取的是风向大于 deg-46 度和风向小于 deg 的数据。
return np.array(tmp)
RoseWind_Speed() 函数返回一个包含八个平均风速值的 NumPy 数组。该数组将作为先前定义的 showRoseWind() 函数的第一个参数，这个函数是用来绘制极区图的。
showRoseWind(RoseWind_Speed(df_ravenna),'Ravenna',max(hist))
image.png
所示的风向频率玫瑰图表示风速在 360 度范围内的分布情况。
网友评论