|
! H: ? x2 O, }/ H4 [* k
气象海洋领域,常常会涉及到大规模的数据处理,比如高时空分辨率的模式数据、雷达数据和卫星数据等。 {- d; R: G8 H+ C1 x3 F* Y
% O! _* X1 {* y. |0 | 大部分情况,我们一般只会涉及小规模数据处理,对计算效率并不会太过追求。但是当数据量变大时,低效的数据处理所耗费的时间非常明显,因此高效的数据处理方式尤为重要。
# ]" t3 K0 y. w, A. V$ E9 R$ w. J3 O2 \4 _; J4 }
本篇以拟合一个高维数据的正太分布参数为例,介绍如何使用xarray+dask加速数据处理。 $ a, a3 ^9 g) H; a$ Z1 T
1 q; U1 f+ L; I8 `1 B) I+ ?# W 数据维度X[time, lev, lat, lon],需要对三维空间每一点,沿着时间维度做正太分布拟合(正太分布拟合只是作为例子,这里可以定义你需要的操作函数)。
' _6 |6 f* l5 t; H! f 其中几个关键点解释一下:
3 F0 j2 a2 s) f- d: X3 }8 G; p (1)首先定义拟合正太分布的函数
0 Y5 K2 v8 C8 X7 @: W8 S def norm_fit(data): ^% [, C S1 \0 Y0 T2 r2 ^
loc, scale = st.norm.fit(data)
$ I9 w/ L8 o7 T2 C% B7 M# m return np.array([loc, scale]) ; I- m, M8 ^4 r& Y* S
这里需要注意的是,拟合的函数,输出参数,需要打包为一个数组。并且它的维度需要和后面aplly_ufun定义的输出维度一致。 1 o. P3 i* v1 J* _
(2)xarray分块 0 A9 T# ~1 A [* u# A$ n$ |
x = xr.DataArray(rs.randn(500, 20, 1500, 1500), dims=["time", "lev", "lat","lon"])* V9 {* P4 J9 u4 E2 Q5 r
x = x.chunk({"lev":5, "lat":100,"lon":100})
, j4 l( [# `2 H xarray结合dask可以将一个大型数组分成一个个数据块(chunk),需要注意的是我们需要沿着时间维操作,拟合需要整个时间维度的数据,因此时间维time不能分块,只能对其他维度分块。 1 p. U. I# B$ ^% r
(3)xarray.apply_ufunc函数
. O4 d& `# K4 D( S result = xr.apply_ufunc(norm_fit,
( I, m# _( w h( y) Y x,' a4 q: d: ?$ r/ H3 N
input_core_dims=[["time"]],% {- f" A9 X% b, k
output_core_dims=[["paras"]],
+ X" J* Z8 w. ^) c7 b) @ dask="parallelized",) l! z. X! K4 o1 j2 R# a, H
output_dtypes=[np.float],8 l+ s6 j5 i; x
dask_gufunc_kwargs={output_sizes:{"paras":2}}
3 C: ]( }) {# G* f- l )
- l( j6 T7 ~3 D( U apply_ufunc函数具体可以参考官网教程,这里只说使用时遇到的困难,即如何定义输出维度:输出维度是用output_core_dim定义的,将输出的拟合参数(期望和标准差)定义为paras维度,维度的大小为2,通过output_sizes参数设置。这样我们输入[time, lev, lat, lon]的数据,在每个空间点对时间维度拟合之后,输出的数据为[lev, lat, lon, paras]。(PS: 这里感谢深雨露大佬的指点) . e0 f# {3 ]9 L2 F& \ \$ L
" P) B! q! H9 t$ N4 O) I( ]
以下是全部代码: 5 Z) t; H! \* i1 s/ U5 g8 C8 A# D
from scipy import stats as st7 n4 H/ f1 P. ~
$ l, N1 O+ D4 [0 C import xarray as xr
" g; b; A* ]3 e$ ]( r/ i import dask& ]) b6 j: w$ N4 W; h$ B
import numpy as np
7 H5 \+ I! {1 U: Z+ r from dask.diagnostics import ProgressBar& i2 Z$ r v5 d* K) k
. q8 l4 h7 g: n* b
def norm_fit(data):- g& p/ a5 b* F: x# R' [* u s
loc, scale = st.norm.fit(data). L) p- X0 H1 L" g1 _( o. Z% \7 d4 ]6 L
return np.array([loc, scale])0 U0 o* k# `& ?+ u0 Z: F W% T
/ E# W4 @0 I' C0 u( j% `4 l rs = np.random.RandomState(0). d9 t; W+ T. ?
x = xr.DataArray(rs.randn(500, 20, 1500, 1500), dims=["time", "lev", "lat","lon"])
z) [3 F! E# F$ O/ l x = x.chunk({"lev":5, "lat":100,"lon":100})
& ^$ _- ]! {0 M8 ?0 o" R+ R8 i' w8 O5 M5 L9 [! j
#使用apply_ufunc计算,并用dask的并行计算& t( O( _0 | u* m8 K
result = xr.apply_ufunc(norm_fit,
4 t. o5 G" V/ S x,/ ]9 I8 t: q7 ~9 @" \* o
input_core_dims=[["time"]],
6 b1 k7 ]4 Y! T+ e2 B output_core_dims=[["paras"]],8 h: E6 F" G, o/ ^; A+ H6 s
dask="parallelized",
7 a. X" a8 q4 t output_dtypes=[np.float], e1 m4 F0 [* d. M
dask_gufunc_kwargs={output_sizes:{"paras":2}}
* l$ ~( U1 U- W; L ): g( m Y3 N5 h' ]8 O% s
- S& h" }# I% }; [% _/ Y7 v* v
#compute进行真实的计算并显示进度1 ^: j! W# T" ~- C# [6 `
with ProgressBar():( h5 s9 @& E4 f8 ?- {3 L. k* K2 c
result = results.compute()
+ R" w O3 N9 W; H3 Q9 z+ I5 H s0 i6 W- {* w
#结果冲命名保存到nc文件! _6 p* O. G" Z: X
result = result.rename("norm_paras")
% j! C9 I+ Q) r2 n3 L result = result.to_netcdf("norm_fit_paras.nc",compute=False)3 N. U6 S' g+ z3 A, t$ j6 p( t
with ProgressBar():
" N2 E9 e* n" ]/ b+ d! l7 v result.compute()
, ?+ x& \/ U# r 转自:气海同途
# R0 m+ H7 C1 h3 b 关注【Ai尚研修科研技术平台】公众号,查看更多课程安排及免费数据资料
/ K6 n( ?) K$ ]$ W ! i1 J# @" p/ r9 Y' [1 E
推荐: - e& l# D* |* @
1、Python语言在地球科学领域中的应用实践技术应用
4 T+ Q; N$ t L1 R$ M 2、【夏令营】针对课题组人员AI培养计划:“开启AI科研之路” - r+ [# e/ M2 N; A3 {" `1 C. n
3、Python在气象与海洋中的实践技术应用精品课程
4 @- R$ X: A- e$ M4 @2 Y 4、Python在WRF模型自动化运行及前后处理中的实践技术应用
, Q6 D, H) p: |& e# d4 b 5、全套Python机器学习核心技术与案例分析实践应用视频
' ]) \' J5 f4 B4 D3 h, n5 y: u 6、全套区域高精度地学模拟-WRF气象建模、多案例应用与精美制图精品课程 ' m) _! t ?3 B( q
7、WRF DA资料同化系统理论、运行与与变分、混合同化新方法技术应用视频教程
' d. @6 N/ {) t; ~. H% h8 ~& {( c- W* \3 r; D: S
8 ]& I4 z% t* x) c! T9 i! z; w# L1 `( ^. N
0 f( [( G4 u; M* U8 R0 Z# E
|