netCDF 网格文件：使用 2D 值从 1D 数组中提取信息

Question

我正在尝试在 Python 3 中使用地形/测深信息（基本上是一个包含 x [十进制度经度]、y [十进制度纬度] 和 z [米] 的网格）。

网格文件的扩展名为 .nc，因此是一个 netCDF 文件。通常我会在通用映射工具之类的映射工具中使用它，而不必担心 netCDF 文件的工作原理，但我需要在 Python 脚本中提取特定信息。现在这只是将数据集限制在某些经度/纬度范围内。

但是，现在我对如何获取特定 x 和 y 值的 z 信息有点迷茫。以下是我目前对数据的了解

import netCDF4

#----------------------
# Load netCDF file
#----------------------
bathymetry_file = 'C:/Users/te279/Matlab/data/gebco_08.nc'
fh = netCDF4.Dataset(bathymetry_file, mode='r')

#----------------------
# Getting information about the file
#----------------------
print(fh.file_format)

NETCDF3_CLASSIC

print(fh)

根组（NETCDF3_CLASSIC 数据模型，文件格式 NETCDF3）： 标题：GEBCO_08 网格 来源：20100927 尺寸（尺寸）：边（2），xysize（933120000） 变量（尺寸）：float64 x_range（边），float64 y_range（边），int16 z_range（边），float64间距（边），int32尺寸（边），int16 z（xysize） 组：

print(fh.dimensions.keys()) 

odict_keys(['side', 'xysize'])

print(fh.dimensions['side'])

: name = 'side', size = 2

print(fh.dimensions['xysize'])

：名称 = 'xysize'，大小 = 933120000

#----------------------
# Variables
#----------------------
print(fh.variables.keys()) # returns all available variable keys

odict_keys(['x_range', 'y_range', 'z_range', 'spacing', 'dimension', 'z'])

xrange = fh.variables['x_range'][:] 
print(xrange)

[-180。 180.] # 包含整个世界经度的值 -180 到 180

yrange = fh.variables['y_range'][:]
print(yrange) 

[-90。 90.] # 包含整个世界纬度的值 -90 到 90

zrange = fh.variables['z_range'][:]

[-10977 8685] # 包含世界的深度/地形范围

spacing = fh.variables['spacing'][:] 

[ 0.00833333 0.00833333] # x 和 y 的间距。如果乘以 x 和 y 范围，则等于维度

dimension = fh.variables['dimension'][:] 

[43200 21600] # 对应于 z 的形状，如果它是我希望的二维数组（目前是 9333120000 的一维数组 - 即 43200*21600）

z = fh.variables['z'][:] # currently an 1D array of the depth/topography/z information I want
fh.close

基于此信息，我仍然不知道如何访问 z 以获取特定的 x/y（经度/纬度）值。我想基本上我需要将 z 的一维数组转换为对应于经度/纬度值的二维数组。我只是不知道该怎么做。我在一些帖子中看到人们试图将 1D 转换为 2D 数组，但我无法知道他们从世界的哪个角落开始以及他们如何进步。

我知道有一个 3 岁的 similar post，但是，我不知道如何为我的问题找到一个模拟的“扁平数组的索引”——或者如何准确地处理它。有人可以帮忙吗？

Answer 1

您需要首先读取z 的所有三个维度（纬度、经度、深度），然后提取每个维度的值。以下是一些示例。

# Read in all 3 dimensions [lat x lon x depth] 
z = fh.variables['z'][:,:,:] 

# Topography at a single lat/lon/depth (1 value):
z_1 = z[5,5,5]

# Topography at all depths for a single lat/lon (1D array):
z_2 = z[5,5,:] 

# Topography at all latitudes and longitudes for a single depth (2D array):
z_3 = z[:,:,5] 

请注意，您为纬度/经度/深度输入的数字是该维度中的索引，而不是实际纬度。您需要事先确定要查找的值的索引。

Answer 2

我刚刚在this post 中找到了解决方案。抱歉，我以前没有看到。这是我的代码现在的样子。感谢 Dave（他在上面的帖子中回答了自己的问题）。我唯一需要做的就是尺寸必须保持整数。

import netCDF4
import numpy as np

#----------------------
# Load netCDF file
#----------------------
bathymetry_file = 'C:/Users/te279/Matlab/data/gebco_08.nc'
fh = netCDF4.Dataset(bathymetry_file, mode='r')

#----------------------
# Extract variables
#----------------------
xrange = fh.variables['x_range'][:] 
yrange = fh.variables['y_range'][:]
zz = fh.variables['z'][:] 
fh.close()

#----------------------
# Compute Lat/Lon
#----------------------
nx = (xrange[-1]-xrange[0])/spacing[0]   # num pts in x-dir
ny = (yrange[-1]-yrange[0])/spacing[1]   # num pts in y-dir
nx = nx.astype(np.integer)
ny = ny.astype(np.integer)
lon = np.linspace(xrange[0],xrange[-1],nx)
lat = np.linspace(yrange[0],yrange[-1],ny)

#----------------------
# Reshape the 1D to an 2D array
#----------------------
bathy = zz[:].reshape(ny, nx)

所以，现在当我查看 zz 和 bathy 的形状时（以下代码），前者是一个长度为 933120000 的一维数组，后者是一个维度为 43200x21600 的二维数组。

print(zz.shape)
print(bathy.shape)

下一步是使用索引来正确访问测深/地形数据，正如 N1B4 在他的帖子中描述的那样