高分辨率海洋模式LLC4320的在线读取

LLC4320模式简介

LLC4320是MITgcm 1/48°全球海洋模式，其诞生的直接原因是为即将到来的SWOT任务提供高分辨率的全球海洋模拟。

LLC4320模式在多个方面具有开创性，特别是其高空间分辨率（全球分辨率在 1 到 2 公里之间）、潮汐驱动、高频（每小时）输出，其海面高度信号中包含了内潮、内波、地转平衡等信号。除了在SWOT相关工作取得应用之外，该模式也在海洋亚中尺度、内波等研究方向得到了广泛应用。

- U4 k% c3 m6 K, g

该模式主要特征是：

9 j) g! ]. f3 D- w/ ~5 t

• • 全球覆盖（包含极地）
• • 垂向90层
• • 分辨率1/48°
• • 全球海洋分成13个face，每个face的网格数为4320*4320
• • 时间采样是1小时，总计时间维度9030
• • 共14 个月（2011 年 9 月至 2012 年 11 月）
• • 数据量巨大，PB级别
• • 以MDS自定义二进制数据格式存储，为MITgcm独有
• • 模型网格复杂，为lat-lon-cap (LLC) 曲线网格 ，很难在常规地图投影中可视化。
  

在数据发布之初，该数据集存储在高度安全的NASA超级计算机上，只有获得NASA资助的研究人员才能访问。

后来，NASA Ames研究中心创建数据共享网站（https://data.nas.nasa.gov/ecco/），开放了LLC4320数据。任何人都可以通过互联网访问数据。在此网站上，您可以单击下载单个大小40GB的二进制文件。除非您知道如何解码其中的内容，否则这些文件毫无用处。
( I. P2 W6 X$ e7 }* l; X# X

0 M+ k( H6 ]3 a3 A% v% z2 S0 b+ Y

xmitgcm.llcreader

xmitgcm 是一个 python 包，可以将 MITgcm 二进制 MDS 文件读入 xarray 结构。通过 dask，xmitgcm可以实现并行计算。

代码库: https://github.com/MITgcm/xmitgcm

文档: https://xmitgcm.readthedocs.io/en/latest/

博客: https://medium.com/pangeo/petaby ... portal-81e3c5e077be

为了使二进制数据方便利用，Ryan Abernathey等开发了xmitgcm的python包，其中llcreader用于读取这些二进制文件。该模块使用xarray和dask从ECCO数据门户网站在线访问数据，使模式大数据的操作变得轻而易举。

7 i7 o: j) Y0 |$ X1 s

海面温度读取示例

以海面温度读取为例，展示其基本操作。用到了如下库：

• • xmitgcm: 提供llcreader
• • xarray: 基本数据结构和操作
• • dask: 大数据并行和lazy计算
• • sholoviews: 交互式的图像展示
  

1 导入库

