5

我正在尝试从 THREDDS 提供的三角形网格模型中获取数据子集。我希望能够指定一个 LAT/LON 边界框并从该框内获取数据。数据 URL 为:

http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_MET_FORECAST.nc

使用网格数据,很容易从 THREDDS 服务器中对数据进行子集化。有谁知道获得由 THREDDS 提供服务的三角形网格的子域的最佳方法是什么?

对于网格数据,我使用 Ferret 作为我的 OPeNDAP 客户端,并且我能够编写下载过程的脚本。尽管我可以使用 Matlab、Python 或其他工具,但我想在这里做类似的事情。

谢谢,

史蒂夫

4

2 回答 2

5

Opendap 和 NetCDF 不允许使用不规则索引进行提取。您只能请求开始、停止和跨步。

并且因为这是一个三角形网格,所以不能保证同一区域中三角形的节点具有相似的索引。因此,如果您只想获取边界框中的那些节点,则必须一个接一个地请求它们。这很慢。因此,在许多情况下,确定最小和最大索引并请求整个块中的一块,然后根据需要提取索引会更快。

这是python中两种方法的示例比较。在此示例中,提取包含所有索引的子集比在每个索引上循环提取时间序列快约 10 倍:

import netCDF4
import time
import numpy as np

url='http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc'
nc = netCDF4.Dataset(url)
ncv = nc.variables
lon = ncv['lon'][:]
lat = ncv['lat'][:]
tim = ncv['time'][:]

# find indices inside box
box = [-71.4,41,-70.2,41.5]
ii = (lon>=box[0])&(lon<=box[2])&(lat>=box[1])&(lat<=box[3])
# jj will have just indices from inside the box:
jj = np.where(ii)[0]

如果我们遍历每个索引,它会很慢:

# loop over indices, extracting time series
time0=time.time()
zi = np.zeros((len(tim),len(jj)))
k=0
for j in jj:
    zi[:,k]=ncv['zeta'][:,j]
    k +=1
print('elapsed time: %d seconds' % (time.time()-time0))

elapsed time: 56 seconds

但是如果我们在每个时间步循环遍历范围,它会快得多:

time0=time.time()
zi2 = np.zeros((len(tim),len(jj)))
jmin=jj.min()
jmax=jj.max()

for i in range(len(tim)):
    ztmp = ncv['zeta'][i,jmin:jmax+1]
    zi2[i,:] = ztmp[jj-jmin]
print('elapsed time: %d seconds' % (time.time()-time0))

elapsed time: 6 seconds

当然,您的结果可能会因非结构化网格的大小、子集中点的接近程度、要提取的点数等而有所不同。但希望这能让您了解所涉及的问题。

于 2014-11-07T10:31:14.393 回答
3

如果您有非常严格的性能限制,Rich 的回答提供了最好的方法。但是,正如他所说,结果可能因网格和子集区域而异。

如果您一次只对 1 个时间步感兴趣(没有双关语),那么从 dap 服务器按原样获取整个空间域会产生边际成本和更简单的代码:

time0 = time.time()
zi3 = ncv['zeta'][0, :]
zi3 = zi3[jj]
print('elapsed time: %d seconds' % (time.time() - time0))

elapsed time: 0 seconds

在分析时(即使在预分配数组时),Rich 所做的粗略最小/最大子集的速度大约是 2 倍,给定这个网格和特定的空间子集。NECOFS 网格在非结构化网格世界中相对较小。例如,您可能会遇到大西洋规模的 ADCIRC 网格问题。

后记:无论如何,当您获得 lat 和 lon 以确定您的子集索引时,您正在抓取整个网格。

于 2014-11-07T14:55:17.257 回答