0%

上节已备好训练数据,本节介绍CNN和DNN网络设置,跳过物理公式,直接使用备好的重力数据来预测水深,然后使用检核点水深来验证模型的精度。本节是这个AI重力水深反演系列的大结局,我们将知道以下问题的答案:

  • 加入物理机制计算的短波重力能提高AI精度吗?
  • 高维度的CNN是不是比低维的DNN好呢?
  • 人工智能是否比传统方法效果好呢?

人工智能的海底地形反演系列在这儿暂时大结局了,但这不意味着真的结束。尽管人工智能在海洋遥感、物理海洋和其他许多地学领域取得了巨大的进展,但在大地测量、海洋测绘领域的应用可能才刚开始。其背后的原因是深刻而复杂的,这是不是测绘界的传统不无关系呢?测绘的主要工作之一就是测,也就是采集真实数据,并且非常重视数据精度和真实性这个概念,《测绘法》里面对数据真实性有明确的法律要求,篡改数据可能是犯罪行为,故对AI这种虚幻的、算出来的数据是有潜意识的否定态度;若AI精度达到多波束精度,那海洋测绘从业者可能失业,整个领域将遭遇百年未有之大变局;原因可能还有更多,欢迎私信讨论。

FC-DNN

模型构建

本例使用的深度学习框架是谷歌tensorflow,如果读者仅应用深度学习,而不研究深度学习本身,使用内置的keras包即可,因为它比较友好,上手轻松。另外,使用GPU可以提高训练速度,比如2080TI的英伟达显卡训练中,DNN每个Epoch仅3s,CNN为6s,对于本案例来说,由于数据量不大,CPU应该也不会太耗费时间。

FC-DNNFully-Connected (FC) Deep Neural Network (DNN),在输入和输出层的中间有一系列全连接层,本例使用的DNN网络结构如下:

这个FC-DNN结构的中间层设计了三个全连接层,全连接顾名思义就是每个节点都互相连接,三个层的神经元数量分别是20、10、5。通过这些个连接层和神经元,经过后面的最优拟合训练,就可以构建起了重力场(input)和海底地形(output)之间的非线性关系。最终在预测水深的时候,只需要输入重力场和位置信息就可以输出对应位置的水深。

DNN结构绘图代码使用了latexhttps://github.com/GenericAltimetryTools/AIBathy/blob/main/figs/fig4/main.tex

这个FC-DNN结构对应的Python代码:

1
2
3
4
5
model = models.Sequential()
model.add(layers.Dense(20, input_shape=(None,6), activation='relu')) # change shape for case
model.add(layers.Dense(10, activation='relu'))
model.add(layers.Dense(5, activation='relu'))
model.add(layers.Dense(1))

上面的代码中input_shape=(None,6)表示输入的元素个数,需要根据具体输入参数的类型来确定,比如如果仅仅输入经纬度和重力异常,那这儿就是3activation='relu'表示激活函数为relu

Keras内置的深度学习模型确实很方便,对于非计算机方向的纯用户非常好用,仅几行代码就完成了一个模型的搭建。

设置完毕后,开始训练:

1
2
3
4
5
6
7
8
9
10
opt = optimizers.Adam(learning_rate=0.001)
Patience = 3
reduce_lr = ReduceLROnPlateau(monitor='root_mean_squared_error', factor=0.5, patience=Patience, mode='min',
min_delta=0.001, min_lr=1e-6, verbose=1)
model.compile(
optimizer=opt,
loss='mse',metrics='RootMeanSquaredError')
model.fit(train_case3, y_trian, epochs=30
, batch_size=128
,callbacks=reduce_lr)

训练过程是确定最优拟合的过程。这里面有一些常见的设置,比如optimizer优化器、loss损失函数、epochs训练的次数等等。 这些知识点可以轻易在网上获得,因为本篇幅过长,这儿就不细说了。

网格重力到网格水深的预测

Python代码:

1
2
3
grid_ai=model.predict(pre_case3, batch_size=120)
grid_depth = scalerytrain.inverse_transform(grid_ai)
depth_ai=np.concatenate((free_short.values[:,0:2], grid_depth), axis=1)

