python - 使用 netcdf 创建向量到数组中

Question

我对 python 相当陌生，并且发现堆栈溢出是最好的资源之一，现在我希望有人可以帮助我解决我认为是一个相当基本的问题。

我正在寻找从 netCDF 文件中提取的经纬度和降雨数据列表来创建陆地掩码。我需要从 netcdf 文件中获取数据进行排列，以便删除降雨量值为“-9999”的行。 (表示没有数据，因为它在海洋上)。我可以访问该文件，我可以创建一个网格，但是当涉及到插入降雨数据进行最终检查时，我得到了奇怪的形状，并且逻辑测试没有成功。有人可以看一下这段代码并告诉我你的想法吗？

from netCDF4 import Dataset
import numpy as np

f=Dataset('/Testing/Ensemble_grid/1970_2012_eMAST_ANUClimate_mon_evap_v1m0_197001.nc')

lat = f.variables['latitude'][:]
lon = f.variables['longitude'][:]
rainfall = np.array(f.variables['lwe_thickness_of_precipitation_amount'])
lons, lats = np.meshgrid(lon,lat)
full_ary = np.array((lats,lons))
full_lats_lons = np.swapaxes(full_ary,0,2)
rain_data = np.squeeze(rainfall,axis=(0,))
grid = np.array((full_lats_lons,rain_data))
full_grid = np.expand_dims(grid,axis=1)
full_grid_col = np.swapaxes(full_grid,0,1)
land_grid = np.logical_not(full_grid_col[:,1]==-9999.)

Answer 1

这是一种替代方法，它只需创建一个新的 2D 变量 landmask，其中每个网格单元为 0(海洋)或 1(陆地)。 (我喜欢使用 1 和 0 土地掩码，因为您可以将其转换为 bool numpy 数组，并以这种方式快速进行土地平均。)

import netCDF4
import numpy as np

ncfile = netCDF4.('/path/to/your/ncfile.nc', 'r')
lat = ncfile.variables['lat'][:]
lon = ncfile.variables['lon'][:]
# Presuming here that rainfall is 2D, if not, just read in the first time step, i.e. [0,:,:]
rain = ncfile.variables['lwe_thickness_of_precipitation_amount'][:,:] 
ncfile.close()

nlat, nlon = len(lat), len(lon)
# Populate a 2D landmask array, where 1=land and 0=ocean
landmask = np.zeros([nlat, nlon], dtype='int')
for y in range(nlat):
    for x in range(nlon):
        if rain[y,x]!=-9999: # We're at a land point
             landmask[y,x] = 1

# Now you can write out the landmask into a new netCDF file
filename_out = './landmask.nc'
ncfile_out = netCDF4.Dataset(filename_out, 'w')
ncfile_out.createDimension('lat', nlat)
ncfile_out.createDimension('lon', nlon)
lat_out = ncfile_out.createVariable('lat', 'f4', ('lat',))
lon_out = ncfile_out.createVariable('lon', 'f4', ('lon',))
landmask_out = ncfile_out.createVariable('landmask', 'i', ('lat', 'lon',))
setattr(lat_out, 'units', 'degrees_north')
setattr(lon_out, 'units', 'degrees_east')
setattr(landmask_out, 'description', '1=land 0=ocean')
lat_out[:] = lat
lon_out[:] = lon
landmask_out[:,:] = landmask[:,:]
ncfile_out.close()