0%

上节说到人工智能水深反演的单波束质控,这一步骤十分重要,最后的结论表明质控后的反演结果精度可提高三倍。单波束在AI反演水深的角色是label,训练数据集则是重力数据,本节介绍从卫星测高到重力异常的反演方法和数据处理。

需具备的编程知识点

  • GMT
  • Python
  • Jupyter Lab
  • Fortran(测高重力反演源码)

卫星测高数据反演重力异常

在我们的AI水深反演中,重力异常以及相关的衍生产品如垂直重力梯度、垂线偏差等是主要的训练数据集合。通常这些产品可以从Sandwell或Anderson(DTU)学术网站下载。Sandwell和DTU产品被广泛使用,其认可度高,可以拿来直接使用,因此我们无需增加自己的重力反演计算工作。

这里我们选择Sandwell的产品,因为它提供了上述所有产品。

下载地址:https://topex.ucsd.edu/pub/global_grav_1min/

如果读者希望自己重复一遍重力反演过程,也是可以的,这可以提高自己对测高重力反演的认识,我们简单概述测高重力反演:

海水和海底玄武岩的密度差异导致了微弱的重力异常信号,这种重力异常信号将海水拽向海山处,进而使海面高度产生起伏(也称大地水准面)。测高反演重力的基本原理即利用这种相关性,它的反演方法在上世纪90年代左右就有了,目前已成熟稳定。方法大致分两类,一类是基于大地水准面起伏(SSH方法),一种是基于垂线偏差(Sea surface slope,SSS方法)。SSH方法通过逆斯托克斯公式和最小二乘配置构建了重力异常和水准面起伏的关系,SSS方法通过逆Vening–Meinesz公式和Laplace方程构建了垂涎偏差和重力异常的关系,这两种方法均基于海面高度信息和重力信号的物理关系。尽管这些公式较长,难以完全理解,不过台湾交通大学Huang教授共享了一套简便的程序代码,并且在国科大授课期间分发给学员学习使用,从事本领域的读者应该容易得到这些源代码。如果您没有这些源码,也可回复本文索取,或者扫描文末二维码加入AI反演大赛微信群索取。

总而言之,重力反演并不是AI水深的重点和难点。这里我们暂时省略这部分工作,而是使用Sandwell公布的重力产品。

等一段时间,小编将邀请国内重力大咖讲解AI反演重力,据说效果可好了。

重力数据预处理

因为下载的模型是全球模型,首先使用pygmt提取研究区域重力网格数据,分别是重力异常、垂直重力梯度、东西垂线偏差、南北垂线偏差。为保持后面的边界插值准备,这里的region略大一点。为避免网格不对应,这里还特意重新网格化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
free = pygmt.grdcut(grid="../GeoModel/Free_Air_Gravity_Anomalies/Free_Air_Gravity_Anomalies.nc", region="-15.5/5.5/-4.5/4.5")
free.to_netcdf("free2.grd")

vgg = pygmt.grdcut(grid="../GeoModel/Vertical_Gravity_Gradient/Vertical_Gravity_Gradient.nc", region="-15.5/5.5/-4.5/4.5")
vgg_1m=pygmt.grd2xyz(grid=vgg)
pygmt.surface(region=reg, spacing="1m",outgrid="vgg.grd",data=vgg_1m)


nvd = pygmt.grdcut(grid="../data/north_32.1.nc", region="-15.5/5.5/-4.5/4.5")
nvd_1m = pygmt.grdsample(grid=nvd, spacing="1m")
nvd_1m=pygmt.grd2xyz(grid=nvd)
pygmt.surface(region=reg, spacing="1m",outgrid="nvd.grd",data=nvd_1m)

evd = pygmt.grdcut(grid="../data/east_32.1.nc", region="-15.5/5.5/-4.5/4.5")
evd_1m=pygmt.grd2xyz(grid=evd)
pygmt.surface(region=reg, spacing="1m",outgrid="evd.grd",data=evd_1m)

重力异常、VGG和东西、南北垂线偏差如下图,绘图代码略,这里使用的调色板是polar

短波重力计算

重力信号中实际和地形相关的是短波分量,传统的GGM方法就使用了短波重力来反演水深,短波重力可以通过布格板公式计算,这里有一个最优密度差的计算,需要通过迭代确定。通过计算,我们发现几内亚湾的最优密度差是1。

其思路是:提取实测单波束点位的重力异常,使用布格板公式计算离散点的短波重力,使用重力异常减去短波重力得到离散点的长波重力,对长波重力做网格化,使用网格化的重力异常减去长波重力,得到网格化的短波重力。

python代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
control_depth_free=pygmt.grdtrack(points=control.values,grid="free2.grd",no_skip=True)
y_train1 = control_depth_free.values[:,0:3] # lon lat depth
x_train1 = control_depth_free.values[:,3] # gravity
x_train1 = np.reshape(x_train1, (len(x_train1), 1))

x_lat = np.reshape(y_train1[:, 1], (len(y_train1[:, 1]), 1))
x_lon = np.reshape(y_train1[:, 0], (len(y_train1[:, 0]), 1))

rou=1.0
d=-8000
control_short=(y_train1[:,2]-d)*2*3.1415*6.67259*(10**-8)*rou*100000 # short wavelength gravity ano
control_short = np.reshape(control_short, (len(control_short), 1))
control_long=x_train1-control_short
short=np.concatenate((y_train1[:,0:2],control_short),axis=1)
long=np.concatenate((y_train1[:,0:2],control_long),axis=1)
reg="-15.5/5.5/-4.5/4.5"

long_bmean = pygmt.blockmean(data=long, region=reg, spacing="1m")
pygmt.surface(region=reg, spacing="1m",outgrid="long.grd",data=long_bmean)
temp=pygmt.grd2xyz(grid="free2.grd")
temp_bmean = pygmt.blockmean(data=temp, region=reg, spacing="1m")
pygmt.surface(region=reg, spacing="1m",outgrid="free3.grd",data=temp_bmean)

long0 = xr.open_dataset("long.grd")
free0 = xr.open_dataset("free3.grd")
short0 = free0-long0
short0.to_netcdf("short.grd")

下图是短波重力和海底地形的对比效果,如果不标注名称,似乎看不出来差异,这也说明短波重力和海底地形的相关性太强了。在AI中加入短波重力,预计将大幅度提供反演精度。问题是AI能否由于GGM?希望在我们即将举办的AI算法大赛中有人使用AI打败GGM。

问题:地址模型磁异常、洋壳年龄、沉积层厚度等能否应用于AI水深?

小结

下节预告

  • AI训练前的数据预处理

AI水深反演大赛交流预热群

对AI水深反演和卫星测高重力有兴趣的伙伴可以扫码入群交流。

前面的一篇卫星测高综述研究表明,与之交叉的海洋大地测量有轻微衰退迹象,在其他地学领域火热朝天的人工智能未能在大地测量、海洋测绘方面得到良好发展,是一个奇怪的现象。这其中必然有一定的历史和现实原因,在这里我们暂不做这方面的讨论。从本节开始,小编推出人工智能的海洋测量应用计划。

凡是通过经验关系得到的,AI必然也可以得到,而且由于AI的非线性特征,它还可能表现更好。我们第一个AI任务是卫星测高重力的地形反演。大家知道,通过卫星测高可以得到重力,通过重力可以得到海底地形,在地形反演环节应用了地形和重力之间的经验关系,因此AI必然也是可以在地形反演发挥作用的,但是效果能否超越最好的传统方法存疑。

接下来,我分如下几个步骤展示AI地形反演:

  • 单波束数据质控
  • 重力数据处理
  • 训练数据的多维数组构建
  • CNN、DNN网络设计
  • 预测和评估

本文先介绍数据质控,一个系列下来进展可能比较慢。热心读者朋友可以前往GAT组织下载完整的AI代码:https://github.com/GenericAltimetryTools/AIBathy

准备知识

读者在重复本系列工作之前,需要掌握一下编程技能:

  • Python编程
  • PyGMT、GMT
  • Xarray

海域选择

目前WoS可检索的AI重力地形反演仅1篇,作者是中国地质大学的Annan和wan,他们的研究区域位于非洲西海岸的几内亚湾。作者用到的深度学习网络为CNN,输入了重力异常、垂直重力梯度、垂线偏差和经纬度,标签为船测的单波束水深。作者得到了比传统方法较好的结果,那么这篇文献也鼓励我们继续开展AI领域的海洋大地测量应用。

R. F. Annan and X. Wan, “Recovering Bathymetry of the Gulf of Guinea Using Altimetry-Derived Gravity Field Products Combined via Convolutional Neural Network,” Surv Geophys, Jul. 2022, doi: 10.1007/s10712-022-09720-5.

区域绘图

我们为了进行方法和结果的比较,也选择几内亚湾进行试验。研究区域的位置示意和单波束轨迹图如下:

上图是通过GMT绘制,另外还用到了一些基本的awk的数据筛选操作:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#!/usr/bin/env bash
gmt gmtset FORMAT_GEO_MAP = dddF MAP_FRAME_WIDTH=1p
gmt gmtset FONT_ANNOT_PRIMARY 7p,Helvetica,black FONT_LABEL 7p,Helvetica,black