通过输入网格重力pre_case3就可以直接得到归一化的网格水深,进一步还需要还原为实际水深值。

评估

1
2
3
4
5
pre = model.predict(test_case3, batch_size=120)
pre_raw = scalerytrain.inverse_transform(pre)
m = np.abs(np.mean((y_test1-pre_raw)))
std = np.std(y_test1-pre_raw)
print(std,m)

训练完成后,使用predict函数预测检核点test_case3的水深,并将归一化的预测结果还原成真实的水深,最后就可以使用实测水深对其进行统计评估了。

或者也可以使用GMT和测试数据对导出的网格进行评估,理论上效果是一样的:

1
2
3
gmt grdtrack check_feizhou.txt -Gdepth_case1.nc | awk '{print $3,$4}' > draper.txt
awk '{print $1-$2}' draper.txt| gmt gmtmath STDIN -Sl MEAN =
awk '{print $1-$2}' draper.txt| gmt gmtmath STDIN -Sl STD =

几点结果

  • 短波重力作为输入可以明显提高模型精度,和检核点水深差异的标准偏差为73.5m;不加入短波重力,若使用重力异常、重力梯度和垂线偏差为训练数据,则STD为227m。
  • 加入短波重力之后,再加入其他数据作为输入,得到结果几乎一样,即仅使用短波重力一个量就可以得到高精度的水深预测值。
  • 使用同样的训练数据和检核数据,GGM和真值差异的STD为85m,即FC-DNN精度超过了GGM。

!注意,为了和已发表论文(Annan和Wan,2022)的结果比较,上述结果并未进行第一节所展示的数据质控,如果进行数据质控FC-DNN的精度大约在30m,几乎接近当今先进多波束的测量精度!!

CNN

CNN用卷积对图像做操作,将所使用的卷积核看做未知参数,在训练网络的过程中求出最优参数。可能会有人疑问,重力和水深不是图像,怎么用CNN呢?的确,和经典数字识别、猫狗识别完全不同,重力数据是一个网格文件,实测水深是一些船测航迹,本质上没有一连串的同样大小的图像。

使用CNN的思路大概是:每一个水深点有一个经纬度坐标,这个坐标的四周组可成一个4×4大小的二维矩阵坐标,矩阵坐标上可以安置任意的input,比如重力异常、重力梯度、垂线偏差和经纬度,这样形成6个通道(类似图像的RGB通道)的三维张量,然后每一个水深点坐标上都可以组成一个同样大小的二维矩阵,这就形成了一些列的伪图像(四维张量),从而可以满足CNN的训练。

模型架构

这里我们复现Annan和Wan的网络结构:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
model = models.Sequential()
model.add(layers.Conv2D(64,(3,3), input_shape=(4,4,6),activation='relu',padding="same", name="conv1",data_format = 'channels_last'))
model.add(layers.BatchNormalization())
model.add(layers.MaxPooling2D(pool_size=(2,2)))
model.add(layers.Conv2D(128, (3,3),activation='relu',padding="same", name="conv2",data_format = 'channels_last'))
model.add(layers.BatchNormalization())
model.add(layers.MaxPooling2D(pool_size=(2,2)))
# model.add(layers.BatchNormalization(axis=-1))
model.add(layers.Conv2D(256, (3,3),activation='relu',padding="same", name="conv3",data_format = 'channels_last'))
model.add(layers.BatchNormalization())
model.add(layers.Conv2D(512, (3,3),activation='relu',padding="same", name="conv5",data_format = 'channels_last'))
model.add(layers.BatchNormalization())
model.add(layers.Dropout(0.2))
model.add(layers.Flatten())
model.add(layers.Dense(512, activation='relu'))
model.add(layers.Dense(1))

上面有一些CNN架构的知识点,比如池化层,卷积层等等,读者可以轻松从网上查到,这里也不介绍这些内容啦。

训练、预测和评估

读者最好前往开源仓库查看完整代码,并自己运行一下,提高对这些代码的理解。

