收藏本站 劰载中...网站公告 | 吾爱海洋论坛交流QQ群:835383472

[数据处理] 从nc文件中提取风速数据并且进行时间序列分析

[复制链接]
                                   本文目的
  • 介绍了如何从nc文件中,提取风速数据;
  • 介绍如何将风速数据转换成时间序列;
  • 简单的时间序列的趋势拆解(首发)。) Y( l0 U" \  K- f- o) `9 a! B

    : i3 y! k# F7 M' w8 _
代码链接

代码我已经放在Github上面了,免费分享使用,https://github.com/yuanzhoulvpi2 ... ree/main/python_GIS


" J5 {. T' `5 f3 w3 j! T  W5 C, V

过程介绍

6 c3 i$ l$ }9 t; b0 K6 f( n  E+ S

/ f4 o7 g7 [5 X8 {3 w0 h% Y  d
9 A/ n8 {: v7 r( i  W7 e
1. 导入包
; [2 T4 Y/ i3 s! a7 s2 C8 W7 @8 n
[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

* g& J$ S; t  f) k! K1 r% ~2 B
4 u7 B. p) m2 D& ~

' |/ u- S/ P$ I9 ^6 E+ ^# H. w2.导入数据 处理数据7 r0 C& R/ w  t2 E3 ^
6 h! J9 [5 g: w0 F5 _

* `* y3 v1 Q5 l0 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]
6 R% `! b4 w' j! r7 U) Z
/ [+ k/ ?! @/ K, U
3. 计算风速数据
$ N" l# |  D( K# y6 j% |8 d. _
+ E/ J# s& A* D- f9 P: A6 |# N* P

8 L8 X6 @4 G/ i# g( D2 ]" K9 @0 M  A0 a
[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
6 b/ ]$ d: o+ w& _- R- g( `& m
6b7fd110a68e6d3fd40460ccdd7a810b.png

, {. p4 l9 \- R
. c) L9 T1 o8 f0 X5 ?- g, h4 l

' E4 p- O$ L  Y# ^& \4. 年度数据可视化
  Q& M! J: t: O, S1 @3 a$ q9 _  I. D  O/ w, }" K) A: X, K, B

$ E- B/ J$ @( ]% L
[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}$")
" M! c0 p/ Y. ?' E" r# G5 n, l
952d93a401a01cd1fa10be892b8b64d6.png
0 l+ z9 \) _. \0 t* h# w, [
4 k4 ?8 [3 |4 D# I

- l# B# G# H" ~1 w5. 月维度数据可视化
4 v4 `% w9 G4 U1 k
[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")
3 P6 {* H" P# l* q9 p% p. n/ R
a520cff3361647efbb668c89005a5570.png

+ N! I) N9 u. A

3 B0 N5 s1 G  D3 E) K' I

" D! O4 h4 W1 p0 d! q! c6.天维度数据可视化
/ X" E! y) C& ]3 u
  • 计算天数据. m5 n- H8 K# d) Q. E. {

    5 O/ s' ~- U" b* z. C
[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()
  • 可视化% Y$ J  N! _1 [5 _

    & n) R- o# |7 j. h
[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')

" ~; G* x% `* k5 X$ \  v
2 ^; k$ q2 J" D* v' v. g
& |# k: a/ _4 D
053571827f212c867e38f40c8aa49ca5.png
3 i4 h( H& e/ ]: V0 H% b3 j
1.天维度数据做趋势拆解8 ^0 X2 j8 z; X9 _# G" ~" B

6 M6 r/ B( p/ k0 C% c
[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')
* T; p  `' z" Z0 t: K7 p
cd8468c3910ecbcfac542ed3328df432.jpeg                 " M! F$ H% e' {7 G$ }# [0 g8 N
+ x% [" J4 A5 J8 C
  i/ ^8 @* N4 i4 G$ v8 r5 s

5 s' }# q! ^; i3 |+ ^# P5 i# c/ q& d# z; [' S
回复

举报 使用道具

相关帖子

全部回帖
暂无回帖,快来参与回复吧
懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
尖叫的土豆
活跃在3 天前
快速回复 返回顶部 返回列表