0

我想实现计算每日 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
4

1 回答 1

0

您的代码存在一些问题，包括在循环部分中计算每日最低温度的读入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)
于 2016-07-13T14:43:27.743 回答