使用OPeNDAP提取Sentinel-6卫星高度计数据

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

导入库

导入必要库：

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'

轨道周期：

cycle = 25

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

print(url)

保存文件名到文件results.csv：

!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")

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

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：

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

下载：

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_:

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函数绘轨迹点图：

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网格文件：

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变量，变量的值为True或False，可用于分离陆地和海洋：

1	ecco_grid.maskC.plot()

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

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：

 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三个变量的网格信息。

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格式：

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，维护高质量数据集