从掩蔽的二维数组中提取平均值

Question

我想从纬度/长度/电导率网格中提取12ºx12º区域，并计算该区域的平均电导率值。我可以在lat / long网格上成功应用蒙版，但不知何故，相同的过程对导电网格不起作用。

我尝试使用for循环屏蔽，现在我正在使用numpy.ma.masked_where函数。我可以成功地绘制掩蔽结果（即：我可以看到在绘制全局地图时提取该区域），但计算出的平均电导率值对应于非掩蔽数据。

我做了一个我想做的简单例子：

x = np.linspace(1, 10, 10)
y = np.linspace(1, 10, 10)

xm = np.median(x)
ym = np.median(y)

x = ma.masked_outside(x, xm-3, xm+3)
y = ma.masked_outside(x, ym-3, ym+3)
x = np.ma.filled(x.astype(float), np.nan)
y = np.ma.filled(y.astype(float), np.nan)

x, y = np.meshgrid(x, y)

z = 2*x + 3*y

z = np.ma.masked_where(np.ma.getmask(x), z)

plt.pcolor(x, y, z)
plt.colorbar()

print('Maximum z:', np.nanmax(z))
print('Minimum z:', np.nanmin(z))
print('Mean z:', np.nanmean(z))

我的代码是：

def Observatory_Cond_Plot(filename, ndcfile, obslon, obslat, obsname, date):

files = np.array(sorted(glob.glob(filename))) #sort txt files containing the 2-D conductivitiy arrays]

filenames = ['January', 'February', 'March', 'April', 'May', 'June', 
             'July', 'August', 'September', 'October', 'November', 'December'] #used for naming output plots and files

for i, fx in zip(filenames, files):

    ndcdata = Dataset(ndcfile) #load netcdf file

    lat = ndcdata.variables['latitude'][:] #import latitude data

    long = ndcdata.variables['longitude'][:] #import longitude data

    cond = np.genfromtxt(fx)

    cond, long = shiftgrid(180., cond, long, start=False) 

    #Mask lat and long arrays and fill masks with nan values

    lat = ma.masked_outside(lat, obslat-12, obslat+12)
    long = ma.masked_outside(long, obslon-12, obslon+12)
    lat = np.ma.filled(lat.astype(float), np.nan)
    long = np.ma.filled(long.astype(float), np.nan)

    longrid, latgrid = np.meshgrid(long, lat)

    cond = np.ma.masked_where(np.ma.getmask(longrid), cond)
    cond = np.ma.filled(cond.astype(float), np.nan)

    condmean = np.nanmean(cond)

    print('Mean Conductivity is:', condmean)
    print('Minimum conductivity is:', np.nanmin(cond))
    print('Maximum conductivity is:', np.nanmax(cond))

之后，其余代码只是绘制数据

我的结果是：

平均电导率为：3.5241649673154587最小电导率为：0.497494528344129最大电导率为：5.997825822915771

但是，从tmy图中可以清楚地看出，该区域的电导率不应低于3.2 S / m。此外，打印lat，long和cond网格：

龙:

[[nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 ...
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]]

年：

[[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]
...
[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]
[nan nan nan ... nan nan nan]]

条件：

[[       nan        nan        nan ...        nan        nan        nan]
[       nan        nan        nan ...        nan        nan        nan]
 [2.86749432 2.86743283 2.86746221 ... 2.87797247 2.87265508 2.87239185]
 ...
 [       nan        nan        nan ...        nan        nan        nan]
 [       nan        nan        nan ...        nan        nan        nan]
 [       nan        nan        nan ...        nan        nan        nan]]

而且看起来面具不能正常工作。

Answer 1

问题是np.ma.filled的调用将解除long变量的掩码。 np.meshgrid也不保留面具。

您可以在创建后直接保存蒙版，也可以从蒙版创建网格网格。我相应调整了你的例子。可以看出，所有版本的numpy mean都考虑了面具。我不得不调整上限（改为2），因为平均值是相等的。

x = np.linspace(1, 10, 10)
y = np.linspace(1, 10, 10)

xm = np.median(x)
ym = np.median(y)

# Note: changed limits
x = np.ma.masked_outside(x, xm-3, xm+2)
y = np.ma.masked_outside(x, ym-3, ym+2)
xmask = np.ma.getmask(x)
ymask = np.ma.getmask(y)

x, y = np.meshgrid(x, y)
xmask, ymask = np.meshgrid(xmask, ymask)

z = 2*x + 3*y


z1 = np.ma.masked_where(np.ma.getmask(x), z)
z2 = np.ma.masked_where(xmask | ymask, z)
print(z1)
print(z2)

print('Type z1, z2:', type(z1), type(z2))
print('Maximum z1, z2:', np.nanmax(z1), np.nanmax(z2))
print('Minimum z1, z2:', np.nanmin(z1), np.nanmin(z2))
print('Mean z1, z2:', np.mean(z1), np.mean(z2) )
print('nan Mean z1, z2:', np.nanmean(z1), np.nanmean(z2) )
print('masked Mean z1, z2:', z1.mean(), z2.mean())

Answer 2

请注意，任何一种简单的平均计算（求和除以总计），如np.mean，如果你在lat-lon网格上求平均值，就不会给你正确的答案，因为当你向极点移动时，面积会发生变化。您需要采用加权平均值，按cos（lat）加权。

如你所说，你有netcdf格式的数据，我希望你允许我使用实用程序气候数据操作符（cdo）从命令行建议一个替代解决方案（在ubuntu上你可以用sudo apt install cdo安装）。

提取感兴趣的区域：

cdo sellonlatbox,lon1,lon2,lat1,lat2 infile.nc outfile.nc

然后你可以计算出正确的加权平均值

cdo fldmean infile.nc outfile.nc

你可以像这样将两者一起管道：

cdo fldmean -sellonlatbox,lon1,lon2,lat1,lat2 infile.nc outfile.nc