从二维掩蔽阵列中提取均值

Question

我想从lat/long/conductivity网格中提取一个12ºx12º区域，并计算该区域的平均电导率值。我可以成功地在lat/long网格上应用遮罩，但不知何故相同的过程对导电网格不起作用。你知道吗

我试过用for循环屏蔽，现在我用妈妈，蒙面的在哪功能。我可以成功地绘制遮罩结果（即：绘制全局地图时可以看到区域被提取），但计算出的平均电导率值与非遮罩数据相对应。你知道吗

我举了一个简单的例子来说明我想做什么：

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.平均值如果在纬度网格上求平均值，则不会给出正确答案，因为面积会随着向两极移动而变化。你需要取一个加权平均值，用cos（lat）加权。你知道吗

正如您所说，您拥有netcdf格式的数据，我希望您允许我从命令行使用实用程序climate data operators（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

Answer 2

问题是np.ma.filled的调用将取消long变量的掩码。而且np.meshgrid不会保留掩码。你知道吗

可以在创建之后直接保存遮罩，也可以从遮罩创建网格。我相应地修改了你的例子。可以看到的是，所有版本的numpymean都考虑了掩码。我不得不调整上限（改为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())