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

xmitgcm下载LLC海洋模式数据

[复制链接]

上一节,我们介绍了xmitgcm的LLCreader入门操作,总体上这个工具可以通过网络自动加载门户网站(ECCO)的模式数据,实现对模式LLC4320(PB级)和LLC2160的读取。但还有许多问题尚未解决:
  • • 能否保存数据到本地
  • • 能否按照经纬范围选取数据
  • • 能否选择连续的一段时间(如70天,每小时一个)
  • • 能否做较复杂计算$ H9 V/ ?6 m! j
本节先介绍怎么通过LLCreader将数据下载本地,并使用GMT绘图。
. ~. R9 S) E7 {0 j
步骤加载llcreader: A: ^6 n$ h! Z! E
import xmitgcm.llcreader as llcreader
# w/ z( y- ~" _7 B! ]$ Z* [
5 D5 @. s. ~2 [/ t- [- x# [& s初始化模型
因为LLC4320数据体量较大,这里先以LLC2160为例。我们选择Theta进行导出测试。
model = llcreader.ECCOPortalLLC2160Model()
5 k2 j# x5 M1 r2 R
7 ^- i( z5 A, R  c9 v# C( Eds_sst = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon',5 A. Y# h& z) ]; ]. q8 i+ j1 P
                          iter_start=model.iter_start,. s0 d$ ~5 X7 V" v' i3 j
                          iter_stop=(model.iter_start + model.iter_step),
1 N6 j2 I" M3 U+ W9 E                          read_grid=True)
2 l4 g6 G2 D* U& hds_sst.Theta                          
e17c398cadb9e820fdac5f77c3133fe4.png
ds_sst的数据结构为Xarray的dataset,ds_sst.Theta的数据结构为Xarray的dataarray。可以看到Theta的数据大小是231M,这只是一个1小时快照的单层全球数据,体积不算太大。
58abfadff7631b95c513bd2043e38d04.png
, a* U" O) i* B: B' `
上图为其coordinates,这是dataarray数据结构的一个关键术语,其中黑色粗体i,j,k,time表示维度坐标(Dimension coordinate),其他为无维度坐标(Non-dimension coordinate)。

' C& X2 F, E, b) C3 C' O+ B% M0 s1 `
绘图查看
4 D7 z1 {: w0 g; @. M( ?; U
/ y/ z5 h$ s' wtemp.Theta.plot()
f270b69e1b4897a0cd47fc89f200ecb6.png
维度坐标i,j作为其x和y轴,另一维度time和k(深度0)出现在标题上。
因为LLC为特殊的网格结构,这里暂时不能直接实现地理绘图。后面会对此进行单独介绍。
% K' N0 Y) J& `  ~. a- D
保存 .to_netcdf()( ]% y# ]' V2 h) Q0 ^8 c0 |
to_netcdf()是Xarray内置的函数,可将数据保存为nc格式。注意数据保存和llcreader函数无关,llcreader的主要作用是在线获得LLC模式数据,并转化成Xarray特有的数据结构,后面所有的操作都是使用Xarray,因此有必要对Xarray做全面的学习。后面会对此介绍。
ds_sst.Theta.to_netcdf('myfile.nc')
上述命令可以将Theta参数保存。通过panoply可以绘图展示SST参数的全球分布:
91cc1704eec295e82dba5dc4b4ea2558.png
删除不重要的坐标* y; p0 H" r4 Q
但可以发现保存的文件大小是2.5G,远远超过了其自身大小,这是什么原因呢?
27bd7c12fcde88ccd1396cf073f4d373.png

+ t0 R$ v3 x7 `! ^+ d原来是把无维度坐标全部保存了,这些坐标其实也是数据,因此占用很大空间。一共有11个这样的2D坐标,再加上Theta,数据体积就达到2.5G了。2 {3 A! [; Z" O! L, L
4 j; `+ F9 c) B5 _( I  u  r: w+ |+ q
d403e8d05dabb82f32cb37999f6e03f8.png
这些无维度坐标可以删去:
ds_sst.Theta.reset_coords() # 将无维度坐标转为变量(不删除)
: {" A8 A  T# D( \; K% t/ q0 G* `ds_sst[["Theta"]].reset_coords(drop=True) # 将会删去无维度坐标
976ce67a006158504eefcf096d561940.png
若想保留某些坐标,比如经纬度XC,YC:
ds_sst.Theta.reset_coords()[["XC", "YC","Theta"]]
上面使用了[[ ]]做变量的筛选。
然后保存:
ds_sst[["Theta"]].reset_coords(drop=True).to_netcdf('myfile2.nc')% p1 l' h! D/ t. D, M# k
ds_sst.Theta.reset_coords()[["XC", "YC","Theta"]].to_netcdf('myfile3.nc')
这样,数据体积就小了,其中myfile2.nc仅包含Theta,大小213M。

3 v& I% i4 j- H3 E2 s
GMT绘图# H  Z7 ^" G' L# a6 d) C
GMT绘图比ython自带函数更加强大,既然已经下载到本地,则看一下GMT绘图效果:
gmt begin ex02 png
& O8 s6 ]% ?: r" h2 L# E    gmt set MAP_ANNOT_OBLIQUE 05 G; E" T. u9 g- m) e6 }1 N% t
    gmt subplot begin 2x1 -A+JTL -Fs16c/10c -M0 -R1000/2000/2500/3500 -JX15/13c -B100 -T"Aghulhas Rings", h- ]3 m) a5 c) G$ i5 Q. H9 I. C
        gmt subplot set 0,0 -Ce3c
+ |9 A6 X6 F& v$ \, M% V2 T9 j        gmt grd2cpt myfile3.nc?Theta -Cjet -Z
  H0 U; k6 b. z' k/ B' G        gmt grdimage myfile3.nc?Theta -I+a90  ; z0 Y1 l: m! P2 q
        gmt colorbar -DJRM+o1.5c/0+mc -I0.3 -Bx2+lTheta -By+ldegree* Z0 i- D; |# M3 a8 n; W5 t

/ F, t& S8 p) j. ^/ n. P        gmt subplot set 1,0 -Ce3c! Z! e# f  i. A2 l# i1 D
        gmt grd2cpt myfile4.nc?Theta -Crainbow -Z. d$ ^2 l/ i1 j9 w
        gmt grdimage myfile4.nc?Theta
5 [6 h' i; v' M8 [7 \) a$ v        gmt colorbar -DJRM+o1.5c/0+e+mc -Bx2+lTheta -By+ldegree3 l: Y/ A! _0 W- o
    gmt subplot end8 Z% P' A! v$ d+ Q  ]
gmt end show : y  |! C9 \1 x! O: U4 P- e% p
c8c573203e6afdb5b6b0dd2d6cb60576.png
台湾海域
1c9c67d88f9b8f004fc59c26c772f279.png
下节预告
  • • 范围的筛选
  • • 时间维度的连续导出
    ; o# l8 l* ?! ^" z2 j) k0 A

9 c" X8 U) K4 Q
回复

举报 使用道具

相关帖子

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