经纬度子集的 netcdf4 提取

Posted 2023-03-11

【中文标题】经纬度子集的 netcdf4 提取

【英文标题】：netcdf4 extract for subset of lat lon

【发布时间】：2015-05-22 01:17:01

【问题描述】：

我想提取一个相当大的 netcdf 文件的空间子集。来自Loop through netcdf files and run calculations - Python or R

from pylab import *
import netCDF4

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.1989.nc')
# print variables
f.variables.keys()
atemp = f.variables['air'] # TODO: extract spatial subset

如何仅提取与某个州（比如爱荷华州）对应的 netcdf 文件的子集。爱荷华州有以下边界纬度：

经度：89° 5' W 至 96° 31' W

纬度：40° 36' N 到 43° 30' N

【参考方案1】：

如果您在 Linux 或 macOS 中工作，则可以使用 nctoolkit (https://nctoolkit.readthedocs.io/en/latest/) 轻松处理：

import nctoolkit as nc
data = nc.open_data('/usgs/data2/rsignell/models/ncep/narr/air.2m.1989.nc')
data.crop(lon = [-(96+31/60), -(89+5/6)], lat = [40 + 36/60, 43 + 30/60])

【参考方案2】：

需要对 lonbounds 部分做些小改动（数据为东度数），因为数据中的经度值范围是 0 到 359，所以负数在这种情况下不起作用。 latli 和 latui 的计算也需要切换，因为值从北到南，从 89 到 -89。

latbounds = [ 40 , 43 ]
lonbounds = [ 260 , 270 ] # degrees east
lats = f.variables['latitude'][:] 
lons = f.variables['longitude'][:]

# latitude lower and upper index
latli = np.argmin( np.abs( lats - latbounds[1] ) )
latui = np.argmin( np.abs( lats - latbounds[0] ) ) 

# longitude lower and upper index
lonli = np.argmin( np.abs( lons - lonbounds[0] ) )
lonui = np.argmin( np.abs( lons - lonbounds[1] ) )  

【参考方案3】：

要反映 N1B4 的响应，您还可以使用气候数据操作员 (cdo) 在一行中完成：

cdo sellonlatbox,-96.5,-89,40,43 in.nc out.nc

因此要循环一组文件，我会在 BASH 脚本中执行此操作，使用 cdo 处理每个文件，然后调用您的 python 脚本：

#!/bin/bash

# pick up a list of files (I'm presuming the loop is over the years)
files=`ls /usgs/data2/rsignell/models/ncep/narr/air.2m.*.nc`

for file in $files ; do 
   # extract the location, I haven't used your exact lat/lons
   cdo sellonlatbox,-96.5,-89,40,43 $file iowa.nc

   # Call your python or R script here to process file iowa.nc
   python script
done 

我总是尝试“离线”处理我的文件，因为我发现它不太容易出错。 cdo 是 ncks 的替代品，我并不是说它更好，我只是发现它更容易记住命令。 nco 通常更强大，当 cdo 无法执行我希望执行的任务时，我会使用它。

【参考方案4】：

如果你喜欢 pandas，那么你应该考虑看看 xarray。

import xarray as xr

ds = xr.open_dataset('http://geoport.whoi.edu/thredds/dodsC/usgs/data2/rsignell/models/ncep/narr/air.2m.1980.nc',
                     decode_cf=False)
lat_bnds, lon_bnds = [40, 43], [-96, -89]
ds.sel(lat=slice(*lat_bnds), lon=slice(*lon_bnds))

【讨论】：

ds.where 将只是 mask，我建议改为 ds.sel(lat=slice(*lat_bnds), lon=slice(*lon_bnds))。

【参考方案5】：

请注意，这可以在命令行上使用NCO's ncks 更快地完成。

ncks -v air -d latitude,40.,43. -d longitude,-89.,-96. infile.nc -O subset_infile.nc

【参考方案6】：

Favo 的回答有效（我假设；尚未检查）。更直接和惯用的方法是使用 numpy 的 where 函数来查找必要的索引。

lats = f.variables['latitude'][:] 
lons = f.variables['longitude'][:]
lat_bnds, lon_bnds = [40, 43], [-96, -89]

lat_inds = np.where((lats > lat_bnds[0]) & (lats < lat_bnds[1]))
lon_inds = np.where((lons > lon_bnds[0]) & (lons < lon_bnds[1]))

air_subset = f.variables['air'][:,lat_inds,lon_inds]

【讨论】：

使用此注释，我收到以下错误 IndexError: Index cannot be multidimensional 这可以通过将 [0] 添加到以 lat_inds 和 lon_inds 开头的行来解决。

【参考方案7】：

这很容易，您必须找到纬度和经度的上下界的索引。 您可以通过查找与您要查找的值最接近的值来做到这一点。

latbounds = [ 40 , 43 ]
lonbounds = [ -96 , -89 ] # degrees east ? 
lats = f.variables['latitude'][:] 
lons = f.variables['longitude'][:]

# latitude lower and upper index
latli = np.argmin( np.abs( lats - latbounds[0] ) )
latui = np.argmin( np.abs( lats - latbounds[1] ) ) 

# longitude lower and upper index
lonli = np.argmin( np.abs( lons - lonbounds[0] ) )
lonui = np.argmin( np.abs( lons - lonbounds[1] ) )  

然后只是对变量数组进行子集化。

# Air (time, latitude, longitude) 
airSubset = f.variables['air'][ : , latli:latui , lonli:lonui ] 

注意，我假设经度维度变量以东度为单位，空气变量具有时间、纬度、经度维度。

【讨论】：

感谢收藏！这很棒

有没有一种简单的方法也可以将 airSubset 吐出为 netcdf 文件？

使用 netcdf4-python 库最直接的方法是创建一个新的 netcdf 文件，添加维度、变量名称、其属性并保存数据数组。那是大约 4~5 行代码。另一个不错的 python 库是 IRIS (scitools.org.uk/iris)，它有很好的方法来绘制、插值、重新网格和保存 netcdf 文件的子集。