Как применить анализ корреляции Пирсона ко всем парам пикселей DataArray в качестве матрицы корреляции?

Question

Я сталкиваюсь с серьезными трудностями в создании корреляционной матрицы (пиксель за пикселем) одного Netcdf с размерами ('lon', 'lat', 'time'). Мое окончательное намерение состоит в том, чтобы сгенерировать то, что называют Картой электросвязи.

Эта Карта составлена ​​из коэффициентов корреляции. Каждый пиксель имеет значение, которое представляет наибольшее значение корреляции (в модуле), найденное в матрице корреляции, среди всех пар пикселей в массиве данных.

Поэтому для создания моей карты электросвязи вместо циклического повторения каждогодолгота ('долгота') и каждая широта ('широта'), а затем, проверяя все возможные комбинации корреляции, для которых была более высокая величина, я думал о применении функции xr.apply_ufunction с обернутой корреляционной функцией внутри.

Несмотря на все мои усилия, я все еще не понимаю, что действительно происходит за кулисами в xr.apply_ufunc. Все, что мне удалось сделать, это как одна Результирующая матрица со всеми пикселями, равными 1 (идеальная корреляция).

См. Код ниже:

---

import numpy as np

import xarray as xr

def correlation(x, y):

    return np.corrcoef(x, y)[0,0] # to return a single correlation index, instead of a matriz


def wrapped_correlation(da, x, coord='time'):
    """Finds the correlation along a given dimension of a dataarray."""


    from functools import partial

    fpartial = partial(correlation, x.values)

    return xr.apply_ufunc(fpartial, 
                          da, 
                          input_core_dims=[[coord]] , 
                          output_core_dims=[[]],
                          vectorize=True,
                          output_dtypes=[float]
                          )



# testing the wrapped correlation for a sample data:

ds = xr.tutorial.open_dataset('air_temperature').load()

# testing for a single point in space.

x = ds['air'].sel(dict(lon=1, lat=92), method='nearest')

# over all points in the DataArray

Corr_over_x = wrapped_correlation(ds['air'], x)

Corr_over_x# notice that the resultant DataArray is composed solely of ones (perfect correlation match). This is impossible. I would expect to have different values of correlation for each pixel in here

# if one would plot the data, I would be composed of a variety of correlation values (see example below):

Corr_over_x.plot()

Это важный актив дляметеорологи и исследования по дистанционному зондированию. Это позволяет оценить потенциальные геофизические закономерности в данной области изучения.

Я благодарю вас за ваше время и надеюсь, что вы скоро услышите.

С уважением,

Answer 1

Во-первых, вам нужно использовать np.corrcoef(x, y)[0,1]. В конце концов, вам вообще не нужно использовать partial, см. Ниже:

def correlation(x1, x2):

    return np.corrcoef(x1, x2)[0,1] # to return a single correlation index, instead of a matriz


def wrapped_correlation(da, x, coord='time'):
    """Finds the correlation along a given dimension of a dataarray."""

    return xr.apply_ufunc(correlation, 
                          da, 
                          x,
                          input_core_dims=[[coord],[coord]] , 
                          output_core_dims=[[]],
                          vectorize=True,
                          output_dtypes=[float]
                          )

Answer 2

Мне удалось решить мой вопрос. Сценарий стал немного длинным. Тем не менее, он делает то, что предполагалось ранее.

Код адаптирован из этой ссылки

Так как здесь слишком долго показывать фрагмент, я публикуюссылка на мою учетную запись Github, в которой алгоритм (организованный в пакете с именем Teleconnection_using_xarray_data) можно проверить здесь .

В пакете есть два модуля с похожими результатами.

Первый модуль (teleconnection_with_connecting_pathways) работает медленнее, чем второй (teleconnection_via_numpy), но он позволяет оценивать соединительные пути между частичными картами телесвязи.

Второй, возвращает только результирующую карту телеконнекциибез соединительных линий (геопандас-линии строк), хотя это намного быстрее.

Не стесняйтесь сотрудничать. Если возможно, я бы хотел объединить оба модуля, обеспечивающие анализ скорости и пути прохождения, в алгоритме Teleconnection.

С уважением,

Филипп Лил