上一节，我们介绍了xmitgcm的LLCreader入门操作，总体上这个工具可以通过网络自动加载门户网站（ECCO）的模式数据，实现对模式LLC4320（PB级）和LLC2160的读取。但还有许多问题尚未解决： - • 能否保存数据到本地
- • 能否按照经纬范围选取数据
- • 能否选择连续的一段时间（如70天，每小时一个）
- • 能否做较复杂计算+ ~% G, V3 |: C; t/ S- X
本节先介绍怎么通过LLCreader将数据下载本地，并使用GMT绘图。 ) U' h# O* m& U1 E. T
步骤加载llcreader
2 s& G# D" p6 ~% c. G' H6 Cimport xmitgcm.llcreader as llcreader
$ @/ ?: Z' v$ l1 V- ]
v/ } Q- v9 l! w" d/ |3 G$ e0 x初始化模型因为LLC4320数据体量较大，这里先以LLC2160为例。我们选择Theta进行导出测试。 model = llcreader.ECCOPortalLLC2160Model()
' L1 B8 p5 h3 q6 h6 w% q2 H( U8 @+ n( A8 A/ a7 X6 Q
ds_sst = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon',
- l7 p6 q/ X8 f iter_start=model.iter_start,6 F0 }+ G+ t' _0 M! n
iter_stop=(model.iter_start + model.iter_step),
; _5 W1 O# ?5 [# C read_grid=True)
9 j- N1 `( f. E3 ~' `1 m% T# Uds_sst.Theta ds_sst的数据结构为Xarray的dataset，ds_sst.Theta的数据结构为Xarray的dataarray。可以看到Theta的数据大小是231M，这只是一个1小时快照的单层全球数据，体积不算太大。
: v/ A* f. H6 v$ @- C上图为其coordinates，这是dataarray数据结构的一个关键术语，其中黑色粗体i,j,k,time表示维度坐标(Dimension coordinate)，其他为无维度坐标(Non-dimension coordinate)。 1 z! f F; {' C5 s9 K
绘图查看" n$ e' U0 B1 C* f
, c1 u: `) A; m0 R/ ]
temp.Theta.plot()维度坐标i,j作为其x和y轴，另一维度time和k(深度0)出现在标题上。 因为LLC为特殊的网格结构，这里暂时不能直接实现地理绘图。后面会对此进行单独介绍。 / K( l1 }% M! r* J
保存 .to_netcdf()
) e5 t7 |3 ~, L* Ato_netcdf()是Xarray内置的函数，可将数据保存为nc格式。注意数据保存和llcreader函数无关，llcreader的主要作用是在线获得LLC模式数据，并转化成Xarray特有的数据结构，后面所有的操作都是使用Xarray，因此有必要对Xarray做全面的学习。后面会对此介绍。 ds_sst.Theta.to_netcdf('myfile.nc')上述命令可以将Theta参数保存。通过panoply可以绘图展示SST参数的全球分布： 删除不重要的坐标8 k+ e2 E5 J' J* n( Z0 F, j# z0 ?
但可以发现保存的文件大小是2.5G，远远超过了其自身大小，这是什么原因呢？ 2 |4 C0 n( h3 N( m7 ^ }) R& {$ T
原来是把无维度坐标全部保存了，这些坐标其实也是数据，因此占用很大空间。一共有11个这样的2D坐标，再加上Theta，数据体积就达到2.5G了。- r' U" [! T; l
; w: H5 P' y0 m& r& J) ^这些无维度坐标可以删去： ds_sst.Theta.reset_coords() # 将无维度坐标转为变量(不删除)
]. }6 Z0 W! b& s6 U/ Uds_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')2 y. r: l V! Q* p7 i
ds_sst.Theta.reset_coords()[["XC", "YC","Theta"]].to_netcdf('myfile3.nc')这样，数据体积就小了，其中myfile2.nc仅包含Theta，大小213M。
2 ^: W' L7 E; G
GMT绘图
4 N# V) U2 F0 }* D# dGMT绘图比ython自带函数更加强大，既然已经下载到本地，则看一下GMT绘图效果： gmt begin ex02 png) [% a7 _1 Q5 X$ [2 N4 {2 R' v
gmt set MAP_ANNOT_OBLIQUE 0
$ h4 J; v9 _: _& |6 j gmt subplot begin 2x1 -A+JTL -Fs16c/10c -M0 -R1000/2000/2500/3500 -JX15/13c -B100 -T"Aghulhas Rings"
, D* h- F5 P# n7 P2 ^ gmt subplot set 0,0 -Ce3c- n. Z. ]( R2 b8 j; G4 `
gmt grd2cpt myfile3.nc?Theta -Cjet -Z
( w& ^( z1 [# x$ \& S/ J: ] gmt grdimage myfile3.nc?Theta -I+a90 8 `0 d, \( w8 _. P$ K! I- i
gmt colorbar -DJRM+o1.5c/0+mc -I0.3 -Bx2+lTheta -By+ldegree6 }: U! h, y$ d* u g
6 u; M3 p, A5 s+ \. ~. d; D gmt subplot set 1,0 -Ce3c
7 v$ w# S" O& j3 p gmt grd2cpt myfile4.nc?Theta -Crainbow -Z/ X+ w. _+ O J& t+ F
gmt grdimage myfile4.nc?Theta3 R7 N( [ e! e R: E& R, k% m0 E
gmt colorbar -DJRM+o1.5c/0+e+mc -Bx2+lTheta -By+ldegree6 \( H9 i9 H# ]. i
gmt subplot end
8 f& B6 J7 K' c2 U2 Q! Bgmt end show
5 T# C% z5 l- i台湾海域下节预告- • 范围的筛选
- • 时间维度的连续导出
- | V c; d ?) A. G% \
4 R' j; Z; `$ Y$ A2 c- L | 
