netcdf4提取经纬度子集

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

来源

2015-03-19 user308827

那么这是很容易，你必须找到纬度和经度上下界的索引。您可以通过查找与您正在查找的最接近的值来完成此操作。

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 ]

注意，我假定经度尺寸可变以度东，并且空气变量具有时间，纬度，经度的尺寸。

来源

2015-03-19 02:24:05 Favo

Thanks Favo！这太棒了 – user308827 2015-03-20 03:05:25

有没有一种直接的方式将airSubset吐出为netcdf文件呢？ – user308827 2015-03-20 03:07:30

使用netcdf4-python库最直接的方法是创建一个新的netcdf文件，添加维度，变量名称，属性以及保存数据数组。这大约是4〜5行代码。另一个好的python库是IRIS（http://scitools.org.uk/iris），它有很好的方法来绘制，插入，重新编译和保存netcdf文件的子集。 – Favo 2015-03-22 22:06:32

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]

来源

2015-03-19 14:47:04

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

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

来源

2016-02-10 16:28:26 N1B4

如果你喜欢大熊猫，那么你应该考虑检查出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))

来源

2016-05-05 19:14:37 jsignell

ds.where只是[掩码]（http://xarray.pydata.org/en/stable/indexing.html#masking-with-where），我建议'ds.sel（lat = slice（* lat_bnds ），lon = slice（* lon_bnds））''。 – Luke 2016-11-23 04:00:26

要镜像从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无法执行我希望执行的任务时，我会诉诸它。

来源

2017-01-10 10:18:33

我真的不能评论答案，因为我只是一个潜伏者。 @ user308827大多是正确的。由于数值是从1到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]))

来源

2017-05-05 20:57:27

netcdf4提取经纬度子集

回答

相关问题