Извлечение средних значений из замаскированного двумерного массива - PullRequest
Новая
       33

Извлечение средних значений из замаскированного двумерного массива

2 голосов
/ 28 марта 2019

Я хочу извлечь область 12º x 12º из сетки широта / долгота / проводимость и вычислить средние значения проводимости в этой области. Я могу успешно применять маски к сеткам широты и долготы, но почему-то тот же процесс не работает для сетки электропроводности.

Я пытался маскировать циклы 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 См / м. А также печатные сетки 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]]

И, похоже, маска не работает должным образом.

Ответы [ 2 ]

1 голос
/ 28 марта 2019

Проблема в том, что вызов 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())
0 голосов
/ 28 марта 2019

Остерегайтесь того, что любой простой расчет среднего значения (суммирование и деление на общее число), например, np.mean, не даст вам правильного ответа, если вы будете усреднять по сетке широт, поскольку область изменяется по мере вашего движенияк полюсам.Вам нужно взять средневзвешенное значение, взвешиваясь по 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
...