ps=fig1.ps

gmt pscoast -JM4.5i -R-16/10/-6/10 -Dl -A10000/0/1 -Slightblue -Glightyellow -BWSen -Bx5g5 -By2g2 --FONT_TITLE=10p -K -W0.1p --MAP_ANNOT_OBLIQUE=45 > $ps

awk 'NR%10==0 {print $1,$2,$3}' ../data/bathy.xyz | gmt psxy -R -J -Sc0.01 -Ggray -O -K >> $ps

gmt select ../data/bathy.xyz -R-15/5/-4/4 -fg > subset.d
awk 'NR%5==0 {print $1,$2,$3}' subset.d | gmt psxy -R-16/10/-6/10 -J -Sc0.03 -Gred -O -K >> $ps

gmt pscoast -Rg -JG-0/0/4c -Dc -A5000 -Gpink -Swhite -O -X3.3i -Y1.6i -K >> $ps
echo -15 -4 >tmp.d
echo -15 4 >>tmp.d
echo 5 4 >>tmp.d
echo 5 -4 >>tmp.d
echo -15 -4 >>tmp.d
gmt psxy tmp.d -W1p,red -R -J -O >> $ps

gmt psconvert $ps -A -P -Tf

数据质控

单波束数据的下载渠道为美国NGDC,为了加速下载,可以通过框选的方式选择研究区域内数据(或者回复本文索取也可)。这个区域大概又60多个航次,时间较久远,因此质量可能不佳,我们先使用最新的GEBCO模型检查这批单波束数据。

1
2
3
4
5
6
7
8
9
10
ps=fig1.ps

gmt psbasemap -R-15/5/-4/4 -JM5c -BWSen+glightblue -Bx5g5 -By2g2 -K --MAP_ANNOT_OBLIQUE=45 > $ps

gmt grdtrack ../guinea.txt -G../gebco.nc | awk '{print $1,$2,$3-$4}' > draper.txt
gmt makecpt -Crainbow -T-500/500 > colors.cpt
gmt psxy draper.txt -R -J -O -K -Sc0.05 -Ccolors.cpt >> $ps

gmt grdtrack ../guinea.txt -G../gebco.nc | awk '{print $1,$2,$3-$4}' | awk 'sqrt($3*$3)>500 {print $1,$2,$3}' >draper.txt
gmt psxy draper.txt -R -J -O -K -Sc0.1 -Gred >> $ps

通过上面的GMT代码,我们把和GEBCO差异大于500m的点全部显示出来,发现大约有466个误差较大的点。

1
2
$ gmt grdtrack ../guinea.txt -G../gebco.nc | awk '{print  $1,$2,$3,$3-$4}' | awk 'sqrt($4*$4)>500 {print  $1,$2,$3}'| gmt info
<Standard Input>: N = 466 <-14.9608/4.9988> <-3.8134/3.96975> <-5944/-80.2>

进一步可以直方图和回归分析展示单波束水深误差:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#!/usr/bin/env bash
ps=compare.ps

gmt grdtrack ../guinea.txt -G../gebco.nc | awk '{print $3,$4}' >draper.txt
gmt regress -Ey -N2 -Fxm draper.txt | head

gmt psbasemap -R-6000/-100/-6000/-100 -JX5c/5c -P -K -Bxa1000f100+l"ship borne depth/m" -Bya1000f100+l"Predicted depth/m" -BWSne > $ps
gmt regress -Ey -N2 -Fxmc -C95 draper.txt | gmt psxy -R -J -O -K -L+d -Glightorange >> $ps

gmt regress -Ey -N2 -Fxym draper.txt | awk '{printf "> error\n%s %s\n%s %s\n", $1, $2, $1, $3}' | gmt psxy -R -J -O -K -W0.1p,red,- >> $ps
gmt psxy -R -J -O -K draper.txt -Sc0.05c -Gblue >> $ps
gmt regress -Ey -N2 -Fxm draper.txt | gmt psxy -R -J -O -K -W0.5p >> $ps
awk '{print $1-$2}' draper.txt| gmt gmtmath STDIN -Sl MEAN =
awk '{print $1-$2}' draper.txt| gmt gmtmath STDIN -Sl STD =
echo -6000 -6000 "(a)"| gmt pstext -F+f8p,Helvetica,red+jBL -R -J -To -P -K -O -D0.05 -Wwhite -Gwhite >> $ps

gmt grdtrack ../guinea.txt -G../gebco.nc | awk '{print $3-$4}'| gmt pshistogram -Bxaf+l"bathymetry (meter)" -Byaf+l"Counts" -BWSne -R-400/400/0/4000 -JX1.8i/1.5i -Ggray -Z0 -T10 -O -Y6c >> $ps

gmt psconvert $ps -P -Tg -A

质控仅一行代码:

1
gmt grdtrack ../guinea.txt -G../gebco.nc | awk '{print  $1,$2,$3,$3-$4}' | awk 'sqrt($4*$4)<500 {print  $1,$2,$3}'> guinea.txtok

小结

  • AI重力地形反演文献严重偏少,一些问题仍开放待解决
  • 单波束质量有不足,需要质控
  • GMT的track功能强大,结合awk可以实现质控

下节预告

  • 重力数据处理
  • 为活跃AI在海洋测绘和海洋大地测量领域的发展,可能组织一次算法比赛,暂定:2022“安海”卫星测高水深反演人工智能算法大赛,已落实好奖金2万元,相关工作正在准备中。有兴趣的读者可以回复本文索取比赛通知(初稿)。

在地球旋转作用下,海洋表现出永恒的、自然的运动。海洋运动包含环流和波动,在时间和空间的不同尺度上具有一定的自然规律和特点,简单说中尺度海流为平衡的地转运动,可以输运物质能量,小尺度的波动是非地转非平衡运动,起到能量传播作用,不输运物质,波浪只做周期往复的椭圆运动。

当前的海洋观测手段主要是卫星观测、定点观测和漂流浮标观测,得到时间上高频(小时)和空间上高波数(亚中尺度)的全球海洋环流观测基本是不可能的,例如传统高度计分辨率的有效分辨率在100km,时间分辨率为1天。

当前波动(含内波)已经在它的领域得到了大量的独立研究,中尺度和大尺度的环流也在它的领域被大量的独立研究,而介于中尺度和小尺度之间的亚中尺度运动依旧是人类认识海洋的空白区之一。LLC4320的高分辨率允许我们使用它进行高频、高波数时空谱分析,得到亚中尺度谱形状特征,进而了解海洋动能的时空特征和注入、级串、耗散等过程。

本节,我们使用LLC4320云计算快速计算wavenumber-frequency谱(w-f谱),的得到三维的FFT结构,并分离海洋中的平衡和非平衡运动。

参考:https://github.com/pangeo-data/swot_adac_ogcms

思路

气象海洋数据往往是三维或者更多维度结构,对于含有空间和时间的三维网格数据,可以计算每个网格点的时间FFT,也可以计算每个时刻的空间FFT。分开计算的效果往往是下面两张图,这样的谱表达是时空分割的,无法同时知道时间和空间尺度的内在关系:

w-f谱本质上是在时间和空间上先后做了FFT。那么把时间和空间都同时进行FFT,则可以得到物理海洋常见的三维w-f谱,其中两个坐标轴分别表示时间频率和空间波数,颜色表示能量密度大小。

效果往往是下面的图,初看可能很混乱,我大致描述下它的含义。在时间尺度上海洋运动的周期性非常明显,在潮汐周期M2频率上会发现一条高能带,如果平行着Y轴做一条剖面线,则在M2有一个峰,f是科氏力参数,也表示地转频率或者惯性频率,和纬度有关系,它是区分平衡运动和非平衡运动的一个重要参数,图上的两条曲线是内波的不同模态,内波的第十模态(图上无此模态)往往也是区分平衡和非平衡运动的一个频率分界线。当前最新的研究成果表明,从海面观测信息中区分地转和非地转能量的参考正是惯性频率和内波模态的结合。

摘自:C. B. Rocha, S. T. Gille, T. K. Chereskin, and D. Menemenlis, “Seasonality of submesoscale dynamics in the Kuroshio Extension,” Geophysical Research Letters, vol. 43, no. 21, p. 11,304-11,311, 2016, doi: 10.1002/2016GL071349.

代码

下面使用云计算和LLC4320模式数据,介绍主要的计算步骤,由于代码略多,这里仅贴出关键部分,最后附完整代码(notebook)下载地址。

首先导入库(略)

云数据读取

同时导入LLC4320的表面流速UV,以及模式的网格文件。计算全球数据需要消耗太大的计算,这里我们选择美国西海岸的湾流区域为例进行分析,时间选择春季。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from intake import open_catalog
import xmitgcm.llcreader as llcreader
cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean/llc4320.yaml")

u = cat.LLC4320_SSU.to_dask()
v = cat.LLC4320_SSV.to_dask()

ds = xr.merge([u, v])
ds = llcreader.llcmodel.faces_dataset_to_latlon(ds, metric_vector_pairs=[])