import xmitgcm.llcreader as llcreader
$ n+ E% G. q% p+ {0 t+ d# V%matplotlib inline
import holoviews as hv
' @  ?, A: \. q- H. a5 k9 Nfrom holoviews.operation.datashader import regrid
& X" `6 J/ x& [8 _% O0 \, }hv.extension('bokeh')
2 o" x3 e2 G' Y  v3 T4 X) t9 P5 w
2 初始化

这里我们使用LLC4320模式数据：

model = llcreader.ECCOPortalLLC4320Model()
model

根据数据分辨率和来源，llcreader 可用模块有：

• • llcreader.ECCOPortalLLC2160Model: LLC2160 accessed via ECCO data portal
• • llcreader.ECCOPortal LLC4320Model: LLC4320 accessed via ECCO data portal
• • llcreader.PleiadesLLC2160Model: LLC2160 accessed on Pleaides filesystem
• • llcreader.PleiadesLLC4320Model: LLC4320 accessed on Pleaides filesystem
• • llcreader.CRIOSPortalASTE270Model: ASTE Release 1 accessed via AWS
• • llcreader.SverdrupASTE270Model: ASTE Release 1 accessed on Sverdrup filesystem at UT Austin
  

3 海表温度参数设置
+ _& s7 l; T3 T  |+ E' O
# d2 a) Q4 L$ H4 ]ds_sst = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon')
& e" o: t! F- D, xds_sst

这里的Theta是模式中固有的海表温度名称。这一行程序执行的lazy模式，数据并没有存储在本地内存，也不会进行计算。该变量的大小接近10T。

ds_sst.nbytes / 1e12

9.257148163328

如果想查看其他变量的名称：

print(model.varnames)

['Eta', 'KPPhbl', 'oceFWflx', 'oceQnet', 'oceQsw', 'oceSflux', 'oceTAUX', 'oceTAUY', 'PhiBot', 'Salt', 'SIarea', 'SIheff', 'SIhsalt', 'SIhsnow', 'SIuice', 'SIvice', 'Theta', 'U', 'V', 'W']

比如Eta表示海面高度，U,V,W为速度。

get_dataset模块的全部参数设置为

get_dataset(varnames=None, iter_start=None, iter_stop=None, iter_step=None, iters=None, k_levels=None, k_chunksize=1, type='faces', read_grid=True, grid_vars_to_coords=True)

常见操作有：

• • ds = model.get_dataset(varnames=['Eta'])
• • ds = model.get_dataset(varnames=['Salt', 'Theta'], k_levels=[1, 10, 40])
• • ds = model.get_dataset(varnames=['Theta'], k_levels=[0], type='latlon')
  

参考：https://xmitgcm.readthedocs.io/e ... LCModel.get_dataset

; p1 Q5 ?& `" ^4 t; R) N( m$ T

4 动态交互可视化
; X% C  N7 K; H
6 X6 P" `7 U+ M+ Qdataset = hv.Dataset(ds_sst.Theta.isel(k=0).astype('f4'))
- h! d( }  B5 \5 e- E, y: F  S6 N/ yhv_im = (dataset.to(hv.Image, ['i', 'j'], dynamic=True)
                .options(cmap='Magma', width=950, height=600, colorbar=True))

7 ]- s" q4 h! J& H: |8 _# ]%output holomap='scrubber' fps=3
$ W( M0 C2 a( |3 a, O, Hregrid(hv_im, precompute=True)

上图是南非Aghulhas Rings，可以看到强大的洋流和丰富的中小尺度涡旋。下图作为对比是LLC2160的结果，和4320的分辨能力比较有一定差距。

' r4 ]% E, N) o- z+ ~( A涡度计算示例

下面展示LLC4320涡度计算步骤。

model = llcreader.ECCOPortalLLC4320Model()
% v- G/ t5 A4 [* _. Y; _- dprint(model)

# volecity
ds = model.get_dataset(varnames=['U', 'V'], k_levels=[0], type='latlon',
, b. _4 D, R2 X6 H                        iter_start=model.iter_start,
( G( Y" q+ G1 [% o* Y0 y                        iter_stop=(model.iter_start + model.iter_step),
. b$ Q# Z& G3 n! w/ d                        read_grid=True)
# q3 u( G( f& ~3 c  F- `
# Normal gridding
" d; p1 H: `" ?( v0 D6 Uimport xgcm
grid = xgcm.Grid(ds, periodic=['X'])
: ]) x3 M1 G1 B; _# E9 o) z, R8 {
# Calculate vorticity
zeta = (-grid.diff(ds.U * ds.dxC, 'Y', boundary='fill') + grid.diff(ds.V * ds.dyC, 'X'))/ds.rAz
2 U& r9 m) O8 f: ^2 h% ]  Tzeta = zeta.squeeze().rename('vorticity').reset_coords(drop=True)

" w) b3 k- U0 {+ u* J% N2 M7 S- o# load data
zeta.load()
$ w! s, E+ f. N* Y; ?
' R4 J  _# ^8 v5 s# Show
% |" J% r% |% i  ^9 H: {) S/ tdataset = hv.Dataset(zeta)
hv_im = (dataset.to(hv.Image, ['i_g', 'j_g'])
2 G' ^& Z9 k: v5 J1 `! h3 @- C                .options(cmap='RdBu_r', width=950, height=600, colorbar=True, symmetric=True))
& c% U; |8 ^2 I0 _. Z5 @3 {) y3 m
regrid(hv_im, precompute=True)

( r  W# _. a7 S扩展：云

虽然 ECCO 门户实现了数据自由访问，但它的带宽有限，国内用户往往难以正常访问。虽然它适合交互探索，但如果想实际处理PB级别的数据，它可能无法提供足够的网络支持。

商业云存储（例如 Amazon S3 或 Google Cloud Storage）具备两全其美的优势。它既可公开访问，又具有极大的数据处理能力。

目前大量的地学大数据已经存储于云端，并可以通过Pangeo进行操作，这其中就包含LLC4320模式。

后面我将介绍云计算平台Pangeo。

2 {( [, Y5 `% o) t$ t, u  t. E
' c. n3 T; i- t* I+ E
  w  o0 R% ~/ R; H

0 e( ^4 U- }: r8 x, N' ~9 b