|
+ ~: k& i- L6 \8 d/ a- n) D
气象海洋领域，常常会涉及到大规模的数据处理，比如高时空分辨率的模式数据、雷达数据和卫星数据等。 5 ^! O2 a5 B- _ w P. C$ }5 N
* h, P% z4 c$ U( Q0 E8 o" f Y' T
大部分情况，我们一般只会涉及小规模数据处理，对计算效率并不会太过追求。但是当数据量变大时，低效的数据处理所耗费的时间非常明显，因此高效的数据处理方式尤为重要。
) V$ W% K5 |8 s) ~0 {
l3 P+ x0 \+ |( J$ o 本篇以拟合一个高维数据的正太分布参数为例，介绍如何使用xarray+dask加速数据处理。
: I( t4 `' \/ c
|( Q4 y- q6 m6 q0 Z7 w/ \- x8 ~) Y 数据维度X[time, lev, lat, lon]，需要对三维空间每一点，沿着时间维度做正太分布拟合(正太分布拟合只是作为例子，这里可以定义你需要的操作函数)。 ) {: ?+ p8 N9 P$ m8 i
其中几个关键点解释一下： $ I. S% Z+ b* C9 c
（1）首先定义拟合正太分布的函数
, f3 E2 b; z" @/ B# | def norm_fit(data):# V v' V- b7 [9 p$ ]
loc, scale = st.norm.fit(data)) F: G/ x1 O5 ~; h
return np.array([loc, scale])
/ @& `4 [. q% ^7 o 这里需要注意的是，拟合的函数，输出参数，需要打包为一个数组。并且它的维度需要和后面aplly_ufun定义的输出维度一致。 ( t- v: ~- Y; S+ u3 v2 @" }
（2）xarray分块 8 \* s7 R2 B- p- ~2 [
x = xr.DataArray(rs.randn(500, 20, 1500, 1500), dims=["time", "lev", "lat","lon"])
4 L+ z( A6 O- ~! ?. @: x x = x.chunk({"lev":5, "lat":100,"lon":100})
1 V1 f5 P, r1 U xarray结合dask可以将一个大型数组分成一个个数据块（chunk），需要注意的是我们需要沿着时间维操作，拟合需要整个时间维度的数据，因此时间维time不能分块，只能对其他维度分块。 9 k) ]$ T4 b6 x% |4 @& ?6 P) V
（3）xarray.apply_ufunc函数 9 ?! u5 O# _: H0 x
result = xr.apply_ufunc(norm_fit,
' ^4 a' t8 w4 Q, z9 Q' t- S1 s x,
- r1 _! n& j1 @0 H& j input_core_dims=[["time"]],
. j' ^2 H" l: I4 t output_core_dims=[["paras"]],
- y# @' \# Y& `* ]. b dask="parallelized",
, m& T! C# A' s. Z( _ output_dtypes=[np.float],
# h, T! @* x; \5 f dask_gufunc_kwargs={output_sizes:{"paras":2}}
0 ^; t4 K5 v3 t) [8 G ) . R9 G2 ?- Y, m2 N' |6 q" ~5 Z
apply_ufunc函数具体可以参考官网教程，这里只说使用时遇到的困难，即如何定义输出维度：输出维度是用output_core_dim定义的，将输出的拟合参数（期望和标准差）定义为paras维度，维度的大小为2，通过output_sizes参数设置。这样我们输入[time, lev, lat, lon]的数据，在每个空间点对时间维度拟合之后，输出的数据为[lev, lat, lon, paras]。（PS: 这里感谢深雨露大佬的指点） ! P; c1 {' \0 V/ P
3 t+ x; }, Q3 W) v) } 以下是全部代码：
* e' _: j! D8 q' F' r& W% O" W from scipy import stats as st
* t5 o O. W, E
& B4 h" L: E6 I4 H1 h) v+ H/ L import xarray as xr7 ^; E5 \! \8 _9 k9 @
import dask
7 I) C2 Q! y4 e, @( p' r& S. N8 i import numpy as np
4 b; C+ V [: o! O. c. j& `8 \ from dask.diagnostics import ProgressBar0 P; P$ E" a( _
# ^% E9 r" `" H
def norm_fit(data):5 y4 O0 p% M; h) n8 y! j
loc, scale = st.norm.fit(data)
% I2 `; a- Q$ A0 E2 J, A' U: g return np.array([loc, scale])
) C9 r$ K$ `. J! P2 C7 I# c6 N- Z3 V
rs = np.random.RandomState(0)
) m4 {% i5 y8 P& I0 Q4 ~. o x = xr.DataArray(rs.randn(500, 20, 1500, 1500), dims=["time", "lev", "lat","lon"]) A* `8 n* G; Q5 s' I- U* K
x = x.chunk({"lev":5, "lat":100,"lon":100})
2 W% q& i! y7 T% L& m
6 a y' d9 Z7 u1 a) ^4 k #使用apply_ufunc计算，并用dask的并行计算4 K& i( u* `: \ }+ w8 y5 b
result = xr.apply_ufunc(norm_fit,/ ?4 x/ _: l1 O1 A% l
x,. c, l) M1 L9 K m& K$ T
input_core_dims=[["time"]],
W e% r8 J2 Q. C output_core_dims=[["paras"]],, M; k# o: Q. K. o0 ]. a
dask="parallelized", L7 q) e4 K- K& P8 E& [
output_dtypes=[np.float],) W7 B/ @3 D0 e+ @( j' N0 p( t5 ]
dask_gufunc_kwargs={output_sizes:{"paras":2}}
! Z- a5 D! Z$ T. O# V )
! c$ S1 b& |4 e p- _) W' C5 c4 I) N: E5 q1 b! i# G* ?
#compute进行真实的计算并显示进度7 {" Q) q' t5 F
with ProgressBar():9 d! z( G: g* q0 E# V8 G
result = results.compute()
* L* b4 h4 X% \( O0 z3 X5 D; z& f0 K7 d! d/ j
#结果冲命名保存到nc文件4 P/ i0 A6 x! ~1 K# j0 ~! o
result = result.rename("norm_paras")% g$ q% t5 _0 P4 V8 G/ {
result = result.to_netcdf("norm_fit_paras.nc",compute=False)) |/ A8 z# p' _$ G! _
with ProgressBar():
\9 A% [: Q* r. @. [1 K result.compute()
! [; t3 D0 Y, {" v 转自：气海同途 6 R' U) O; H- p* O8 ]9 x3 q, o
关注【Ai尚研修科研技术平台】公众号，查看更多课程安排及免费数据资料
( i5 h9 b7 [# B2 k/ l- B9 Y2 h, H ! d. {0 E# H$ Q. T2 n" k
推荐： 7 k/ D$ ^4 |! O" z# R
1、Python语言在地球科学领域中的应用实践技术应用
) m8 f2 p& C) @. } 2、【夏令营】针对课题组人员AI培养计划:“开启AI科研之路”
. ]: U2 s1 J' d+ J* D 3、Python在气象与海洋中的实践技术应用精品课程 6 |6 e! T! O4 A
4、Python在WRF模型自动化运行及前后处理中的实践技术应用 . r* X' A. {0 s5 B( J* ]
5、全套Python机器学习核心技术与案例分析实践应用视频
' b) T7 V$ U9 { 6、全套区域高精度地学模拟-WRF气象建模、多案例应用与精美制图精品课程 ) t; Y; ]5 u2 \8 y) {
7、WRF DA资料同化系统理论、运行与与变分、混合同化新方法技术应用视频教程
/ c" H8 y+ I; u: w, }; c
& s; b0 o. \4 r7 | u' J* T/ D% q- W3 X. D
1 ~- u% U4 a- c9 A/ r9 E1 C+ S+ `$ J7 r* L/ e Q3 F
| 