coords = cat.LLC4320_grid.to_dask().reset_coords()
coords = llcreader.llcmodel.faces_dataset_to_latlon(coords)

llcw0=ds.sel(time=slice('2012-02-01','2012-04-30')).isel(j=slice(9555,10198-115),j_g=slice(9555,10198-115),
i=slice(15355+170,15845),i_g=slice(15355+170,15845))

载入网格信息,主要用于计算空间间隔:

1
2
3
4
5
6
7
dxC_sel = coords.dxC.isel(i_g=slice(15355+170,15845), j=slice(9555,10198-115))
dyG_sel = coords.dyG.isel(i_g=slice(15355+170,15845), j=slice(9555,10198-115))
dyC_sel = coords.dyC.isel(i=slice(15355+170,15845), j_g=slice(9555,10198-115))
dxG_sel = coords.dxG.isel(i=slice(15355+170,15845), j_g=slice(9555,10198-115))
XC_sel = coords.XC.isel(i=slice(15355+170,15845), j=slice(9555,10198-115))
YC_sel = coords.YC.isel(i=slice(15355+170,15845), j=slice(9555,10198-115))
dyG_sel= coords.dyG.isel(i_g=slice(15355+170,15845), j=slice(9555,10198-115))

合并:

1
2
llcw2 = xr.merge([llcw0,  dxC_sel, dyC_sel, dxG_sel,
dyG_sel,XC_sel, YC_sel,dyG_sel])

数据预处理

  • 网格化
1
2
3
4
5
6
7
8
9
10
11
12
13
14
metrics = {
('X',): ['dxC', 'dxG'], # X distances
('Y',): ['dyC', 'dyG'], # Y distances
# ('Z',): ['drC'], # Z distances
# ('X', 'Y'): ['rA', 'rAs', 'rAw'] # Areas
}

gridllc = Grid(llcw2, periodic=[],
coords={
# 'Z':{'center':'k','outer':'k_p1'},
'Y':{'center':'j','left':'j_g'},
'X':{'center':'i','left':'i_g'}},
metrics=metrics
)

在海洋模式中不同变量彼此偏移,并位于不同的位置,例如LLC4320关于位置有许多的坐标轴,i,i_g,j,j_g等均表示这类意思。这些网格为xarray数据模型操作带来了两难境地,这意味着它们不能使用一致的坐标来表示,xgcm 的一个目标是通过插值进行坐标统一。

  • 坐标轴规则插值

U的坐标轴原本是(j,i_g),下面将转换为(j,i):

1
gridllc.interp(llcw2.U,'X',boundary='extend') # 坐标轴的变换
  • 几何化

FFT操作的对象是距离单位,而不是网格单位,因此需要把网格近似转换为距离,并得到重组的结构化xarray,便于后续的FFT输入。

1
2
3
4
5
6
7
8
9
10
11
12
with ProgressBar():
U_llc = xr.DataArray(llcw2.U, dims=['time','YC','XC'],
coords={'time':np.arange(len(llcw2.U.time))*3600, # second unit
'YC':np.arange(0,Ny1*dy1,dy1),
'XC':np.arange(0,Nx1*dx1,dx1)} # Disctance in m
).chunk({'time':100,'YC':-1,'XC':-1})

V_llc = xr.DataArray(llcw2.U, dims=['time','YC','XC'],
coords={'time':np.arange(len(llcw2.U.time))*3600,
'YC':np.arange(0,Ny1*dy1,dy1),
'XC':np.arange(0,Nx1*dx1,dx1)}
).chunk({'time':100,'YC':-1,'XC':-1})

时空

下面我们先展示湾流地区精美的KE时空结构变化,(文末可下载全部代码)。那么怎么挖掘里面的科学信息呢?谱分析或许是很好的方法之一。

谱分析

下面是最核心的FFT计算,可以看到它的本质是连续在空间和时间上进行FFT,注意dim表示FFT处理的时空维度。得到谱信息后,还需要要做各项同性的规则化处理(代码略)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
%%time
# FFT. 先进行空间的,在进行时间的,最后得到的两个都有
Fu_llc = xrft.fft(xrft.fft(U_llc.fillna(0.),
dim=['YC','XC'], window='hann', detrend='constant',
true_phase=True, true_amplitude=True
).chunk({'time':-1,'freq_YC':200,'freq_XC':200}),
dim=['time'], window='hann', detrend='constant',
true_phase=True, true_amplitude=True
)

