0

我需要计算 GRASS GIS python 脚本中两个栅格值之间的 TheilSen 回归斜率和截距。此示例中的两个栅格（xtile 和 ytile）都具有相同的尺寸 250x250 像素并包含 nodata（空）值。到目前为止，我只使用了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 值都被读取为零到数组中。在所讨论的栅格中，大多数 nodata （或在 GRASS 中命名的 null ）像素，在回归中应该被忽略，即。如果栅格 x 或 y 中的任何值是 nodata，则不应在回归计算中使用相应的 x,y 数据对。是否可以在数组中定义 nodata 值，以便以所描述的方式直接忽略这些值，或者是否需要首先从数组对中过滤掉 nodata 对？

谢谢你。

4

1 回答 1

0

我不确定回答我自己的问题是否合适，但也许对其他有类似问题的人有用。

我终于让它工作了。工作（即它计算某些东西）修改示例解决了原始问题中提到的两个问题：

>>> 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)

说明：在回归之前，我过滤掉了原始栅格中的所有零值（也可能是负值，不应该有，如果它们是，这将意味着有缺陷的数据 - 数据的物理意义是量化反射率）计算。stats.theilslopes 不需要使用 reshape，因为我用测试数据进行了实验验证，但它使两个数组的过滤变得更加容易。

现在，我仍然不确定为什么 stats.theilslopes 需要过滤才能完成而没有错误（尽管数据中的所有零结果都是错误的）。可能是，过滤掉零只是使集合足够小以适合内存，但我不这么认为。我认为大多数零值使得无法计算 x,y 点对的中值斜率，因为如果大多数点具有相同的 x,y 值，那么它们的大多数对也有未定义的斜率，中值约为多数。但它可能完全是另一回事。

另外，由于我不是非常熟练的 Python 程序员，也许我这样做的方式不是最有效的方式。这可以由其他人纠正。

最后一点，如果 y 是因变量，x 是自变量，我可能将 x 和 y 数据颠倒过来。直觉上我觉得它们应该是 y,x 的顺序，但是我在所有的教程和文档中看到它总是 x,y。我将其保留在原始问题中，因为在某些情况下您可以搜索逆公式 x=f(y) 回归线，这不是我试图解决的问题的一部分。

于 2020-02-20T11:31:46.983 回答