Python气象数据分析

记录《Python数据分析实战》一书中关于意大利北部沿海地区气象数据分析的练习。

此次分析的目的是验证靠海对气候的影响，因此，选取10个城市分析他们的天气数据，其中5个城市距离海100公里以内，另外5个城市距离海100~400公里距离。此外，为了避免山区气候对天气数据造成影响，选取的城市均来自平原地区。

1. 加载数据集

# 导入模块
import numpy as np
import pandas as pd
import datetime

#导入数据可视化模块
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from dateutil import parser

# 加载数据集
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')

2. 温度数据分析

a. 气温走势

分析一天中气温变化趋势，以米兰为例

# 分析一天中气温变化趋势，x轴为时间，y轴为气温
y1 = df_milano['temp']
x1 = df_milano['day']

# 由于原始数据中【day】列中时间格式为：2015/6/27  1:54:14 PM
# 用dateutil.parser将日期数据转换成datetime格式
day_milano = [parser.parse(x) for x in x1]

# 调用subplot函数， fig是图像对象，ax是坐标轴对象
fig, ax = plt.subplot()

# 将x轴坐标刻度标签逆时针旋转70度，方便查看
plt.xticks(rotation=70)

# 设置时间格式为24：00（时：分）('%H:%M')
hours = mdates.DateFormatter('%H:%M')

# 设置x轴时间显示的格式为以上设置的格式
ax.xaxis.set_major_formatter(hours)

# 作折线图
ax.plot(day_milano,y1,'r')

由上图可知，气温走势接近正弦曲线，从早晨开始，气温逐渐上升，最高温出现在下午6点钟，其中30度以上的气温出现在下午3点到7点之间，随后气温逐渐下降，在次日早上6点达到最低。

b. 多城市气温走势对比

选取三个最靠近海的城市：ravenna, faenza, cesena和三个最不靠海的城市：milano, asti, torino

ax.plot(x,y)作折线图对比。

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']

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)

ax.plot(day_ravenna, y1, 'r', day_faenza, y2,'r', day_cesena, y3, 'r', day_milano, y4, 'g', day_asti, y5, 'g', day_torino, y6, 'g')

三个靠海城市用红色折线表示，三个不靠海城市用绿色折线表示。

由上图可知，靠海最近的三个城市最高气温比靠海最远的三个城市低很多，同时最低气温则差别很小。

因此，可以将10个城市的最高气温、最低气温与离海的距离作对比，进一步研究。

c. 横轴为距离，纵轴为气温，画10个城市折线图对比

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 = [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, 'r--o')

由上图可以看到：

1. 在这个时间段内，海洋对气温有一定的影响

2. 进一步，海洋对气温的影响衰减快，在20~70公里影响显著，在70公里以外的城市都处于气温高位。

d. 海洋对气温的影响存在一个强影响区（近海的5个城市）和弱影响区（远海的5个城市），因此可以拟合两条变化直线表示不同气温趋势。

使用sklearn.svm里的SVR拟合直线