https://github.com/GenericAltimetryTools/AIBathy/blob/main/Africa_cnn.ipynb

训练

Python代码:

1
2
3
4
5
6
7
8
9
10
11
12
opt = optimizers.Adam(learning_rate=0.001)
Patience = 3
early_stop = EarlyStopping(monitor='loss', mode='min', patience=Patience)
reduce_lr = ReduceLROnPlateau(monitor='root_mean_squared_error', factor=0.5, patience=Patience, mode='min',
min_delta=0.001, min_lr=1e-7, verbose=1)
model.compile(
optimizer=opt,
loss='mse',metrics='RootMeanSquaredError')
# tf.config.experimental_run_functions_eagerly(True)
model.fit(x_train, y_train, epochs=40
, batch_size=128
,callbacks=reduce_lr)

预测

Python代码:

1
2
pre_depth = model.predict(x_pre,batch_size=120)
grid_depth = scalerdepth.inverse_transform(pre_depth)

评估

方法和FC-DNN一致。可以使用GMT对输出的水深网格进行评估,或者使用keras内置函数在训练结束后进行评估。 此外我们还可以通过差异直方图、回归分析和功率密度谱等方式对结果进行评估分析,具体代码已上传到Github的fig文件夹。

几点结果

  • CNN使用卷积,有固定窗口大小,自带滤波效果,因此预测水深略平滑。使用功率谱分析也表明在短波尺度它不如FC-DNN,但在中间尺度,它的能量高于FC-DNN。
  • CNN的精度大概是120 m,这和Annan和Wan的结论基本一致。
  • CNN效果不如GGM,并且差距较大。
  • CNN效果不如FC-DNN,这说明高维度的、更复杂的网络并不一定能得到更好的预测结果。

总结

问题答案:

  • 物理机制计算的短波重力能提高精度吗?
  • 高维度的CNN是不是比更低维的DNN好呢?
  • 人工智能是否比传统方法效果好呢?

致谢

深度学习的核心代码由张雨元同学编写,对此小编表示万分感谢。

预告

年底了,人们都更加忙碌起来。距离明年的基金时间也不远了,因此下一系列会关注自然科学基金:

  • 国家自然科学基金的海洋遥感、海洋测绘课题资助分析,总体趋势如何?热点领域在哪?
  • 卫星测高在哪些领域获得资助,趋势如何?
  • 据传,测绘口凡带着深度学习的基金会轻易被毙?是真的吗?原因是什么?

上两节分别介绍了单波束数据质控和卫星测高的重力数据处理,我们发现经过质控的单波束作为训练数据,可以明显提高水深反演的精度,重力异常信号中的短波分量和地形非常相关,加入短波重力有可能大幅提高水深反演精度。

今天介绍AI反演水深最后的前奏:张量构建和归一化处理。

张量的构建

人工智能的数据表示方式是存储在Numpy里的多维数组,数据类型是ndarray,在AI领域它还有一个更具B格的名字:张量(tensor)。一般来说,当前所有机器学习系统都使用张量作为基本数据结构,张量在AI领域是非常重要的概念,以至于Google的TensorFlow都直接以它来命名。

下表是不同维度(TensorFlow中也成为秩、阶数、度数或者是n维)的张量含义:

维度(秩) 数学实体
0 标量(只有大小)
1 向量 (有大小和方向)
2 矩阵 (由数构成的表)
3 3 维张量 (由数构成的方体)
n n 维张量 (你可以自行想象一下)

在水深反演中,我们首先也构建这样的多维张量。因为DNN和CNN的张量略有不同,这里分开介绍它们的构建方式。

DNN

假设人工智能训练数据的输入包括lat,lon,ga,vgg,evd,nvd,sg,输出只有一个od:

  • lat:latitude
  • lon:longitude
  • ga:gravity anomaly
  • vgg:vertical gravity gradient
  • evd:east-west vertical deflection
  • nvd:north-south vertical deflection
  • sg:short wavelength gravity anomaly
  • od:ocean depth

