Python选取NetCDF变量指定时间深度范围数据及报错解决

• 现有一套日分辨率NetCDF海洋观测数据，为单个经纬度点位的观测结果，目标变量维度为('time', 'z')，其中z代表海洋深度
• 需求为筛选2019-01-01至2020-12-31时间范围、100m至1500m深度范围的变量数据并绘图，深度筛选逻辑已实现，时间筛选代码运行报错

import warnings
warnings.filterwarnings('ignore')

import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl

path_in = "./"
file_ds = xr.open_dataset(path_in + 'ocean.nc')

lon = file_ds.variables['lon'][:]
lat = file_ds.variables['lat'][:]
li = file_ds.variables['depth'][0,:,0,0] 
time = file_ds.variables['time'][:]
var = file_ds.variables['var'][:,:,0,0]

index = np.intersect1d(np.where(li<-100),np.where(li>=-1500))
tindex = np.intersect1d(np.where(time<(['2020-01-01'],dtype='datetime64[ns]'),np.where(time>=['2020-12-31'],dtype='datetime64[ns]'))

var_100_1500 = var[tindex,index] 

fig, ax = plt.subplots(figsize=(15, 9))
fillprof = ax.contourf(time[tindex],li[index],var_100_1500.T,levels=40)

clb=plt.colorbar(fillprof, orientation="vertical", pad=0.02)
clb.ax.tick_params(labelsize=18) 

tindex = np.intersect1d(np.where(time<(['2020-01-01'],dtype='datetime64[ns]'),np.where(time>=['2020-12-31'],dtype='datetime64[ns]'))
                                                               ^
SyntaxError: invalid syntax

1. 语法错误：时间索引行括号嵌套混乱，缺少np.array构造调用，参数位置错误，导致解释器无法识别语法
2. 筛选逻辑写反：原代码写的筛选条件是时间早于2020-01-01 且 时间晚于等于2020-12-31，不存在符合该条件的时间值，和需求的2019-2020范围完全相反
3. 时间值构造方式冗余：numpy可直接识别标准ISO格式的日期字符串为datetime64类型，不需要手动指定dtype

方案1：修正原有numpy索引逻辑 #

只需要修改tindex的计算逻辑，修正括号配对、调整筛选条件即可：

# 修正时间筛选逻辑，注意条件是>=起始日期，<=结束日期
start_time = np.datetime64('2019-01-01')
end_time = np.datetime64('2020-12-31')
tindex = np.intersect1d(
    np.where(time >= start_time),
    np.where(time <= end_time)
)

注意：如果时间维度带时分秒，建议将结束时间设为np.datetime64('2021-01-01')，筛选条件改为time < end_time，避免漏掉2020-12-31当天的记录

方案2：使用xarray原生筛选（推荐） #

xarray自带的.sel()方法支持直接用日期字符串做维度切片，不需要手动计算索引，代码更简洁不易出错，完整修正代码如下：

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

path_in = "./"
file_ds = xr.open_dataset(path_in + 'ocean.nc')

# 直接选取单点、时间范围、深度范围，xarray自动匹配维度
var_sel = file_ds['var'].sel(
    time=slice('2019-01-01', '2020-12-31'),
    # 若深度以海面为0、向下为负值，和原代码逻辑一致则改为 depth=slice(-1500, -100)
    depth=slice(100, 1500)
).squeeze() # squeeze去掉长度为1的经纬度维度

# 直接取筛选后的值用于绘图
time_plot = var_sel.time.values
depth_plot = var_sel.depth.values
var_plot = var_sel.values

fig, ax = plt.subplots(figsize=(15, 9))
fillprof = ax.contourf(time_plot, depth_plot, var_plot.T, levels=40)

clb=plt.colorbar(fillprof, orientation="vertical", pad=0.02)
clb.ax.tick_params(labelsize=18)
plt.show()

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