grass.script + scipy theilsen наклон регрессии и перехват между значениями двух растров - PullRequest
Новая
       71

grass.script + scipy theilsen наклон регрессии и перехват между значениями двух растров

0 голосов
/ 18 февраля 2020

Мне нужно вычислить наклон регрессии TheilSen и перехватить значения двух растров в скрипте GRASS GIS python. В этом примере два растра (xtile и ytile) имеют одинаковые размеры 250x250 пикселей и содержат значения nodata (null). До сих пор я использовал один только Grass.script, поэтому я новичок в scipy. Я попытался прочитать некоторые учебные пособия, и на основе этого появился следующий код, который я пробовал в командной строке: 

>>> from scipy import stats
>>> import grass.script as grass
>>> import grass.script.array as garray
>>> x = garray.array(mapname="xtile")
>>> y = garray.array(mapname="ytile")
>>> res_list = stats.theilslopes(y, x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/scipy/stats/_stats_mstats_common.py", line 214, in theilslopes
    deltax = x[:, np.newaxis] - x
MemoryError

Очевидно, это было бы не так просто. Редактировать: Я удалил свои мысли о проблеме размерности массива, я ошибся. Теперь кажется, что размер массива 250x250 - это слишком много. Это так? Есть идеи, как этого избежать?

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

>>> print x
[[  0.   0.   0. ...   0.   0.   0.]
 [  0.   0.   0. ...   0.   0.   0.]
 [  0. 402.   0. ...   0.   0.   0.]
 ...
 [  0.   0.   0. ...   0.   0.   0.]
 [  0.   0.   0. ...   0.   0.   0.]
 [  0.   0.   0. ...   0.   0.   0.]]

, все значения узлов из растра считываются как нули в массиве. В рассматриваемом растре есть большинство пикселей nodata (или null, как указано в GRASS), которые следует игнорировать в регрессии ie. если какое-либо значение в растрах x или y является nodata, соответствующая пара данных x, y не должна использоваться при вычислении регрессии. Можно ли определить значения узлов в массивах так, чтобы они игнорировались напрямую описанным способом, или необходимо сначала отфильтровать пары узлов из пары массивов?

Спасибо.

1 Ответ

0 голосов
/ 20 февраля 2020

Я не уверен, что уместно ответить на мой собственный вопрос, но, возможно, это было бы полезно для других с подобной проблемой.

Мне удалось наконец заставить его работать. Рабочий (т.е. что-то вычисляет) модифицированный пример решает обе проблемы, упомянутые в исходном вопросе: 

>>> from scipy import stats
>>> import grass.script as grass
>>> import grass.script.array as garray
>>> x = garray.array(mapname="xtile").reshape(-1)
>>> y = garray.array(mapname="ytile").reshape(-1)
>>> # Now the null raster values are changed to zeroes. 
>>> # Let us filter them out by pairs of x, y values for any pair which contains a zero value.
>>> xfiltered = np.array([])
>>> yfiltered = np.array([])
>>> i = 0
>>> for xi in np.nditer(x):
...    if x[i] > 0:
...       if y[i] > 0:
...          xfiltered = np.append(xfiltered, [x[i]])
...          yfiltered = np.append(yfiltered, [y[i]])
...    i += 1
... 
>>> # Compute the regression.
>>> res_list = stats.theilslopes(yfiltered, xfiltered)
>>> res_list
(0.8738738738738738, -26.207207207207148, 0.8327338129496403, 0.9155844155844156)

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

Теперь я до сих пор не уверен, зачем нужна была фильтрация для stats.theilslopes для завершения sh без ошибок (хотя со всеми нулями в данных результаты все равно будут неверными). Возможно, что фильтрация нулей просто сделала набор достаточно маленьким, чтобы поместиться в памяти, но я не верю в это. Я думаю, что большинство нулевых значений делает невозможным вычисление медианного наклона пар точек x, y, поскольку, если большинство точек имеют одинаковые значения x, y, то и большинство их пар имеют неопределенный наклон, а медиана равна большинство. Но это может быть что-то совсем другое.

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

Последнее замечание: у меня, вероятно, обратные данные x и y, если y являются зависимыми, а x независимыми переменными. Интуитивно я чувствовал, что они должны быть в порядке y, x, но я вижу, что во всех руководствах и документах это всегда x, y. Я оставил его, как и в первоначальном вопросе, поскольку в некоторых случаях вы можете искать обратную линию формулы x = f (y), и это не является частью проблемы, которую я пытался решить.

...