遍历netcdf文件并运行计算-Python或R

Question

这是我第一次使用netCDF，我正在努力解决这个问题。

我有多个版本3 netcdf文件（NOAA-NARR-air.2m，全年日均）。每个文件的时间跨度在1979年至2012年之间。它们是349x 277个网格，分辨率约为32km。数据是从here下载的。

维度是时间（从1800年1月1日起的小时数），我感兴趣的变量是空气。我需要用温度来计算累积天数<；0。例如

    Day 1 = +4 degrees, accumulated days = 0
    Day 2 = -1 degrees, accumulated days = 1
    Day 3 = -2 degrees, accumulated days = 2
    Day 4 = -4 degrees, accumulated days = 3
    Day 5 = +2 degrees, accumulated days = 0
    Day 6 = -3 degrees, accumulated days = 1

我需要将这些数据存储在一个新的netcdf文件中。我对Python比较熟悉，对R也比较熟悉。什么是循环遍历每一天的最佳方式，检查前几天的值，并在此基础上，将一个值输出到一个具有完全相同维度和变量的新netcdf文件。。。。或者只需在原始netcdf文件中添加另一个变量，并输出我要查找的结果。

最好将所有文件分开还是合并？我把它们和ncrcat结合在一起，效果很好，但是文件是2.3gb。

谢谢你的意见。

我目前在python方面的进展：

import numpy
import netCDF4
#Change my working DIR
f = netCDF4.dataset('air7912.nc', 'r')
for a in f.variables:
  print(a)

#output =
     lat
     long
     x
     y
     Lambert_Conformal
     time
     time_bnds
     air

f.variables['air'][1, 1, 1]
#Output
     298.37473

为了帮助我更好地理解这一点，我使用的是哪种类型的数据结构？在上面的例子中，['air']是键吗，[1,1,1]也是键吗？得到298.37473的值。我怎样才能循环通过[1,1,1]？

Answer 1

我知道，对于2013年的这篇文章来说，这已经很晚了，但我只想指出，公认的解决方案并不能解决所提出的问题。这个问题似乎希望每个连续的温度周期的长度低于零（注意这个问题中，如果温度超过零，计数器会重置），这对于气候应用（例如农业）很重要，而公认的解决方案只给出一年中温度低于零的总天数。如果这确实是mkmitchell想要的（它已经被接受为答案），那么可以在cdo的命令行中完成，而不必担心NETCDF的输入/输出：

 cdo timsum -lec,273.15 in.nc out.nc

所以循环脚本是：

files=`ls *.nc` # pick up all the netcdf files in a directory
for file in $files ; do
    # I use 273.15 as from the question seems T is in Kelvin 
    cdo timsum -lec,273.15 $file ${file%???}_numdays.nc
done 

如果您希望得到整个期间的总数，则可以将要小得多的numdays文件分类：

cdo cat *_numdays.nc total.nc 
cdo timsum total.nc total_below_zero.nc 

但同样，这个问题似乎需要每个事件的累积天数，这是不同的，但不是由公认的答案提供的。

Answer 2

这是一个R解决方案。

infiles <- list.files("data", pattern = "nc", full.names = TRUE, include.dirs = TRUE)

outfile <- "data/air.colddays.nc"     

library(raster)

r <- raster::stack(infiles) 
r <- sum((r - 273.15) < 0)

plot(r)

Answer 3

您可以使用netCDF4中非常好的MFDataset特性将一堆文件视为一个聚合文件，而无需使用ncrcat。所以你的代码应该是这样的：

from pylab import *
import netCDF4

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.19??.nc')
# print variables
f.variables.keys()

atemp = f.variables['air']
print atemp

ntimes, ny, nx = shape(atemp)
cold_days = zeros((ny,nx),dtype=int)

for i in xrange(ntimes):
    cold_days += atemp[i,:,:].data-273.15 < 0

pcolormesh(cold_days)
colorbar()

这里有一种写文件的方法（可能有更简单的方法）：

# create NetCDF file
nco = netCDF4.Dataset('/usgs/data2/notebook/cold_days.nc','w',clobber=True)
nco.createDimension('x',nx)
nco.createDimension('y',ny)

cold_days_v = nco.createVariable('cold_days', 'i4',  ( 'y', 'x'))
cold_days_v.units='days'
cold_days_v.long_name='total number of days below 0 degC'
cold_days_v.grid_mapping = 'Lambert_Conformal'

lono = nco.createVariable('lon','f4',('y','x'))
lato = nco.createVariable('lat','f4',('y','x'))
xo = nco.createVariable('x','f4',('x'))
yo = nco.createVariable('y','f4',('y'))
lco = nco.createVariable('Lambert_Conformal','i4')

# copy all the variable attributes from original file
for var in ['lon','lat','x','y','Lambert_Conformal']:
    for att in f.variables[var].ncattrs():
        setattr(nco.variables[var],att,getattr(f.variables[var],att))

# copy variable data for lon,lat,x and y
lono[:]=f.variables['lon'][:]
lato[:]=f.variables['lat'][:]
xo[:]=f.variables['x'][:]
yo[:]=f.variables['y'][:]

#  write the cold_days data
cold_days_v[:,:]=cold_days

# copy Global attributes from original file
for att in f.ncattrs():
    setattr(nco,att,getattr(f,att))

nco.Conventions='CF-1.6'
nco.close()

如果我尝试在the Unidata NetCDF-Java Tools-UI GUI中查看结果文件，似乎可以： 还要注意，这里我刚刚下载了两个数据集进行测试，所以我使用

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.19??.nc')

作为一个例子。对于所有数据，您可以使用

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.????.nc')

或者

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.*.nc')