如何在 Python 中为每个日常文件计算 netCDF

Question

我想实现计算每日 netCDF4 文件温度最小值的事件数。我有一个如下所示的代码，但它一直告诉我索引越界。netCDf4 文件是一个具有 349x277x2920 的三维数组。第三个维度是时间，温度是每三个小时拍摄一次，所以总共有365天的试探。我想计算每日最小值，其中像素低于266，计数然后计算所有365天的每日计数一起，代码是这样的：

from osgeo import gdal
from osgeo import osr
import numpy as np
import netCDF4 as net

finalgrid=np.zeros(shape=[277,349],dtype=int)
infile=net.Dataset("W:/air.2m.2014.nc")
value=infile.variables['air'][:,:,]
infile.close()
dailymean=np.full([277,349],300,dtype=float)

for k in range(0,2921,8):
    emptygrid=np.zeros(shape=[277,349],dtype=int)
    for i in range(0,278):
        for j in range(0,350):
            try:
                if value[k][i][j]>1:
                    dailymean[i][j]=min(value[k][i][j],value[k+1][i][j],value[k+2][i][j],value[k+3][i][j],value[k+4][i][j],value[k+5][i][j],value[k+6][i][j],value[k+7][i][j])
            except:
                continue
    emptygrid[dailymean<=266.0]=emptygrid[dailymean<=266.0]+1
    finalgrid=finalgrid+emptygrid
    print k

driver=gdal.GetDriverByName('gtiff')
print finalgrid.max()
outDs=driver.Create('W:/2014aircount.tif',349,277,1,gdal.GDT_Float64)
outBand=outDs.GetRasterBand(1)
outBand.WriteArray(finalgrid,0,0)
ds=gdal.Open("W:/air.2m.2014.nc")
gt = ds.GetGeoTransform()
outDs.SetGeoTransform(gt)
outDs.SetProjection(ds.GetProjection())
#outDs.SetProjection(ds.GetProjection())
outBand.FlushCache()
ds=None

Answer 1

您的代码存在一些问题，包括读入air和air切片以计算循环部分中的每日最小温度。 下面是一个替代方法，它也只涉及时间维度上的 1 个 for 循环。

import netCDF4
import numpy as np
ncfile = netCDF4.Dataset('/nas/home/air.2m.2014.nc', 'r')
# Read in all three dimensions (lat x lon x time) below 
temp = ncfile.variables['air'][:,:,:] 
# Compute the dimension sizes
nlat, nlon, ntim = np.shape(temp)
# Compute the number of days, knowing that the time step is 3-hourly (or 8 snapshots per day)
ndays = int(ntim/8) 
# Create an empty 3D array (time x lat x lon) that will store daily min temps 
daily_min_temp = np.zeros([ndays, nlat, nlon], dtype='f4')
# Compute daily min temps across the grid
for day in range(ndays):
    daily_temp = temp[:,:,day*8:day*8+8]
    daily_min_temp[day,:,:] = np.min(daily_temp, axis=2)