因为在DNN中,没有要求输入必须是图像,所以它接受的输入张量可以是二维矩阵。假设有10000个水深训练点,那么输入张量大小为(10000,7),输出张量为(10000,1)。注意位置信息并不单独占据两个维度,而是作为输入保存在了张量内部。

在python中使用numpy的基本函数可完成张量构建:

1
2
input =  np.concatenate((lon,lat,ga,vgg,evd,nvd,sg), axis=1)
output = np.reshape(free_short[:,2],(len(free_short[:, 2]),1))

CNN

CNN的操作对象必须是图像,因此张量必须是4维度的。CNN的张量构建稍微复杂一些,消耗的时间也长一点。整体思路是:根据经纬度确定index,然后通过index确定矩形范围,再提取重力异常等数据,增加数组维度的深度。

Index确定

使用get_indexer函数得到船测水深点的经纬度在重力网格数组中的存储位置:

1
2
con_lat_index=free_15s.indexes["lat"].get_indexer(control.values[:,1],method="nearest")
con_lon_index=free_15s.indexes["lon"].get_indexer(control.values[:,0],method="nearest")

比如我们想得到index=126处的5×5图像:

1
free_15s.isel(lat=slice(126-2,126+3,1),lon=slice(1056-2,1056+3,1))

根据这个逻辑写循环增加数据参数,假设CNN输入是5×5大小的图像:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
temp2 = []
for i in range(len(control["lat"])):
temp = free_15s.isel(lat=slice(con_lat_index[i]-2,con_lat_index[i]+3,1),lon=slice(con_lon_index[i]-2,con_lon_index[i]+3,1))
temp_vgg = vgg_15s.isel(lat=slice(con_lat_index[i]-2,con_lat_index[i]+3,1),lon=slice(con_lon_index[i]-2,con_lon_index[i]+3,1))
temp_e = east_15s.isel(lat=slice(con_lat_index[i]-2,con_lat_index[i]+3,1),lon=slice(con_lon_index[i]-2,con_lon_index[i]+3,1))
temp_n = north_15s.isel(lat=slice(con_lat_index[i]-2,con_lat_index[i]+3,1),lon=slice(con_lon_index[i]-2,con_lon_index[i]+3,1))
temp_shortg = shortg_15s.isel(lat=slice(con_lat_index[i]-2,con_lat_index[i]+3,1),lon=slice(con_lon_index[i]-2,con_lon_index[i]+3,1))
lat1 = free_15s.lat[con_lat_index[i]-2:con_lat_index[i]+3]
lon1 = free_15s.lon[con_lon_index[i]-2:con_lon_index[i]+3]
LON,LAT = np.meshgrid(lon1,lat1)
temp1 = np.concatenate([np.reshape(LON,(5,5,1)),np.reshape(LAT,(5,5,1)),np.reshape(temp.values,(5,5,1)),np.reshape(temp_vgg.values,(5,5,1)),np.reshape(temp_e.values,(5,5,1)),np.reshape(temp_n.values,(5,5,1)),np.reshape(temp_shortg.values,(5,5,1))],axis=2)
temp2.append(temp1)
temp2 = np.array(temp2)
temp2.shape

最终得到的张量空间是(10000, 5, 5, 7)

使用index来确定位置,可以加速张量构建。

归一化

归一化是将一个变量的取值缩放到[0~1]之间,也称数据标准化(Normalization)。数据标准化处理是AI的一项基础工作,不同评价指标往往具有不同的量纲和量纲单位,这样的情况会影响到数据分析的结果,为了消除指标之间的量纲影响,需要对数据进行归一化处理,解决数据指标之间的可比性问题。

优点:提高梯度下降求解最优解的速度;归一化有可能提高精度

在DNN中,我们使用了MinMaxScalertransform函数完成归一化:

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
28
# normalization
y_train2 = np.reshape(y_train1[:, 2], (len(y_train1[:, 2]), 1)) # depth