# 用线型回归算法得到两条直线，分别表示不同的气温趋势
from sklearn.svm import SVR
# 两条直线分别是：离海较近的五个城市温度变化，离海较远的五个城市温度变化
dist1 = dist[:5]
dist2 = dist[5:10]
# 此处dist1输出为[8,14,37,47,71]为1维，shape为(1，N)，而fit()和predict()函数都只接收二维shape为(N,1)的DF
dist1 = [[x] for x in dist1]
dist2 = [[x] for x in dist2]
print(dist1)
temp_max1 = temp_max[:5]
temp_max2 = temp_max[5:]
print(temp_max1)
# 调用SVR函数，使用线型拟合
# 因为目的是线型拟合，拟合越精确越好，不需要考虑泛化能力，因此C可以设置很大
svr_lin1 = SVR(C=1000,kernel='linear')
svr_lin2 = SVR(C=1000, kernel='linear')
# 加入数据进行拟合
svr_lin1.fit(dist1, temp_max1)
svr_lin2.fit(dist2, temp_max2)
# 拟合后，对新的数据集进行预测，方便生产拟合直线的参数
# 因为predict()只接收shape为(N,1)的数据集，用reshape()将(1,N)变为(N,1)
xp1 = np.arange(10,100,10).reshape((-1,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')
# 显示拟合图
fig

由上图可知，离海60公里以内，气温上升速度很快，由28度升到31度，随后增速减缓60-400公里才上升1度。

要想知道这个海洋影响由强转弱转变点，需算出两直线公式y = ax +b，两直线相交点即为影响转变点。

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

[[0.04794118]]
[27.65617647]
[[0.00401274]]
[29.98745223]

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的坐标函数
# fsolve(func(x), x0) fsolve求解func(x)=0时x的值，其中x0为初始值，若求解的x有多个值，寻找接近x0的值。
def findIntersection(fun1,fun2,x0):
    return fsolve(lambda x: fun1(x)-fun2(x),x0)
result = findIntersection(line1, line2, 0.0)
# 格式符%d 返回十进制整数
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]

可以得到，海洋对气温影响显著的区域为离海距离53公里内。

e. 分析最低气温

# 规定x轴和y轴取值范围
plt.axis((0,400,15,25))
plt.plot(dist,temp_min,'bo')

由此，可知在8月26日这一天之中，海洋对各地最低气温并不大，似乎并不像最高温那样，海洋对气温有个缓和的作用。但要知道，此数据来自夏日的某天，或许冬日的某天海洋对最低温的影响会有所不同。

3. 湿度数据分析

a. 湿度整体趋势对比

# 选取离海最近的三个城市湿度数据和离海最远的三个城市湿度数据对比
# 读取湿度数据
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)

# parse时间数据
# 将时间从str类型转换为标准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]

# 设置坐标轴中时间显示格式为时：分（Y,m,d,H,M,S)(年月日时分秒)
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')

由上图整体趋势可以看出，似乎离海较近三个城市湿度普遍比厉害较远的三个城市湿度高20%。

b. 进一步，分析其湿度极值与离海远近距离的关系

# 最大湿度与距离关系
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')

由上图可以看出，无论是最大湿度还是最小湿度，离海较近的城市湿度均高于离海较远的城市。

但无法说明湿度与距离之间存在线性关系或其他曲线关系（需要更多城市样本量支撑）

4. 风向频率玫瑰图

在气象数据中，与风有关的数据有：风速和风向

若直接使用散点图

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

如上图使用笛卡尔坐标系作图，对于风向数据并不直观，可以考虑使用极坐标

a. 绘制极区图

首先创建一个直方图，将360度分为八个面元，每个面元为45度，把所有数据分到这八个面元中

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

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

以上histogram()函数返回的结果中，hist表示落在每个面元的数据点的数目，bins表示将【0，360】分割成8个面元的9条边界

以ravenna为例

N = 8
# 数据角度
theta = np.arange(0.+np.pi/8,2*np.pi+np.pi/8,2*np.pi/8)
# 数据极径
radii = np.array(hist)
# 绘制极区图坐标系,此处的[0,0,1,1]为极坐标图所在位置下面坐标limit[xmin,ymin,xmax,ymax]
plt.axes([0,0,1,1],polar=True)
# 定义每个扇区的RGB值（R,G,B），x越大，对应的颜色越接近蓝色
colors = [(1-x/max(hist), 1-x/max(hist),0.75) for x in radii]
plt.bar(theta,radii,width=(2*np.pi/N),bottom=0.0,color=colors)
plt.title("Ravenna",x=0.2,fontsize=20)

由上图可知，城市Ravenna一天中大部分时间的风都是西偏南，其次是西偏北方向。

定义一个风向玫瑰图函数，可以方便绘制不同城市的风向图:

def showRoseWind(values, city_name, max_value):
    N = 8
    theta = np.arange(0.+np.pi/8,2*np.pi+np.pi/8,2*np.pi/8)
    # 数据极径
    radii = np.array(values)
    # 绘制极区图坐标系
    plt.axes([0,0,1,1],polar=True)
    # 定义每个扇区的RGB值（R,G,B），x越大，对应的颜色越接近蓝色
    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))

# 查看都灵Torino风向
hist, bin = np.histogram(df_torino['wind_deg'],8,[0,360])
print(hist)
showRoseWind(hist, 'Torino',max(hist))

由上图可以看出，都灵在这一天中主要风向集中在北风，其次是西偏北。

b. 统计完了一天中风向的集中性，那么风速有什么影响呢

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)

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

由上图可知，一天中，城市Ravenna风速最大在西偏北平均4.5级，其次为西偏南平均3级。