上一节，我们介绍了xmitgcm的LLCreader入门操作，总体上这个工具可以通过网络自动加载门户网站（ECCO）的模式数据，实现对模式LLC4320（PB级）和LLC2160的读取。但还有许多问题尚未解决： - • 能否保存数据到本地
- • 能否按照经纬范围选取数据
- • 能否选择连续的一段时间（如70天，每小时一个）
- • 能否做较复杂计算! T- Z% K# L" I/ S" a* i* M; V
本节先介绍怎么通过LLCreader将数据下载本地，并使用GMT绘图。
# N. j/ j# o5 F
步骤加载llcreader
+ R8 S5 D X- S6 R* `. c Iimport xmitgcm.llcreader as llcreader
9 \' j$ P4 |- ^1 B9 o z% g; i6 P3 S: G. m- T& k% D; V3 A
初始化模型因为LLC4320数据体量较大，这里先以LLC2160为例。我们选择Theta进行导出测试。 model = llcreader.ECCOPortalLLC2160Model(): v W4 e. U1 @$ @: R K" @. ~' a- e
1 N, |9 P9 ~5 V% [& W- F C2 Hds_sst = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon',+ v* c6 t0 \1 S0 F( r5 B0 E; ?8 @
iter_start=model.iter_start,
4 J3 \6 Q" _/ ~; _# G iter_stop=(model.iter_start + model.iter_step),
! g0 v9 V8 Q3 c: w( P! G read_grid=True)& L; P6 @$ C, q# e
ds_sst.Theta ds_sst的数据结构为Xarray的dataset，ds_sst.Theta的数据结构为Xarray的dataarray。可以看到Theta的数据大小是231M，这只是一个1小时快照的单层全球数据，体积不算太大。
& M1 u. j0 b, t! X0 m上图为其coordinates，这是dataarray数据结构的一个关键术语，其中黑色粗体i,j,k,time表示维度坐标(Dimension coordinate)，其他为无维度坐标(Non-dimension coordinate)。 1 x, F; x2 y: U4 m$ M
绘图查看/ x& d; T9 X. m7 J
) ~$ o7 j3 [( o6 F( I$ Ytemp.Theta.plot()维度坐标i,j作为其x和y轴，另一维度time和k(深度0)出现在标题上。 因为LLC为特殊的网格结构，这里暂时不能直接实现地理绘图。后面会对此进行单独介绍。
' C* D9 q( k; d* p& o8 ~% S0 J" K0 U' D
保存 .to_netcdf()
2 J! K; S) p" x3 y) _ ^to_netcdf()是Xarray内置的函数，可将数据保存为nc格式。注意数据保存和llcreader函数无关，llcreader的主要作用是在线获得LLC模式数据，并转化成Xarray特有的数据结构，后面所有的操作都是使用Xarray，因此有必要对Xarray做全面的学习。后面会对此介绍。 ds_sst.Theta.to_netcdf('myfile.nc')上述命令可以将Theta参数保存。通过panoply可以绘图展示SST参数的全球分布： 删除不重要的坐标
1 q9 b1 r0 j: E; O& ]但可以发现保存的文件大小是2.5G，远远超过了其自身大小，这是什么原因呢？ # U% h2 E x! U0 s9 c
原来是把无维度坐标全部保存了，这些坐标其实也是数据，因此占用很大空间。一共有11个这样的2D坐标，再加上Theta，数据体积就达到2.5G了。
6 D7 z4 S0 o" t9 d; z$ ^" i/ }: V! W3 A% Y; h% w, H- S! q
这些无维度坐标可以删去： ds_sst.Theta.reset_coords() # 将无维度坐标转为变量(不删除)
9 p" R/ R8 B( `3 ~3 J' |; nds_sst[["Theta"]].reset_coords(drop=True) # 将会删去无维度坐标若想保留某些坐标，比如经纬度XC,YC： ds_sst.Theta.reset_coords()[["XC", "YC","Theta"]]上面使用了[[ ]]做变量的筛选。 然后保存： ds_sst[["Theta"]].reset_coords(drop=True).to_netcdf('myfile2.nc')
& X" ?+ c' L B& G, H! Hds_sst.Theta.reset_coords()[["XC", "YC","Theta"]].to_netcdf('myfile3.nc')这样，数据体积就小了，其中myfile2.nc仅包含Theta，大小213M。 0 {' H6 q5 W% i
GMT绘图7 G x+ n; w6 W/ o2 E1 O! y$ _! v
GMT绘图比ython自带函数更加强大，既然已经下载到本地，则看一下GMT绘图效果： gmt begin ex02 png1 _4 m# R5 j4 b9 A+ x
gmt set MAP_ANNOT_OBLIQUE 0
/ Z; X4 a! f' q1 E gmt subplot begin 2x1 -A+JTL -Fs16c/10c -M0 -R1000/2000/2500/3500 -JX15/13c -B100 -T"Aghulhas Rings"
8 X W5 l/ B) u gmt subplot set 0,0 -Ce3c
+ Q- W5 B; K7 k9 d7 R3 @ gmt grd2cpt myfile3.nc?Theta -Cjet -Z
- W6 V9 O5 q0 j$ Z& C gmt grdimage myfile3.nc?Theta -I+a90 % x# A7 h2 f, u: @2 c' [7 a0 i
gmt colorbar -DJRM+o1.5c/0+mc -I0.3 -Bx2+lTheta -By+ldegree
6 ~. |0 N9 C+ P, t- k5 j
5 I H6 x: q M# N gmt subplot set 1,0 -Ce3c
0 P' O; ~2 z6 y6 u' T gmt grd2cpt myfile4.nc?Theta -Crainbow -Z
. ?/ D# x/ N# Z* y9 s% n) Z0 L7 S gmt grdimage myfile4.nc?Theta
* U* O% {+ w, A; i; T8 R5 k$ Q" {& n gmt colorbar -DJRM+o1.5c/0+e+mc -Bx2+lTheta -By+ldegree2 j v; N' W$ v# P7 }
gmt subplot end7 a, P2 Z4 `5 C( [# P
gmt end show " W; P% F, e/ B" O0 ?4 ~# c1 ^0 a6 ~
台湾海域下节预告- • 范围的筛选
- • 时间维度的连续导出1 d& j! v# O- Y7 U3 x+ h H, O8 K; E5 G
1 V2 G- k; ~( Q. P
| 