scalerlon = MinMaxScaler().fit(x_lon)
scalerlat = MinMaxScaler().fit(x_lat)
scalershort = MinMaxScaler().fit(control_short)
scalerga = MinMaxScaler().fit(control_ga)
scalervgg = MinMaxScaler().fit(control_vgg)
scalerevd = MinMaxScaler().fit(control_evd)
scalernvd = MinMaxScaler().fit(control_nvd)
scalersed = MinMaxScaler().fit(control_sed)
scalerage = MinMaxScaler().fit(control_age)
scalermoho = MinMaxScaler().fit(control_moho)
scalermag = MinMaxScaler().fit(control_mag)
scalerytrain = MinMaxScaler().fit(y_train2)

lat_train = scalerlat.transform(x_lat)
lon_train = scalerlon.transform(x_lon)
short_train = scalershort.transform(control_short)
ga_train = scalerga.transform(control_ga)
vgg_train = scalervgg.transform(control_vgg)
evd_train = scalerevd.transform(control_evd)
nvd_train = scalernvd.transform(control_nvd)
sed_train = scalersed.transform(control_sed)
age_train = scalerage.transform(control_age)
moho_train = scalermoho.transform(control_moho)
mag_train = scalermag.transform(control_mag)
y_trian = scalerytrain.transform(y_train2)

在CNN中,我们使用一样的方法完成归一化:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
m,n,k,p = np.shape(temp2)
scalerlon = MinMaxScaler().fit(np.reshape(temp2[:,:,:,4],(m*n*k,1)))
scalerlat = MinMaxScaler().fit(np.reshape(temp2[:,:,:,5],(m*n*k,1)))
scalerfree = MinMaxScaler().fit(np.reshape(temp2[:,:,:,0],(m*n*k,1)))
scalervgg = MinMaxScaler().fit(np.reshape(temp2[:,:,:,1],(m*n*k,1)))
scalernvd = MinMaxScaler().fit(np.reshape(temp2[:,:,:,2],(m*n*k,1)))
scalerevd = MinMaxScaler().fit(np.reshape(temp2[:,:,:,3],(m*n*k,1)))
# scalermae = MinMaxScaler().fit(np.reshape(temp2[:,:,:,6],(m*n*k,1)))
scalerdepth = MinMaxScaler().fit(np.reshape(control.values[:,2],(-1,1)))
lon_train = scalerlon.transform(np.reshape(temp2[:,:,:,4],(m*n*k,1)))
lat_train = scalerlat.transform(np.reshape(temp2[:,:,:,5],(m*n*k,1)))
free_train = scalerfree.transform(np.reshape(temp2[:,:,:,0],(m*n*k,1)))
vgg_train = scalervgg.transform(np.reshape(temp2[:,:,:,1],(m*n*k,1)))
nvd_train = scalernvd.transform(np.reshape(temp2[:,:,:,2],(m*n*k,1)))
evd_train = scalerevd.transform(np.reshape(temp2[:,:,:,3],(m*n*k,1)))
# mae_train = scalerevd.transform(np.reshape(temp2[:,:,:,6],(m*n*k,1)))
x_train = np.reshape(np.concatenate([free_train,vgg_train,nvd_train,evd_train,lon_train,lat_train],axis=1),(m,n,k,p))
y_train = scalerdepth.transform(np.reshape(control.values[:,2],(-1,1)))

注意对实测水深label也必须做归一化处理。

小结

  • 张量是各种维度的向量和矩阵的统称。CNN的张量维度大,比较消耗算力,DNN张量是二维矩阵,计算简单。
  • 归一化是为了解决不同量纲的数据差异问题,提高AI训练质量。

下节预告

完成全部的人工智能的重力水深反演,并展现出高精度的AI水深反演结果。

  • DNN和CNN网络设计
  • 水深预测和评估

俺最近更新很慢,希望早点上手的伙伴可查看完整代码:https://github.com/GenericAltimetryTools/AIBathy

人工智能水深反演比赛的企业总冠名权征集

为提高奖金数额,促进参赛积极性,在保留原奖金的情况下,组委现向全社会征集比赛总冠名权,请有意向的海洋遥感、海洋测绘、人工智能等领域的企业私信联系。

上节说到人工智能水深反演的单波束质控,这一步骤十分重要,最后的结论表明质控后的反演结果精度可提高三倍。单波束在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,维护高质量数据集