本文目的- 介绍了如何从nc文件中,提取风速数据;
- 介绍如何将风速数据转换成时间序列;
- 简单的时间序列的趋势拆解(首发)。
6 [# ~$ U- c& d, H1 O9 k& Y7 }! i+ E& V) N
, H, `# [' Q) j# v3 ]) p 代码链接代码我已经放在Github上面了,免费分享使用,https://github.com/yuanzhoulvpi2 ... ree/main/python_GIS。 , |6 J5 g1 m+ E+ j* U0 ^
过程介绍
0 f3 L, q. k7 L5 _2 T* G. k C! Z% @2 s
+ v) v- I9 T: T' X# _6 S! ~ D: c
1. 导入包1 }% L$ j4 b, R* Q+ M2 O& d3 Z! B
% C+ i: S" C6 @, N5 `- N0 b% ~[Python] 纯文本查看 复制代码 # 基础的数据处理工具
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt # 可视化
import datetime # 处理python时间函数
import netCDF4 as nc # 处理nc数据
from netCDF4 import num2date # 处理nc数据
import geopandas as gpd # 处理网格数据,shp之类的
import rasterio # 处理tiff文件
from shapely.geometry import Point # gis的一些逻辑判断
from cartopy import crs as ccrs # 设置投影坐标系等
from tqdm import tqdm # 打印进度条
from joblib import Parallel, delayed # 并行
import platform # 检测系统
tqdm.pandas()
# matplotlib 显示中文的问题
if platform.system() == 'Darwin':
plt.rcParams["font.family"] = 'Arial Unicode MS'
elif platform.system() == 'Windows':
plt.rcParams["font.family"] = 'SimHei'
else:
pass / p$ ]9 c4 a7 O8 ?! }! d% U; c% Q
5 I4 B$ U: x& ^/ b: l. {- S
: D; ]- e' v# |* V2.导入数据 处理数据7 a& ~6 @! ?6 k$ W2 _0 c* Q
" I, M/ ~6 e5 Y6 g X$ }7 u( h5 A
) |4 M! G: }% _7 d% H9 S5 J[Python] 纯文本查看 复制代码 # 导入数据
nc_data = nc.Dataset("./数据集/GIS实践3/2016_2020.nc")
# 处理数据
raw_latitude = np.array(nc_data.variables['latitude'])
raw_longitude = np.array(nc_data.variables['longitude'])
raw_time = np.array(nc_data.variables['time'])
raw_u10 = np.array(nc_data.variables['u10'])
raw_v10 = np.array(nc_data.variables['v10'])
# 提取缺失值,并且将缺失值替换
missing_u10_value = nc_data.variables['u10'].missing_value
missing_v10_value = nc_data.variables['v10'].missing_value
raw_v10[raw_v10 == missing_v10_value] = np.nan
raw_u10[raw_u10 == missing_u10_value] = np.nan
# 处理时间
def cftime2datetime(cftime, units, format='%Y-%m-%d %H:%M:%S'):
"""
将nc文件里面的时间格式 从cftime 转换到 datetime格式
:param cftime:
:param units:
:param format:
:return:
"""
return datetime.datetime.strptime(num2date(times=cftime, units=units).strftime(format), format)
clean_time_data = pd.Series([cftime2datetime(i, units=str(nc_data.variables['time'].units)) for i in tqdm(raw_time)])
clean_time_data[:4] ) v; R6 G8 s, q8 ~$ p- F% m7 f7 B
" X: r5 g- y5 G7 d
3. 计算风速数据
) {2 [: S( X1 S! L1 ]5 h8 ^( P# d5 G+ H3 i5 ?$ ]& i
5 e% p( V4 Y! P$ M& f( g[Python] 纯文本查看 复制代码 windspeed_mean = pd.Series([np.sqrt(raw_v10[i,:, :] ** 2 + raw_u10[i, :, :]**2).mean() for i in tqdm(range(clean_time_data.shape[0]))])
time_windspeed = pd.DataFrame({'time':clean_time_data,'mean_ws':windspeed_mean})
time_windspeed
J& e8 @9 k+ G9 L- g5 @9 O
q/ Z- R$ t2 V0 D# |0 O4 e) o, y q+ Q$ R7 L; h
" d0 I; W6 S8 z( i" K' v& R6 Q a
4. 年度数据可视化$ ], L5 U% ]) F! ]* M, X4 m' |
W9 f2 i9 x. ]$ V
9 X& U5 V! y9 d8 S( Z[Python] 纯文本查看 复制代码 year_data = time_windspeed.groupby(time_windspeed.time.dt.year).agg(
mean_ws = ('mean_ws', 'mean')
).reset_index()
# year_data
with plt.style.context('fivethirtyeight') as style:
fig, ax = plt.subplots(figsize=(10,3), dpi=300)
ax.plot(year_data['time'], year_data['mean_ws'], '-o',linewidth=3, ms=6)
ax.set_xticks(year_data['time'])
#
#
for i in range(year_data.shape[0]):
ax.text(year_data.iloc[/size][/font][i][font=新宋体][size=3]['time']+0.1, year_data.iloc[/size][/font][i][font=新宋体][size=3]['mean_ws'], str(np.around(year_data.iloc[/size][/font][i][font=新宋体][size=3]['mean_ws'], 2)),
bbox=dict(boxstyle='round', facecolor='white', alpha=0.5))
#
for i in ['top', 'right']:
ax.spines[/size][/font][i][font=新宋体][size=3].set_visible(False)
ax.set_title("各年平均风速")
ax.set_ylabel("$Wind Speed / m.s^{-1}$")
% }1 a( f/ o- T8 u: Z' Z7 T
4 e. _/ ?4 g; `) s
+ j/ d1 Y, g T: o J. I# X. s! `2 ^
5. 月维度数据可视化
, c& P- R6 B; I1 x7 d5 C[Python] 纯文本查看 复制代码 month_data = time_windspeed.groupby(time_windspeed.time.dt.month).agg(
mean_ws = ('mean_ws', 'mean')
).reset_index()
with plt.style.context('fivethirtyeight') as style:
fig, ax = plt.subplots(figsize=(10,3), dpi=300)
ax.plot(month_data['time'], month_data['mean_ws'], '-o',linewidth=3, ms=6)
ax.set_xticks(month_data['time'])
_ = ax.set_xticklabels(labels=[f'{i}月' for i in month_data['time']])
for i in range(month_data.shape[0]):
ax.text(month_data.iloc[/size][/font][i][font=新宋体][size=3]['time'], month_data.iloc[/size][/font][i][font=新宋体][size=3]['mean_ws']+0.05, str(np.around(month_data.iloc[/size][/font][i][font=新宋体][size=3]['mean_ws'], 2)),
bbox=dict(boxstyle='round', facecolor='white', alpha=0.5))
for i in ['top', 'right']:
ax.spines[/size][/font][i][font=新宋体][size=3].set_visible(False)
ax.set_title("各月平均风速")
ax.set_ylabel("$Wind Speed / m.s^{-1}$")
fig.savefig("month_plot.png") . \7 i, G. [) _) h( w: b
* S! H& g! P4 [! _! N4 D- }
& X' m+ e* @" r; K 1 Z" n9 |7 K& Z3 @% s+ W7 r
6.天维度数据可视化4 d1 K. g1 ^6 m/ y2 H
- 计算天数据" @- a( |2 [; }6 D0 f1 `5 T Y- R
# B- G c( M& U* f$ U, `
[Python] 纯文本查看 复制代码 day_data = time_windspeed.groupby(time_windspeed.time.apply(lambda x: x.strftime('%Y-%m-%d'))).agg(
mean_ws = ('mean_ws', 'mean')
).reset_index()
day_data['time'] = pd.to_datetime(day_data['time'])
day_data = day_data.set_index('time')
day_data.head() - 可视化0 b: G, I- X( {: |
% T" n: Y7 l1 d* J/ R# E- }
[Python] 纯文本查看 复制代码 # day_data.dtypes
fig, ax = plt.subplots(figsize=(20,4), dpi=300)
ax.plot(day_data.index, day_data['mean_ws'], '-o')
# ax.xaxis.set_ticks_position('none')
# ax.tick_params(axis="x", labelbottom=False)
ax.set_title("每天平均风速")
ax.set_ylabel("$Wind Speed / m.s^{-1}$")
ax.set_xlabel("date")
fig.savefig('day_plot.png')
% p1 [2 U4 }7 N9 u
! n4 k& `" a3 f- c- y4 E R
/ E6 A1 E0 d! k3 o' r3 `! O) [4 V7 I( r2 Y& k
1.天维度数据做趋势拆解. }, l+ _1 `/ k# _5 ?8 a
! P3 t, `3 h8 ?9 ]) z7 ?) C0 j[Python] 纯文本查看 复制代码 # 导入包
from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse
# 乘法模型
result_mul = seasonal_decompose(day_data['mean_ws'], model="multilicative", extrapolate_trend='freq')
result_add = seasonal_decompose(day_data['mean_ws'], model="additive", extrapolate_trend='freq')
font = {'family': 'serif',
'color': 'darkred',
'weight': 'normal',
'size': 16,
}
# 画图
with plt.style.context('classic'):
fig, ax = plt.subplots(ncols=2, nrows=4, figsize=(22, 15), sharex=True, dpi=300)
def plot_decompose(result, ax, index, title, fontdict=font):
ax[0, index].set_title(title, fontdict=fontdict)
result.observed.plot(ax=ax[0, index])
ax[0, index].set_ylabel("Observed")
result.trend.plot(ax=ax[1, index])
ax[1, index].set_ylabel("Trend")
result.seasonal.plot(ax=ax[2, index])
ax[2, index].set_ylabel("Seasonal")
result.resid.plot(ax=ax[3, index])
ax[3, index].set_ylabel("Resid")
plot_decompose(result=result_add, ax=ax, index=0, title="Additive Decompose", fontdict=font)
plot_decompose(result=result_mul, ax=ax, index=1, title="Multiplicative Decompose", fontdict=font)
fig.savefig('decompose.png')
1 f# e$ Y2 v4 g, E- E9 V
% Z/ S/ M, H( }
+ _* R! ~) {, E8 M
! d& {* m1 s. S. i6 `+ _
& S G, W. j$ E g: c7 ~
5 W- y w: [: R8 S |