NetCDF4提取经纬度子集

Question

NetCDF4提取经纬度子集

12

我希望提取一个较大的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

- user308827

7个回答

12

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]

- Spencer Hill

2

使用这个注释会出现以下错误 IndexError: 索引不能是多维的可以通过在以lat_ind和lon_ind开头的行中添加[0]来解决此问题。 - Femkemilene

我也遇到了同样的错误，但是你在哪里添加了 [0] 以使其工作？你能提供一个例子吗？ - MrKingsley

7

如果你喜欢使用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))

- jsignell

2

ds.where仅会进行掩蔽操作，我建议改用 ds.sel(lat=slice(*lat_bnds), lon=slice(*lon_bnds))。 - user2856

当我尝试这个代码时，会在我的数据集中出现关键错误'lat不是有效的维度或坐标'，但实际上我的数据集中有'latitude'和'longitude'。当我将代码行更改为ds.sel(latitude=slice(*latbound), longitude=slice(*lonbound))时，会出现相同的错误。使用print(List.ds.keys())，可以看到纬度和经度是关键字，但仍然出现关键错误。我可以使用索引吗？ - MrKingsley

4

请注意，甚至可以在命令行上使用NCO的ncks更快地完成此操作。 ncks -v air -d latitude,40.,43. -d longitude,-89.,-96. infile.nc -O subset_infile.nc

- N1B4

2

为了模仿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无法执行我想要执行的任务时，我会使用nco。

- ClimateUnboxed

1

需要对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] ) )

- Qingling Ni

0

如果您正在使用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])

- Robert Wilson

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Favo · Accepted Answer

这很简单，您需要找到纬度和经度的上下界索引。可以通过查找最接近您要查找的值来实现。

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 ]

请注意，我假定经度维度变量为东经度，并且空气变量具有时间、纬度和经度维度。

注意，我假设经度维度变量以东经度表示，并且空气变量具有时间、纬度和经度维度。