Fu_lls = Fu_llc.isel(freq_time=slice(len(Fu_llc.freq_time)//2,None)) * 2

Fv_llc = xrft.fft(xrft.fft(V_llc.fillna(0.),
dim=['YC','XC'], window='hann', detrend='constant',
true_phase=True, true_amplitude=True
).chunk({'time':-1,'freq_YC':200,'freq_XC':200}),
dim=['time'], window='hann', detrend='constant',
true_phase=True, true_amplitude=True
)

Fv_llc = Fv_llc.isel(freq_time=slice(len(Fv_llc.freq_time)//2,None)) * 2

绘图

首先设计一个plot函数(代码略),然后调用绘图:

1
2
3
4
5
6
7
8
9
10
enatlf = xr.DataArray(xr.apply_ufunc(gsw.f, llcw2.YC, dask='parallelized').data,
dims=['YC','XC'], coords={'YC':U_llc.YC,'XC':U_llc.XC})
omega = np.sqrt(((enatlf.isel(YC=0))**2).mean()
+ cs**2*(isoFu_llc.freq_r.isel(freq_r=slice(1,None))*2*np.pi)**2
)/2/np.pi
omega = xr.DataArray(np.minimum(omega, m2/3600), dims='freq_r',
coords={'freq_r':isoFu_llc.freq_r.isel(freq_r=slice(1,None))}
)
plot(.5*(isoFu_llc + isoFv_llc),
omega, enatlf,title=r"LLC4320")

这样我们就得到了物理海洋最常用的w-f谱分析图,同时还利用能量积分方法分离了平衡和非平衡运动。关于里面的科学性问题,已有大量论文进行了解读,本文仅技术性描述计算过程。

小结

下结预告

  • 物理海洋是不是有点难?休息休息
  • 海洋测绘的AI应用,以海底地形为例

将全球海洋模式数据存储于云端,并进行云计算是一项几乎不可能的任务。当前全球高分辨率模式数据的空间分辨率可到1/48°,时间分辨率为1小时,考虑深度分层,其数据量可以说是非常巨大,如一年LLC4320数据接近3P。

当前气象海洋云计算的先驱Pangeo已转码了ECCOv4r3低分模式全水深数据和LLC4320高分模式的表层数据,这两种模式均为全球海洋模式。近期,为推动云计算和开放共享,Pangeo物理海洋小组转码了8组海洋模式数据,并存储于云端,和Pangeo云计算无缝衔接。这其中包括FIO-COM海洋模式(First Institute of Oceanography Coupled Ocean Model)。这一批模式数据,主要集中在4个SWOT高度计定标对比区,具有丰富的亚中尺度信号,是研究海洋动力学的优秀数据集合。

关于这批数据的介绍可见:

相关论文:T. Uchida et al., “Cloud-based framework for inter-comparing submesoscale-permitting realistic ocean models,” Geosci. Model Dev., vol. 15, no. 14, pp. 5829–5856, Jul. 2022, doi: 10.5194/gmd-15-5829-2022.

Pangeo Forge

云计算依赖规范化的数据格式,为降低Pangeo社区的数据维护成本,特别提出了众包概念的Pangeo Forge。全世界任何人可以使用Pangeo Forge将开放模式产品转化为Zarr格式,形成统一规范的云数据,存储到云端,这其中的流量和计算费用已由Pangeo所获基金支付。

关于FIO-COM等模式数据的转码和存储,已经由上述论文作者完成。

读取

上述论文给了代码库地址,读者可以在GitHub找到该库,并fork到Pangeo。

8套海洋模式数据的元数据已经被归档,注意文件夹中有一个catalog.yaml文件,里面保存模式信息,如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
FIO-COM32:
description: FIO-COM32
parameters: # User parameters
region:
description: region
type: str
default: "1"
allowed: ["1", "2", "3", "4"]
datatype:
description: '"surface_hourly", "surface_flux_hourly", or "interior_daily"'
type: str
default: "surface_hourly"
allowed: ["surface_hourly", "surface_flux_hourly", "interior_daily"]
season:
description: Feb, Mar, Apr ("fma") or Aug, Sep, Oct ("aso")
type: str
default: "fma"
allowed: ["fma", "aso"]
driver: zarr
args:
urlpath: 's3://Pangeo/pangeo-forge/swot_adac/FIO-COM32/Region0{{ region }}/{{ datatype }}/{{ season }}.zarr'
consolidated: True
storage_options:
anon: True
client_kwargs:
endpoint_url: 'https://ncsa.osn.xsede.org'

这些信息给出了云端的FIO模式数据描述,比如共有四个区域(参考论文图片1),数据集中在春秋季节,数据时间采样主要为一小时。

因为数据稍微复杂,作者编写了validate_catalog模块,可以实现快速读取:

1
2
3
from validate_catalog import all_params
params_dict, cat = all_params()
params_dict.keys()

有一种模式由于版权问题,未存储云端。

dict_keys([‘GIGATL’, ‘HYCOM25’, ‘HYCOM50’, ‘eNATL60’, ‘FESOM’, ‘ORCA36’, ‘FIO-COM32’])

FIO

读取FIO模式:

1
2
3
item = "FIO-COM32"
params = params_dict[item][0]
print(item, params)

这样我们选择了湾流区域的春季数据:

FIO-COM32 {‘region’: ‘1’, ‘datatype’: ‘surface_hourly’, ‘season’: ‘fma’}

将数据通过lazy模式读入,这批数据时间是2018-2-1至2018-5-1,参数包含海表温盐流和SSH:

1
2
3
%%time
print(item, params)
ds = cat[item](**params).to_dask()

载入一个时间点的温度,中尺度现象清晰可见:

1
ds.surface_temp.sel(time='2018-05-01T10:00:00.000000000').plot(size=10)

流场:

1
ds.usurf.sel(time='2018-05-01T10:00:00.000000000').plot(size=10)

KE动能

下面算一算春季这块海域的平均动能,云计算耗时仅2s:

1
2
3
%%time
EKE = 0.5*(ds.usurf**2 + ds.usurf**2).mean('time')
EKE.load()

计算平均海面高度,耗时1s:

1
2
3
%%time
meanSSH = ds.eta_t.mean('time')
meanSSH.load()

绘图:

1
2
3
4
5
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(8,5))
np.log10(EKE).plot.contourf(vmin=-4.5, vmax=1.5, levels=30,cmap='RdBu_r')
meanSSH.plot.contour(levels=30, linewidths=0.75)

锋面

通过温度场,分析这里的锋面(fronts)。主要使用滤波技巧,首先滑动平均,然后减去这个空间平均场,得到亚中尺度的温度变化信息:

1
2
3
4
lwindow = 40
sst_ls = (ds.surface_temp.rolling(yt_ocean=lwindow, center=True).mean()
).rolling(xt_ocean=lwindow, center=True).mean();
sst_ss = ds.surface_temp - sst_ls

绘图对比,可见湾流特征清晰,锋面信息丰富,里面蕴含的科学还等待您的挖掘:

1
2
3
4
5
6
7
8
9
10
11
12
13
plt.figure(figsize=(14,4))
plt.subplot(131)
ds.surface_temp.isel(time=0).plot(vmin=-1, vmax=30, cmap='RdBu_r')
plt.title('Temp Field')

plt.subplot(132)
sst_ls.isel(time=0).plot(vmin=-1, vmax=30, cmap='RdBu_r')
plt.title('Large Scale')

plt.subplot(133)
sst_ss.isel(time=0).plot(vmin=-2, vmax=2, cmap='RdBu_r')
plt.title('Small Scale')
plt.tight_layout()

小结

  • 介绍了FIO海洋模式的读取和云计算
  • 云端计算还是一如既往的方便快捷
  • 介绍了一点物理海洋的知识:KE和锋面

下节重磅推出

  • 云端运行ROMS

基于Pangeo的云生态系统,介绍卫星测高沿轨和网格两种产品的云计算,主要包含数据读取、绘图、时空特征分析等,并略谈并行计算。

沿轨数据

沿轨测高数据的源为Copernicus的L3级别产品,目前云端数据包含了全世界主要的测高任务,其中也有我国HY-2A:

1
'al', 'alg', 'c2', 'e1', 'e1g', 'e2', 'en', 'enn', 'g2', 'h2', 'j1', 'j1g', 'j1n', 'j2', 'j2g', 'j2n', 'j3', 's3a', 's3b', 'tp', 'tpn'

导入库

1
2
3
4
5
6
7
import fsspec
import xarray as xr
import numpy as np
import hvplot
import hvplot.dask
import hvplot.pandas
import hvplot.xarray

HY-2A数据

Pangeo云端的测高数据已经转化为Zarr格式,并通过yaml文件进行元数据的规范描述,例如HY-2的元信息为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
h2:
args:
consolidated: true
storage_options:
requester_pays: True
urlpath: gs://pangeo-cnes/alti/h2
description: Haiyang-2A Global Ocean Along track CMEMS Sea Surface Height L3 product
driver: zarr
metadata:
tags:
- altimetry
- ocean
- observations
url: http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=SEALEVEL_GLO_PHY_L3_REP_OBSERVATIONS_008_062

可使用intake读取卫星测高数据,使用to_dask把数据转为xarray格式,方便运算操作:

1
2
3
from intake import open_catalog
cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean/altimetry.yaml")
ds = cat['h2'].to_dask()

时间序列

选择经纬度和SLA参数,舍弃其他参数,绘图查看。下图显示了HY-2A的100个测量点SLA,如果显示全部的轨迹点时间序列,修改slice内容包含全部时间段即可。:

1
2
ds_ll = ds[['latitude', 'longitude', 'sla_filtered']]
ds_ll.isel(time=slice(0,100,1)).sla_filtered.plot()

接下来计算日均值时间序列,这里使用了一个滑动平均函数rolling,计算量不算太小,但云计算表现不错,用时仅27s:

1
2
tmp=ds_ll.sla_filtered.isel(time=slice(0,28571752,100)).dropna("time").rolling(time=86400, center=True).mean()
tmp.plot()

空间分布

把所有数据载入,并简化coordinate:

1
ds_ll = ds[['latitude', 'longitude', 'sla_filtered']].reset_coords().astype('f4').load()

这样原本在xarray数据结构中是coordinate的经纬度就变成了variable。我们选择14天的HY-2A轨迹数据,展示空间分布图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
tmp=ds_ll.isel(time=slice(86400*14))

import matplotlib.pyplot as plt
tmp.plot.scatter( y="latitude",
x="longitude",
hue="sla_filtered",
s=2,
vmin=-0.4,
vmax=0.4,
levels=9,
cmap="jet",
aspect=2.5,
size=9, )

plt.title(f"sla from HY-2A")
plt.xlim( 0., 360.)
plt.ylim(-67., 67.)
plt.show()

空间变化

因为这里使用的不是网格文件,是轨迹文件,所以不能使用xarray的std进行网格点的变化计算,这里使用一个统计库xhistogram计算器sla的变化,也就是标准差:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from xhistogram.xarray import histogram

lon_bins = np.arange(0, 361, 2)
lat_bins = np.arange(-70, 71, 2)

# helps with memory management
ds_ll_chunked = ds_ll.chunk({'time': '5MB'})

sla_variance = histogram(ds_ll_chunked.longitude, ds_ll_chunked.latitude,
bins=[lon_bins, lat_bins],
weights=ds_ll_chunked.sla_filtered.fillna(0.)**2)

norm = histogram(ds_ll_chunked.longitude, ds_ll_chunked.latitude,
bins=[lon_bins, lat_bins])

# let's get at least 200 points in a box for it to be unmasked
thresh = 200
sla_variance = sla_variance / norm.where(norm > thresh)

在几秒钟后(云计算太快了),完成全部HY-2A数据的计算,然后绘图,我们清晰的看到南极绕极流、南非大回旋、黑潮、湾流等强流特征,这表明我们国家的HY-2A数据质量很好:

1
sla_variance.plot(y='latitude_bin', figsize=(12, 6), vmax=0.07)

网格数据

下面是对网格测高的云计算。网格测高数据的源同样是欧洲Copernicus海洋项目,通过数据格式转化,这一批数据已经以Zarr格式存储在云端,我们使用Pangeo瞬间打开这一批巨大的数据(74G):

1
2
3
from intake import open_catalog
cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean.yaml")
ds = cat["sea_surface_height"].to_dask()

数据内容包含sla和流场。数据时间从1993年到2017年,数据不能更新是Pangeo云计算的一个弊端,不仅对测高如此,对别的云数据也一样。因此Pangeo社区后来发起了Pangeo Forge项目,召集全球学者参加众包数据计划,小编非常希望未来有人专门维护一套高质量的测高云数据。

绘图使用ds.sla.isel(time=0).plot()即可,速度也是秒级,这里不展示。

时间序列

计算每一个网格文件(约9000个)的均值,得到全球海面高度日均时间序列,这一步计算相对消耗时间,大约用时1-2min:

1
sla_timeseries = ds.sla.mean(dim=('latitude', 'longitude')).load()

对数据进一步做年滑动平均,然后绘图:

1
2
3
4
5
6
sla_timeseries.plot(label='full data')
sla_timeseries.rolling(time=365, center=True).mean().plot(label='rolling annual mean')
plt.ylabel('Sea Level Anomaly [m]')
plt.title('Global Mean Sea Level')
plt.legend()
plt.grid()

空间变化

进一步,我们想研究海平面变化的经纬方向空间变化,这需要仅对一个维度进行求均值,这样绘制的图形时间沿横坐标,而SLA沿纵坐标(仿Hovmöller diagram)。

1
2
3
4
sla_hov = ds.sla.mean(dim='longitude').load()
fig, ax = plt.subplots(figsize=(12, 4))
sla_hov.name = 'Sea Level Anomaly [m]'
sla_hov.transpose().plot(vmax=0.2, ax=ax)

可以看到北纬40°有特别显著的强变化区,那里是西边界流导致的黑潮和湾流区,赤道也有一个随时间变化的流场,呈现了正弦分布。这张图中的科学太多了,小编一时无法完整解读,有兴趣的读者可以研究下。

同样我们也能研究SLA在南北向的变化时间序列。此处略。

STD空间分布

因为是网格文件的时间序列,这里可以使用xarray对每一个固定网格点的时间轴做计算,同样计算量大,也略消耗时间,耗时2min。有兴趣的话,此处可以使用dask来并行计算,时间快一倍(1min)。因为dask的启动时间和任务分配时间较长,因此综合来看,小的计算量节约时间有限,甚至还慢。大数据计算耗时在1小时之上可尝试dask并行计算。

1
2
3
4
sla_std = ds.sla.std(dim='time').load()
sla_std.name = 'Sea Level Variability [m]'
ax = sla_std.plot()
_ = plt.title('Sea Level Variability')

Dask并行计算

本次并行计算视频:

小结

  • 测高云计算比较方便,计算速度快,数据无需下载
  • Pangeo生态系统成熟,并行计算可急速大型计算

下节预告

下节将准备一些复杂的云计算:

  • FIO模式大数据处理、混合层计算
  • KE计算和谱分析

Pangeo Forge为气象海洋的云计算和大数据而生,是一种众包(crowd sourcing)科学数据制作方法,目前Pangeo社区负责该技术的开发。

云计算

云计算为科学研究提供了极大便利,使我们可以将算力部署在数据服务器周围,这些服务器往往由Google或者AMS等大厂维护,先进的硬件设施和强大算力使得数据读取速度大大提高,这避免了将数据下载到计算机进行处理。

数据格式

大气海洋领域、测绘、遥感等领域大量使用的Netcdf格式不适合大数据的云计算,为了适用于大数据和云计算,一种新的多维数据结构Zarr应运而生,成为Pangeo气象海洋云数据的规范格式。

规范化的数据是高效云计算的前提。就像Python有无数的大众开发的软件包,可以通过pip等安装,但最稳定可靠的安装是使用Conda Forge,这是一个更规范化的Python源。Pangeo Forge的出现即受其启发,主要目的是维护一个规范化的云数据源。

技术流程

Pangeo Forge的技术细节并不复杂,简单说它利用了GitHub的pull request功能,任何人可以Fork在GitHub上的Pangeo Forge仓库,对meta.yamlrecipe.py两个文件修改,设置数据的源等信息,然后pull request,这样就会在这个仓库深成一个数据制作的申请,机器人会检核申请,若无误则进入数据制作步骤。一旦制作完成,数据即在云端公开可见。

Pangeo Forge不需要用户自己下载上传数据,只要提供数据的源地址,系统将自动爬取数据,进行数据格式的转换,最终存储在Pangeo的云端(Google Cloud服务商)。

全球任何人均可以使用Pangeo Forge制作规范化云数据,并存储在Pangeo的Google云端。这其中的数据下载上传、格式转化所耗费的流量和算力,以及数据所占的云存储将产生不小的费用,均已由Pangeo所获得的基金资助支付,目前对所有人免费开放使用。

初期阶段

Pangeo Forge仍旧在初期阶段,目前完成制作的标准数据并不多,例如现在发布的数据有:

大量的申请得到许可后,进入处理列队,这一过程可能比较缓慢。目前列队中有1000多个任务正在等待处理。

Pangeo Forge实例

这里放一个Pangeo Forge的视频,视频后半段有实例讲解,演讲者即Pangeo创始人Ryan Abernathey:

尝试云数据

云计算的高效到底多震撼呢?

我们可以瞬间打开九千个NOAA的降水文件,几秒钟完成九千个文件的气候态计算,并绘图展示。
以上面视频所提到的降水为例:

1
2
3
4
5
%%time
import xarray as xr

store = 'https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr'
ds = xr.open_dataset(store, engine='zarr', chunks={})

计算降水气候态:

1
2
3
4
5
from dask.diagnostics import ProgressBar
mask = (ds.precip>=0) & (ds.precip <=1000)
with ProgressBar():
pre_clim = ds.precip.where(mask).mean('time').compute()
pre_clim.plot(size=10)

小结

  • Pangeo Forge目的是为了制作规范化的Zarr格式数据
  • 所有步骤均在云上进行,包括数据抓取,格式转换和数据部署
  • 处在项目初期,列队内容多,排队等待时间长
  • 尚不能及时的自动更新数据

Pangeo 云计算的生态系统非常值得尝试,已经较成熟,Pangeo Forge还不是很成熟,想法超前,不建议新人制作维护自己的数据。

下节预告

  • Pangeo Forge 制作FIO-COM32(海洋一所)海洋模式数据
  • 云上处理FIO模式,计算混合层厚度。

大赛简介

为加速AI和卫星测高的跨界研究,推动AI海洋大地测量领域的整体良性发展,海洋遥感与数据共享公众号携手多家海洋研究机构和高校科研单位,联合发起“安海”卫星测高水深反演算法大赛,邀请多为卫星测高和海洋大地测量专家拟定赛题,聚集国内海洋遥感和AI人才,发挥AI辅助大洋水深反演效能,解决船测水深资料不足困境。

赛道信息

本次比赛设置唯一赛道:卫星测高水深反演

背景

海底地形是最基础的海洋环境要素之一,船测单波束和多波束是主要的水深数据来源,然而2019年全球海洋实测水深只占15%,在日本财团支持的联合国海洋科学可持续发展十年项目Seabed 2030支持下,2022年这一数字已达到23.4%,卫星测高反演仍旧是海洋空白区的重要水深来源。目前的反演算法中,多采用物理模型或者经验函数,建立重力信号和地形的相关联系。地形和重力之间的关系不是简单的数学公式可以描述的,存在有复杂的非线性,而当前迅速发展的人工智能技术为该领域提供了新的途径。然而目前人工智能在海洋测绘,特别是海洋大地测量领域,仍旧少有应用。通过人工智能辅助大洋水深反演显得十分必要。

意义

人工智能的大洋水深反演不仅可挖掘重力和水深的复杂非线性关系,也可以在提升分辨率、数据融合等方面提供有价值的参考。

赛题介绍

选手使用船测单波束数据和重力异常模型、重力梯度模型、垂线偏差模型以及相关地质地球物理模型进行训练建模,使用上述模型构建网格的水深模型,并使用船测单波束检核点进行精度评估。

研究数据和验证指标

研究区域为太平洋马里亚纳海沟海域,此处地形复杂,有海沟、海山和平原,水深反演具有挑战性和代表性,前期研究成果较多,是测试算法能力的理想海域。

单波束数据来自美国NGDI数据中心,重力相关模型由Sandwell等提供,地质地球物理模型由earthbyte提供。采用随机抽样的方式,将单波束数据分为训练和测试两部分。以上数据均由大赛主办方提供。

数据验证指标包括标准差STD、偏差BIAS、相关系数和其他误差统计分布特征。

数据描述和操作

单波束数据格式为文本格式,重力等模型均为Netcdf格式,使用Python或者MATLAB均可方便的读取数据。

组委会提供了一套开源的GGM重力地质法水深反演代码,其反演精度可作为选手的参考。该程序基于MATLAB和GMT,要求选手对GMT有一定的了解:

GGM:https://github.com/GenericAltimetryTools/gravity-derived-depth

组委会还将提供Python的数据操作程序,该程序可得到训练集和测试集多维数组,选手可参考使用,快速上手。

相关数据和程序将在比赛开始后再GAT组织获得。

比赛流程

报名

大赛面向全社会开放,个人、高等院校、科研单位、企业等人员均可报名参赛。

  • 报名方式:注册Github,公众号后台发送账号或者用户名,选手将被邀请加入比赛组织GAT,https://github.com/GenericAltimetryTools
  • 选手可以单人参赛,也可以组队,每位选手只能参加一支队伍。
  • 比赛结束之前均可报名,截止时间为10月20日中午12点。

比赛规则

  • 本次比赛不提供算力和存储资源,请选手自备。
  • 比赛提供研究区域的单波束水深和基本的卫星测高重力模型,并提供基础Python代读取数据。选手只需要专注于AI算法。
  • 无特殊情况,选手不可以更改Train和Check数据。
  • 本次比赛,所有程序通过GAT组织开源。
  • 只允许单一算法模型/网络框架,不得使用多个模型投票的方式构建最终模型。
  • 比赛不设现场,均为远程参加。

评分

  • 本次比赛需要选手发送一份程序解释材料,包含核心代码解释和运行结果的误差参数。
  • 评委根据选手技术思路和结果表现进行评分。
  • 评委由算法专家和海洋测绘专家组成,评分包括算法创新性、反演误差和实际可行性组成。

论文发表

  • 本次比赛将整合优秀的AI算法,在相关SCI专刊联合发表学术论文,对得奖选手给与优先署名。
  • 选手也可自行撰写学术论文,组委会将优先推荐发表。

奖项设置

  • 冠军1名,奖金10000元人民币,颁发获奖证书,获得论文署名权。
  • 亚军1名,奖金5000元人民币,颁发获奖证书,获得论文署名权。
  • 季军1名,奖金3000元人民币,颁发获奖证书,获得论文署名权。

评委介绍

  • 海洋测绘专家、AI专家

合作伙伴

  • 指导单位:自然资源部第一海洋研究所?
  • 主办单位:海洋遥感数据共享公众号平台
  • 联合发起方:高校、研究所、
  • 顾问单位:海洋测绘学会?
  • 企业支持:青岛安海(一家海洋测量装备制造小微企业,主要研发北斗波浪浮标、潮汐浮标和漂流浮标,关注:公众号)

咨询

如在报名、比赛中有任何问题,欢迎微信扫码咨询

OPeNDAP简介

OPeNDAP(Open-source Project for a Network Data Access Protocol)可以对远程结构化数据进行变量提取。它一般在这样的场景使用:

假设这里有一些好的数据: http://test.opendap.org/dap/data/nc/sst.mnmean.nc.gz
此URL为全球海面温度的月均值,您可以做的最简单的事情是下载它。但是该URL中存在您不想要的大量数据。了解数据所包含的数据变量后,OPeNDAP允许仅提取您需要的参数。

下面介绍一个简单的例子。打开测试数据网页:

http://test.opendap.org/opendap/data/nc/sst.mnmean.nc.gz.html

勾选变量,并设置采样间隔:

注意,此时在Data URL处产生了新的数据地址:

http://test.opendap.org/opendap/data/nc/sst.mnmean.nc.gz?lat[0:1:88],lon[0:1:179],time[0:1:1856],time_bnds[0:1:1856][0:1:1],sst[0:1:1856][0:1:88][0:1:179]

可直接下载该链接对应的数据,得到sst.mnmean.nc.gz文件,也可以通过Panoply、Python等程序打开这个远程文件,默认为nc格式。

OPeNDAP获得Sentinel-6海面高异常

下面以最新的卫星高度计Sentinel-6 Michael Freilich(又称Jason-CS)为例,介绍如何远程导入SLA参数,并进行加权网格处理。

以下内容从2021年NASA的云计算Hackweek活动获得:

2021-Cloud-Hackathon:https://nasa-openscapes.github.io/2021-Cloud-Hackathon/tutorials/06_S6_OPeNDAP_Access_Gridding.html

操作本例子,需要提前掌握NASA Earthdata的授权方法:

https://nasa-openscapes.github.io/2021-Cloud-Hackathon/tutorials/04_NASA_Earthdata_Authentication.html

Sentinel-6

该卫星是欧洲哥白尼计划的一部分,其首要目标是在2020-2030年期间提供精确的全球海面高度。该卫星主要有效载荷为Poseidon-4合成孔径雷达 (SAR) 双频高度计,它采用 Ku (13.575 GHz) 和 C (5.3 GHz)波段测量距离、sigma0、有效波高和电离层校正等,Sentinel-6 同时进行低分辨率模式 (LRM) 和 SAR 高分辨率模式测量,卫星轨道参数和Jason系列相同。

https://podaac.jpl.nasa.gov/Sentinel-6

导入库

导入必要库:

os
1
2
3
4
5
6
7
import tqdm
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from concurrent.futures import ThreadPoolExecutor
from pyresample.kd_tree import resample_gauss
import pyresample as pr

一个周期的L2级数据

我们选择Sentinel-6的L2准实时数据,通过数据集的ShortName和concept-id可以唯一确定该数据。

1
2
ShortName = 'JASON_CS_S6A_L2_ALT_LR_RED_OST_NRT_F'
concept_id = 'C1968980576-POCLOUD'

轨道周期:

1
2
3
4
5
cycle = 25

url = f"https://cmr.earthdata.nasa.gov/search/granules.csv?ShortName={ShortName}&cycle={cycle}&page_size=200"

print(url)

保存文件名到文件results.csv

1
2
3
4
5
!curl --silent --output "results.csv" "$url"

files = !cat results.csv | tail --lines=+2 | cut --delimiter=',' --fields=5 | cut --delimiter='/' --fields=6

print(files.s.replace(" ", "\n"))

OPeNDAP的使用

选择第一个文件,查看路径:

1
2
3
tmp = files.l[0].split('.')[0]

print(f"https://opendap.earthdata.nasa.gov/collections/{concept_id}/granules/{tmp}.html")

选择需要的参数,一般有时间和经纬度,以及一个观测参数,这里选择海面高度异常。

1
2
3
4
5
6
7
8
9
10
11
12
variables = ['data_01_time',
'data_01_longitude',
'data_01_latitude',
'data_01_ku_ssha']

v = ",".join(variables)

urls = []
for f in files:
urls.append(f"https://opendap.earthdata.nasa.gov/collections/{concept_id}/granules/{f}4?{v}")

print(urls[0])

此时,可以看到OPeNDAP路径表示为:

https://opendap.earthdata.nasa.gov/collections/C1968980576-POCLOUD/granules/S6A_P4_2__LR_RED__NR_067_001_20220903T032524_20220903T052104_F06.nc4?data_01_time,data_01_longitude,data_01_latitude,data_01_ku_ssha

此时可以直接使用该地址下载数据。

批量下载

定义下载函数,这里使用的wget

1
2
3
4
5
6
7
8
def download(source: str):

target = os.path.basename(source.split("?")[0])

if not os.path.isfile(target):
!wget --quiet --continue --output-document $target $source

return target

下载:

1
2
3
4
5
6
7
n_workers = 12

with ThreadPoolExecutor(max_workers=n_workers) as pool:

workers = pool.map(download, urls)

files = list(tqdm.tqdm(workers, total=len(urls)))

虽然原始文件大小约 2.5MB-3.0M,但OPeNDAP子集每个文件约为100KB,因此这将大大加快数据加载速度。

Xarray批量读取nc文件

这里使用xarray的open_mfdataset,注意查看一下是否有空文件!du -sh *.nc4,有的话需要删除,否则报错:

1
2
ds = xr.open_mfdataset("./*.nc4")
print(ds)

重命名变量,去掉前缀data_01_:

1
2
3
4
5
new_variable_names = list(map(lambda x: x.split("_")[-1], variables))

map_variable_names = dict(zip(variables, new_variable_names))

ds = ds.rename(map_variable_names)

轨迹图

使用xarray的scatter函数绘轨迹点图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ds.plot.scatter( y="latitude",
x="longitude",
hue="ssha",
s=5,
vmin=-0.4,
vmax=0.4,
levels=9,
cmap="jet",
aspect=2.5,
size=9, )

plt.title(f"ssha from s6 cycle {cycle}")
plt.xlim( 0., 360.)
plt.ylim(-67., 67.)
plt.show()

网格化

卫星测高领域有很多种网格化办法,比如常用的GMT的surface。这里我们使用的网格化程序为pyresample, 并且还借鉴了ECCO V4r4的网格。

首先下载一个独立的ECCO V4r4网格文件:

1
2
3
4
5
6
7
ecco_url = "https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/ECCO_L4_GEOMETRY_05DEG_V4R4/GRID_GEOMETRY_ECCO_V4r4_latlon_0p50deg.nc"

ecco_file = download(ecco_url)

ecco_grid = xr.open_dataset(ecco_file)

print(ecco_grid)

ECCO有许多深度层,这里我们只选择0层:

1
2
ecco_grid = ecco_grid.isel(Z=0).copy()
ecco_grid

ECCO网格中有mask变量,变量的值为TrueFalse,可用于分离陆地和海洋:

1
ecco_grid.maskC.plot()

使用pyresample.geometry.SwathDefinition进行网格操作:

1
2
3
4
5
6
7
8
ecco_lons = ecco_grid.maskC.longitude.values
ecco_lats = ecco_grid.maskC.latitude.values

ecco_lons_2d, ecco_lats_2d = np.meshgrid(ecco_lons, ecco_lats)

print(ecco_lons_2d.shape, ecco_lats_2d.shape)

tgt = pr.SwathDefinition(ecco_lons_2d, ecco_lats_2d)

剔除NAN:

1
2
3
4
5
6
7
8
 nans = ~np.isnan(ds.ssha.values)

ssha = ds.ssha.values[nans]
lons = ds.longitude.values[nans]
lats = ds.latitude.values[nans]

print(lons.shape, lats.shape)
ssha.shape

之后可见数据大小小了很多,说明NaN数值比较多。(820122到513357)

经度处理,修改为-180:180:

1
2
3
lons = (lons + 180) % 360 - 180

src = pr.SwathDefinition(lons, lats)

使用kd-tree gaussian weighting neighbour方法空间插值,得到result, stddev, counts三个变量的网格信息。

1
2
3
4
5
6
7
8
9
10
11
12
result, stddev, counts = resample_gauss(
src,
ssha,
tgt,
radius_of_influence=175000,
sigmas=25000,
neighbours=100,
fill_value=np.nan,
with_uncert=True,
)

result.shape

定义函数,得到xarray格式:

1
2
3
4
5
6
def to_xrda(data):
return xr.DataArray(data,
dims=['latitude', 'longitude'],
coords={'time': time,
'latitude': ecco_lats,
'longitude': ecco_lons})

确定时间中点:

1
time = np.datetime64(ds['time'].mean().data)

SLA绘图:

1
2
3
grid = to_xrda(result)

grid.sel(latitude=slice(-67.0, 67.0)).plot(vmin=-0.4, vmax=0.4, cmap="jet", figsize=(18, 7))

STD绘图:

1
2
3
stddev = to_xrda(stddev)

stddev.sel(latitude=slice(-67.0, 67.0)).plot(robust=True, cmap="jet", figsize=(18, 7))

点个数绘图:

1
2
3
counts = to_xrda(counts)

counts.sel(latitude=slice(-67.0, 67.0)).plot(robust=True, cmap="jet", figsize=(18, 7))

小结

OPeNDAP有一定的数据提取优势,前提是知道您想要的数据在哪里。

下面夏威夷大学维护的一个海洋数据OPeNDAP网站,包含了Argo、CMIP5、ECCO、ERA5、NCEP、GFDL、HadSST、NOAA_SST、PODAAC、SODA、WOA以及MODIS、卫星高度计、盐度计等卫星测量数据。

http://apdrc.soest.hawaii.edu/dods/

预告

  • 云上处理卫星测高数据,海平面变化时空分析、KE分析
  • 云上操作LLC4320、HYCOM、FESOM、FIO-COM等海洋模式,亚中尺度分析
  • Pangeo Forge,维护高质量数据集

JULE GREGORY CHARNEY,1917年1日1日出生于加利福尼亚,1981年6月16日逝世于波士顿。

Charney的准地转涡度方程允许对大尺度大气和海洋环流进行简明的数学描述,从而实现数值天气预报。 在这和他对该领域的许多基本贡献中,Charney确定了“斜压不稳定性”,这种机制解释了中纬度天气系统的大小、结构和增长率,并且在像我们的海洋这样的旋转分层流体中是一种普遍存在的现象和气氛。他的创新研究为天气系统、水动力不稳定性、大气波传播、飓风、干旱、荒漠化、大气阻塞和洋流等理论提供了见解。

其主要成就有:

  • 1946年博士论文“Dynamics of long waves in a baroclinic westerly current”发表在Journal of Meteorology,当期目录只有这一篇文章,他强调了高层大气中的“长波”对整个大气行为的影响。
  • Jule在跟随罗斯贝做助教时,提出了QG( quasi-geostrophic)理论,可用于描述大尺度的大气和海洋运动;
  • 通过和诺依曼(John von Neumann)合作,他首次将计算机应用于天气预报,开创数值天气预报先驱(1950),Jule被认为是现代动力气象之父;
  • 1979年Charney提交气候变暖分析报告,预测了二氧化碳增加和气候变暖速度的趋势关系,其结果至今有效。

他是美国国家科学院院士,也是挪威和瑞典皇家科学院的外籍院士。他获得的众多奖项包括美国气象学会的迈辛格奖(1949年)和Rossby奖章(1964年)、航空科学研究所的洛西奖(1957年)、皇家气象学会的西蒙斯奖章(1961年)、史密森学会霍奇金斯奖章(1968年)和国际气象组织奖(1971年)。

虽然Jule是著名的气象学家,但在物理海洋领域也有他的重要地位。在MIT时,Stommel和Jule是同事,Stommel曾经想Jule寻求数学上的帮助,Jule在物理海洋领域深受Stommel的影响,并把气象领域的知识应用于海洋研究,如西边界湾流的理论研究。

下文选自美国国家科学院发表的Jule传记。

http://www.nasonline.org/publications/biographical-memoirs/memoir-pdfs/charney-jule-g.pdf

简介

Jule是二战后三十年大气科学领域的主要人物之一。气象学从一门艺术到一门科学的变化很大程度上是由于他的科学眼光和他对该领域的人和项目的彻底承诺。1946年,他与加州大学洛杉矶分校语言学学生Elinor结婚,他们有两个孩子,Elinor前一段婚姻所生的儿子也随了Charney的姓氏,他们的婚姻持续了将近21年。1967年,Jule与Lois Swirnoff结婚,Lois Swirnoff是一位画家和色彩理论家,曾是加州大学洛杉矶分校和哈佛大学的教授,他们的婚姻持续了近十年。Jule在生命最后的几年与摄影艺术家Patricia Peck一起渡过。他于1981年6月16日在波士顿去世,死于肺癌。

从数学到气象

Jule 1917年出生在旧金山,他的父母在本世纪初从白俄移民到美国。

他青年时代的大部分时间都在洛杉矶度过,14岁时他的父母闹矛盾,妈妈带着Jule搬回了纽约。Jule后来回忆说,他不喜欢纽约,但他也记得他是在一个亲戚家里偶然发现微积分书。美国的普通高中都不教微积分,接触这本书后激发了他对科学的兴趣。

他本科就读于加州大学洛杉矶校区,学习数学和物理。1939年,他跟着托马斯开始读研,1940年获得硕士学位后,很快完成了一篇论文《度量曲线空间》。

硕导托马斯有一天邀请了气象学小组J.Holmboe发表演讲,Jule在台下听讲,1941年春Holmboe邀请Jule担任他的助手,并参加了在陆军和海军赞助的ULCA气象培训计划。此时,欧洲战争和太平洋的紧张局势让大学生开始考虑从事更加实用的研究工作。为了寻求建议,Jule拜访了冯·卡门(钱学森导师),冯·卡门建议他在航空工业中从事气象学工作。1941年,Jule成为UCLA气象项目的助教和学生。

气象新星

1941年,美国只有少数大学将气象学作为一门学科。加州大学洛杉矶分校气象学小组的负责人是最近从挪威来的J.Bjerknes,他以描述冷暖锋而闻名。J. Holmboe是一位年轻的挪威人,对这些概念很熟悉,对流体动力学也比较熟悉。

另一方面,M. Neiburger 在MIT师从C-G Rossby。Rossby倾向于将流体动力学应用于大气和海洋的简化模型。在1939年,他顺着Bjerknes最近的观点,即科氏参数随纬度的变化在大尺度环流系统向东迁移中发挥了重要作用,使用了一个纯水平移动的均匀大气简单模型,得出理想化大气中从西向东移动的速度定量公式(Rossby波)。

在接下来的十年中,Jule为气象研究中带来了深刻的变化。他与J.von Neumann合作,展示了新开发的电子计算机如何通过流体动力学方程进行数值积分来预测天气。这种基于物理的程序的基本前提并不新鲜,Bjerknes在20世纪初曾提出过,L.Richardson在第一次世界大战期间也曾部分尝试过。但它已25年无进展。

Jule的大学社交生活很快乐。Jule和其他几个学生合租了一所房子,享受着热闹的社交生活。Jule后来学会了滑雪和打网球,这是他在生命的最后几年都喜欢的运动。Jule是一位能够每晚与船员一起玩扑克的科学家,能够持续赢得他们的钱,而且从不产生丝毫恶意。由于他的数学背景,Jule并不被Bjerknes和Holmboe使用的描述性推理所吸引。在Jule担任助理期间,Neiburger让他接触了Rossby的论文。在1957年Rossby去世前的十年中,Rossby和Charney交换了许多信件,在其中一篇文章中,Rossby描述了他自己的教学方法Perhaps I occasionally sought to give, or inadvertently gave, to the student a sense of battle on the intellectual battlefield. If all you do is to give them a faultless and complete and uninhabited architectural masterpiece, then you do not help them to become builders of their own.这种哲学也成为Rossby论文的特点,似乎对Charney的思想产生了永久性的影响。

大约在1944年或1945年,Charney开始认为自己有资格考虑气象学论文。他逐渐制定了自己的目标,即在中纬度大气中平均西向东气流的不稳定性理论。这是他的选择,没有老师的指导。当考虑到包含非均匀流时,大气扰动方程非常复杂,为了得到一个易于处理的数学问题,Jule发现在推导最终的微分方程时,有必要进行近似。这种推理在当时的任何科学领域中都不常见。在没有任何流体动力学家的帮助下,他完成了这项工作。

经过大量手工计算,Jule能够找到一条零增长率曲线,将短波的不稳定波与较长稳定波分开。尽管当时很少有气象学家达到这样的数学水平,但这篇论文很快被发表。后来的研究表明,完整的解更复杂,但他的解包含了最重要的方面。最重要的是,他的论文使他相信他有能力在气象学方面进行高水平的原创研究。他随后将气象学作为一项永久职业,这对大气科学的发展具有重大意义。

数值天气预报

在1946年博士论文答辩之前的几个月里,他获得了欧洲的国家研究委员会奖学金,并计划访问奥斯陆的H. Solberg(挪威的数学家)和英国剑桥的G. I. Taylor。幸运的是,Jule在途中拜访了芝加哥大学的Rossby,Rossby带领该部门进入了全盛时期。Rossby凭借其非凡的说服力,毫不费劲地说服Jule推迟他的奖学金,并在大学呆了近一年。Jule后来将今年视为他职业生涯中最具形成性的经历。

1946年8月,Rossby安排Jule参加冯·诺依曼在普林斯顿高等研究院举行的一次会议,主题是电子计算机在天气预报中的应用。包括Rossby在内的十几位美国领先的动力气象学家出席了会议。这次会议的唯一重要结果是让Jule Charney认识到,John von Neumann是一个对物理问题有着相当深刻理解的人。

1947年春天,Jule夫妇前往挪威。他们的第一站是卑尔根,这里是自第一次世界大战以来挪威前沿研究的家园。Jule很快就发现如何修改流体动力学方程,将与气象相关的大尺度运动与更快的声波和重力波分离开来,这是Richardson研究中的难题。Jule通过对动量、质量和熵的流体动力学方程中的每一项进行仔细的尺度分析,证明了该公式的合理性,由上述公式得出的预测方程组现在被称为准地转理论(QG理论),准地转理论可能是自第一次世界大战以来气象学和海洋学中最有价值的发现之一。

1948年初,冯·诺依曼邀请Jule领导他的电子计算机项目中的气象学小组,其财政支持来自海军研究办公室。Jule在研究所大院住了三年,这里的研究人员大都是该研究所的一年临时成员,主要是年轻的数学家和物理学家。对于Jule来说,这是一个令人振奋的地方。冯·诺依曼虽然经常离开研究所,但他是Jule的热心听众和积极参与者。Charney夫妇很快会见了奥本海默(一把手)夫妇和研究所的其他常任理事以及普林斯顿大学的教员。

最初,Jule采取了几个主要步骤,准备用计算机预测气流的类型。由于运算复杂导致计算机不可用,Charney将正压方程线性化,并将Rossby频率公式扩展为格林函数。试验表明计算机可以应用于全非线性涡度方程,预测结果很快发表在Rossby的《Tellus》杂志上。

Jule的创作兴趣继续迅速发展。他从20世纪30年代阅读Rossby的论文以及最近与H.Stommel和 Munk(H.Sverdrup也是Jule的论文委员会成员)的熟识中,对现有的海洋大尺度运动理论有所了解。Jule首先应用了Eliassen对Ekman理论的发展,以说明如何将海洋顶部的风应力效应用于海洋准地转内部运动的边界条件。更引人注目的一步是墨西哥湾流理论。在与Stommel进行了大量讨论后,Charney展示了在海洋内部缓慢向西移动的水团中的位涡守恒应如何导致海洋西海岸的狭窄边界流,并在流向速度上具有地转平衡。

Jule已在普林斯顿高等研究院工作了七年,他已三十八岁。1951年底,他被任命为研究院成员,任期五年,但没编制,无法长待。此时,Jule的导师和领导冯·诺依曼已在华盛顿负有重任,很明显计算机项目及其应用特色并不是该研究院的核心主题。一把手奥本海默无法向Jule承诺永久编制。当冯·诺依曼患上癌症时,情况变得更加紧急。这时,Jule意识到需要换工作了,并开始询问大学。1956年夏天,在奥本海默的敦促下,研究院向Jule提供了5000美元的特别个人补助金,Jule便去了MIT。

MIT

在他搬到MIT后,Jule可以卸下他在数值天气预报方面的大部分责任。他与数学系的W.Malkus一起,组织了地球物理流体动力学非正式研讨会。该研讨会每两周在周五下午举行,逐渐有来自哈佛大学和伍兹霍尔的气象学、海洋学、地球物理学和应用数学的人员参加,耶鲁大学、布朗大学和罗德岛大学也经常参加。他们每次都被选择在不同的地点开会,长途汽车旅行自然需要一个小时的社交时间来减压,研讨会持续了22年。

1957年11月,Jule提交了一份报告,强调了全球范围内观测大气的气象卫星、布局细节观测的飞机和雷达和分析数据的电子计算机的新作用。大约在这个时候,L.Berkner考虑建立一个国家大气科学研究中心。Jule在随后的组织会议中给与非常积极地协助。加上领导推动,导致1959年3月成立了大学大气研究联盟,并很快成立国家大气研究中心(NCAR)。

Jule外向的性格和开放的思想很快使他许多科学家和管理人员建立了友好的关系。1960年初,他被任命为总统科学咨询委员会大气科学小组成员。一年后,他与麻省理工学院物理系的B.Rossi进行了讨论,后者为J.Wiesner(肯尼迪总统的科学顾问)提供关于可能和平利用外层空间的建议。Jule安排了一次与其他几位气象学家的会议,会上建议卫星可以改善天气预报。随后,联合国第1721号和第1802号决议要求世界气象组织和国际科学联盟理事会为此制定行动实践和研究计划。1966年,他成为美国国家科学院全球大气研究计划(GARP)委员会的负责人,并一直担任这一要求很高的职位,直到1971年。

Jule夫妇在1972-73学年休假,第一部分在英国剑桥度过。在此期间,Jule对如何通过Rossby波之间的非线性相互作用产生更高频率的重力波运动进行了大量思考,但最终放弃了。

Charney一家随后留在特拉维夫的魏茨曼研究所,Jule在那里研究荒漠化理论。(植被的损失将增加地面反照率,并将更多的太阳辐射反射回太空。地面吸收的太阳辐射的减少将减少对流对空气的局部加热。这反过来将减少空气的向上运动,导致降雨量减少,并有进一步减少植被的趋势。)这是他第一次访问以色列,他将这次旅行描述为一次“令人感动的经历”。

QG理论和波数谱斜率问题

CHARNEY J G. Geostrophic Turbulence[J/OL]. Journal of the Atmospheric Sciences, 1971, 28(6): 1087-1095. DOI:10.1175/1520-0469(1971)028<1087:GT>2.0.CO;2.

1971年,Jule发表了一篇高度抽象的短文,但具有深刻的实际意义。这篇论文讨论了大气大尺度准地转运动的频谱。

苏联Kolmogorov和Obukhov多年前研究表明,三维均匀各向同性湍流的波数谱随-5/3幂次而变化。同时,一些理论家认为二维湍流是一种纯粹的数学抽象,并得出了一个更陡峭的波数为-3次方的定律。这样的气流在观察者看来比传统的风洞模式要平滑得多。Jule提出了一个令人信服的论点,即受准地转方程动力学控制的系统将具有与假设的二维情况类似的频谱,即使其运动是三维的。这一点的实际意义是,根据零散观测构建的天气图在-3定律下比在-5/3定律下更有意义

当前物理海洋领域的波数谱斜率分析,即主要来源于此,它本质上代表了能量的尺度特征。

全球变暖的首次评估

1979年,Charney主持了国家研究委员会的“二氧化碳和气候特设研究小组”。 由此产生的22页报告“二氧化碳与气候:科学评估”是最早对全球变暖的现代科学评估之一。 其主要结论可在第2页找到:“我们估计CO2翻倍导致的最可能的全球变暖接近3°C,误差为±1.5°C。”这种对气候敏感性的估计在过去40年里基本没有变化。

去世

1981年,Jule去世,享年64岁。
为了纪念他,美国气象学会设立了以他名字命名的奖励。

自1973年Skylab-3携带首颗试验性高度计卫星发射升空,已过50年。

50周年之际,国际知名遥感杂志《Remote Sensing》在线发表了中国、希腊等海洋遥感学者撰写的综述论文“Satellite Altimetry: Achievements and Future Trends by a Scientometrics Analysis”。

作者根据SCI数据库Web of Science中收录的近五十年卫星高度计领域的8000余篇科学论文,全面回顾了自上世纪70年代以来,卫星测高走过的发展历程。作者统计了卫星测高的文章年度数量、应用领域分布、研究单位的地理分布特征、卫星测高知识体系网络、作者合作关系、引文共现等,并就前沿领域和今后的发展进行了展望。


图1:测高卫星在轨历史

在过去的50年里,卫星测高研究以惊人的速度发展,从最初每年发表不到10篇文献,发展到现在每年发表超过600篇文献,特别是中国学者发表的文献已超过了美国,跃居世界第一位。


图2:测高卫星相关研究论文时间演变

研究表明卫星测高的研究领域在不断得到拓展,卫星测高早已不限于传统的海洋学和大地测量学,它已经普遍应用于陆地水文、冰冻圈、陆地地形,甚至森林和土壤监测。经过50年发展,测高精度已得到显著提高,空间/时间分辨率也不断缩小,测高研究的尺度不断向精细化发展。


图3:测高卫星研究领域

在全球卫星测高作者集群中,有很多来自我国或华人学者,他们和国内外科学家存在密切的合作,为我国卫星测高的各方面发展做出了突出贡献。


图4:测高卫星作者群体

50年间,卫星测高研究热点在逐渐的演变,从早期的雷达高度计理论研究,逐步进化到SAR高度计、GNSS-R,以及新一代宽幅高度计。趋势分析揭示了当前卫星测高的热点研究领域和新的测高技术密切相关。特别是预计将于2022年发射的SWOT任务,将使得卫星测高从“网格”过渡到“像素”分辨率阶段。研究还表明,和人工智能、深度学习等跨领域合作也是未来发展的重点方向之一。


图5:测高卫星研究的知识节点

文章链接:https://www.mdpi.com/2072-4292/14/14/3332