从NetCDF文件按坐标提取蒸散量时序数据时遇索引越界问题

问题背景 #

我尝试从NetCDF文件中提取指定坐标网格点的潜在蒸散量（Potential evapotranspiration）数据，但在获取索引值时遇到了索引越界错误，坐标数据的维度令我困惑。

NetCDF文件信息 #

NetCDF dimension information:
    Name: x
        size: 656
        type: dtype('float64')
        units: 'm'
        axis: 'X'
        long_name: 'easting of British National Grid (BNG) coordinate system'
        standard_name: 'projection_x_coordinate'
        bounds: 'x_bnds'
    Name: y
        size: 1057
        type: dtype('float64')
        units: 'm'
        axis: 'Y'
        long_name: 'northing of British National Grid (BNG) coordinate system'
        standard_name: 'projection_y_coordinate'
        bounds: 'y_bnds'
    Name: time
        size: 31
        type: dtype('float64')
        units: 'days since 1961-01-01 00:00:00 UTC'
        calendar: 'gregorian'
        axis: 'T'
        standard_name: 'time'
        long_name: 'time in days since 1961-01-01 00:00:00 UTC'
        bounds: 'time_bnds'
    Name: bnds
        size: 2
        WARNING: bnds does not contain variable attributes
NetCDF variable information:
    Name: lat
        dimensions: ('y', 'x')
        size: 693392
        type: dtype('float64')
        long_name: 'latitude'
        standard_name: 'latitude'
        units: 'degrees_north'
    Name: lon
        dimensions: ('y', 'x')
        size: 693392
        type: dtype('float64')
        long_name: 'longitude'
        standard_name: 'longitude'
        units: 'degrees_east'
    Name: time_bnds
        dimensions: ('time', 'bnds')
        size: 62
        type: dtype('float64')
    Name: x_bnds
        dimensions: ('x', 'bnds')
        size: 1312
        type: dtype('float64')
    Name: y_bnds
        dimensions: ('y', 'bnds')
        size: 2114
        type: dtype('float64')
    Name: crsOSGB
        dimensions: ()
        size: 1
        type: dtype('int32')
        grid_mapping_name: 'transverse_mercator'
        longitude_of_central_meridian: -2.0
        false_easting: 400000.0
        false_northing: -100000.0
        latitude_of_projection_origin: 49.0
        scale_factor_at_projection_origin: 0.9996012717
        longitude_of_prime_meridian: 0.0
        semi_major_axis: 6377563.396
        inverse_flattening: 299.3249646
        projected_coordinate_system_name: 'OSGB 1936 / British National Grid'
        geographic_coordinate_system_name: 'OSGB 1936'
        horizontal_datum_name: 'OSGB_1936'
        reference_ellipsoid_name: 'Airy 1830'
        prime_meridian_name: 'Greenwich'
        EPSG_code: 'EPSG:27700'
        unit: 'm'
        towgs84: array([ 375., -111.,  431.,    0.,    0.,    0.,    0.])
    Name: pet
        dimensions: ('time', 'y', 'x')
        size: 21495152
        type: dtype('float32')
        _FillValue: -99999.0
        standard_name: 'water_potential_evaporation_amount'
        standard_name_url: 'http://vocab.nerc.ac.uk/standard_name/water_potential_evaporation_amount/'
        long_name: 'Potential evapotranspiration'
        units: 'kg m-2'
        coordinates: 'lon lat'
        grid_mapping: 'crsOSGB'
        actual_range: array([0.       , 1.7461436], dtype=float32)
        cell_methods: 'time: sum'

我的代码 #

import datetime as dt
import numpy as np
from netCDF4 import Dataset  
import pandas as pd

nc_f = 'pet/chess-pe_pet_gb_1km_daily_20010101-20010131.nc'
nc_fid = Dataset(nc_f, 'r') #read file into a dataset

brora = {'name' : 'Brora', 'lon': -4.1317, 'lat' : 58.1171} #coordinates for target data

lons = nc_fid.variables['lon']
lats = nc_fid.variables['lat']
time = nc_fid.variables['time_bnds'][:]
pet = nc_fid.variables['pet'][:]

lons_idx = np.abs(lons - brora['lon']).argmin()
lats_idx = np.abs(lats - brora['lat']).argmin()
bpet = pet[:, lats_idx,lons_idx]

错误信息 #

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[96], line 4
      2 lats_idx = np.abs(lats - brora['lat']).argmin()
      3 print(lats_idx,lons_idx)
----> 4 bpet = pet[:, lats_idx,lons_idx]
      5 #bpet

File /usr/lib/python3.11/site-packages/numpy/ma/core.py:3216, in MaskedArray.__getitem__(self, indx)
   3206 """
   3207 x.__getitem__(y) <==> x[y]
   3208 
   3209 Return the item described by i, as a masked array.
   3210 
   3211 """
   3212 # We could directly use ndarray.__getitem__ on self.
   3213 # But then we would have to modify __array_finalize__ to prevent the
   3214 # mask of being reshaped if it hasn't been set up properly yet
   3215 # So it's easier to stick to the current version
-> 3216 dout = self.data[indx]
   3217 _mask = self._mask
   3219 def _is_scalar(m):

IndexError: index 600538 is out of bounds for axis 1 with size 1057

解决方案 #

问题核心是维度处理错误：lat和lon是二维数组（y×x），直接调用argmin()会返回数组扁平化后的一维索引，这个数值远大于pet数组y轴的尺寸（1057），导致越界。

需要将一维索引转换为对应的二维坐标(y_idx, x_idx)，修改后的代码如下：

import datetime as dt
import numpy as np
from netCDF4 import Dataset, num2date
import pandas as pd

nc_f = 'pet/chess-pe_pet_gb_1km_daily_20010101-20010131.nc'
nc_fid = Dataset(nc_f, 'r')

brora = {'name' : 'Brora', 'lon': -4.1317, 'lat' : 58.1171}

# 读取变量并加载为numpy数组
lons = nc_fid.variables['lon'][:]
lats = nc_fid.variables['lat'][:]
pet = nc_fid.variables['pet'][:]

# 计算最小差值的一维索引，转换为二维坐标
lon_diff = np.abs(lons - brora['lon'])
min_idx = np.unravel_index(lon_diff.argmin(), lon_diff.shape)

# 提取对应点的潜在蒸散量数据
bpet = pet[:, min_idx[0], min_idx[1]]

# 可选：将时间转换为可读格式并生成DataFrame
time_var = nc_fid.variables['time']
time_dates = num2date(time_var[:], time_var.units, time_var.calendar)
df = pd.DataFrame({'date': time_dates, 'pet': bpet})
print(df.head())

关键说明 #

• lat和lon的维度是(y, x)，np.unravel_index返回的第一个值是y轴索引，第二个是x轴索引，与pet的维度(time, y, x)完全匹配
• 读取变量时加上[:]直接加载为numpy数组，避免延迟加载带来的潜在问题
• 如果需要更精确的坐标匹配，可以使用插值方法（如scipy.interpolate.griddata），但对于规则网格数据，最近邻索引已能满足需求

内容的提问来源于stack exchange，提问作者raherin