Python (numpy) усреднение двух- и трехмерных массивов подсетки - PullRequest
Новая
       20

Python (numpy) усреднение двух- и трехмерных массивов подсетки

1 голос
/ 11 мая 2019

Использование Python для извлечения данных атмосферы в сетке на различных уровнях и преобразования в netCDF <- COMPLETED </p>

Использование Python для поиска данных в сетке для региона, а затем усреднение этих данных по сеткам подрешеток (2x2) <- INCORRECT</p>

Я могу заставить это работать в Octave / Matlab, но я бы хотел оставить все это в Python.Я полагаю, что проблема заключается в синтаксисе индексации и моей неспособности превратить Python в представление в отношении индексации.

Данные: одномерный массив долгот, широт и уровней давления.Долгое время имеют 49 элементов, широты имеют 13 элементов, а уровни имеют 12. Данные, которые я пытаюсь усреднить, представляют собой в первую очередь 2D-матрицу (13x49), а во втором - 3D-матрицу (форма = 12x13x49).

#DEFINE LARGE AREA OF GLOBLE
londim_g  = 49
latdim_g  = 13
lonmin_g  = 60
lonmax_g  = 180
latmin_g  = -60
latmax_g  = -30
dlat=dlon = 2.5
lats_g      = arange(latmin_g,latmax_g+dlon,dlon)
lons_g      = arange(lonmin_g,lonmax_g+dlat,dlat)
LON_G,LAT_G = meshgrid(lons_g,lats_g) #THE SHAPE OF THIS IS A PROBLEM!!
# DEFINE SMALLER REGION
lonmin  = 120;
lonmax  = 130;
latmax  = -40;
latmin  = -50;
N       = 2; #THIS IS NxN SUB-GRID AVERAGE OF SMALLER REGION
ind  = argwhere( (LON_G>=lonmin) & (LON_G<=lonmax) & (LAT_G<=latmax) & (LAT_G>=latmin) )
ri   = ind[:,0]; 
ci   = ind[:,1];
LON = LON_G[ix_(ri,ci)]
LAT = LAT_G[ix_(ri,ci)]
LON = LON[1].reshape(5,5) #THIS IS STEP IS A RESULT OF LON_G,LAT_G BEING MIS-SHAPEN
LAT = LAT[1].reshape(5,5) #THIS IS STEP IS A RESULT OF LON_G,LAT_G BEING MIS-SHAPEN
# AVERAGE on NxN sub-grids such that
#INDEX GRID
# 
# Essentially we averaging each sub-grid within the domain, that is each 2x2, grid points
# IF the following is the domain:
#
#      (ln1,lt1)      (ln2,lt1)     (ln3,lt1)     (ln4,lt1)     (ln5,lt1)
#
#      (ln1,lt2)      (ln2,lt2)     (ln3,lt2)     (ln4,lt2)     (ln5,lt2)
#
#      (ln1,lt3)      (ln2,lt3)     (ln3,lt3)     (ln4,lt3)     (ln5,lt3)
#
#      (ln1,lt4)      (ln2,lt4)     (ln3,lt4)     (ln4,lt4)     (ln5,lt4)
#
#      (ln1,lt5)      (ln2,lt5)     (ln3,lt5)     (ln4,lt5)     (ln5,lt5)
#
# then the first sub-grid is:
#
#      (ln1,lt1)      (ln2,lt1) 
#
#      (ln1,lt2)      (ln2,lt2) 
#
# the next sub-grid is:
#
#      (ln2,lt1)     (ln3,lt1)
#
#      (ln2,lt2)     (ln3,lt2)
#
# So on, and so forth. If we associate each grid point with it's data then compute the average
# value of that sub-grid then we will have an `array', in this of 16 mean values, i.e.:
#
#      (ln1,lt1)      (ln2,lt1)     (ln3,lt1)     (ln4,lt1)     (ln5,lt1)
#               mean1          mean2         mean3         mean4
#      (ln1,lt2)      (ln2,lt2)     (ln3,lt2)     (ln4,lt2)     (ln5,lt2)
#               mean5          mean6         mean7         mean8
#      (ln1,lt3)      (ln2,lt3)     (ln3,lt3)     (ln4,lt3)     (ln5,lt3)
#               mean9          mean10        mean11        mean12
#      (ln1,lt4)      (ln2,lt4)     (ln3,lt4)     (ln4,lt4)     (ln5,lt4)
#               mean13         mean14        mean15        mean16
#      (ln1,lt5)      (ln2,lt5)     (ln3,lt5)     (ln4,lt5)     (ln5,lt5)
#
# We then take the mean of those means to get the mean of domain/region of each level.
# In doing the mean this way the over-lap in averaging towards the interior values provides
# more weight to those values and hence a more statistically significant mean for the
# the region.
#
TROP = trop[ix_(ri,ci)]
TROP = TROP[1].reshape(5,5) #Hmmm, I FEEL LIKE I'M REALLY NOT UNDERSTANDING PYTHON INDEXING
n,m = TROP.shape
TROP_BAR = average(split(average(split(TROP, m // N, axis=1), axis=-1), n // N, axis=1), axis=-1)
print(TROP_BAR)
OMEGA_BAR = zeros(12)
for i1 in range (0,11):
    oms = om[i1]
    OMS = OMS[ix_(ri,ci)]
    OMS = OMS[1].reshape(5,5)
    OMEGA_BAR[i1] = average(split(average(split(, m // N, axis=1), axis=-1), n // N, axis=1), axis=-1)

Средние значения, которые я получаю, не имеют смысла.Поэтому я хотел бы получить средние значения, которые действительно имеют смысл.Заранее спасибо.

Ответы [ 2 ]

1 голос
/ 11 мая 2019

Хотя я не верю, что это самый эффективный способ. Я нашел решение, которое дает мне правильный ответ, поэтому я опубликую его здесь (ниже). Однако мне было бы любопытно, если бы у кого-нибудь было более эффективное решение, кроме циклического прохождения матрицы. Еще раз спасибо заранее. 

def gridavg(mat,n,k): 
if (k > n): 
    return
# row number of first cell in current 
# sub-square of size k x k
avg = zeros(15)
for i in range(n - k + 1): 
    # column of first cell in current  
    # sub-square of size k x k 
    for j in range(n - k + 1): 
        # Calculate mean of current sub-square 
        for p in range(i, k + i): 
            for q in range(j, k + j):
                avg[i+j+p+q] = mean(mat[p][q])
return avg
0 голосов
/ 14 мая 2019

Можете ли вы проверить этот код: 

import numpy as np
def gridavg(testin):
    testin=np.array(testin)
    test_a=0.5*(testin[:-1,:]+testin[1:,:])
    testout=0.5*(test_a[:,:-1]+test_a[:,1:]);
    return testout

К сожалению, это только для матриц 2x2, но это должно быть намного быстрее, так как для усреднения вы используете массивы NumPy и операции с матрицами.

Для более универсального подхода вы можете попробовать это: 

def gridavg(testin,n,k):
    testin=np.array(testin);
    from_end=1-k;
    if k>0:
        sum_a = None;
        for i0 in range(k):
            if sum_a is None:
                sum_a = np.array(testin[i0:(from_end+n),:]);
            else:
                sum_a = sum_a + testin[i0:(from_end+i0+n),:]
        sum_a = sum_a/float(k);

        sum_b = None;
        for j0 in range(k):
            if sum_b is None:
                sum_b = sum_a[:,j0:(from_end+n)];
            else:
                sum_b = sum_b + sum_a[:,j0:(from_end+j0+n)]
        testout = sum_b/float(k);
    return testout

Я пытался со случайным matrice и для k-значений 2 и 3, казалось, работал.

...