grass.script+scipy-theilsen回归斜率和两个光栅值之间的截距

Question

我需要在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）像素，在回归中应忽略这些像素，即，如果光栅x或y中的任何值为nodata，则在回归计算中不应使用相应的x、y数据对。是否可以在数组中定义nodata值，以便以所述方式直接忽略这些值，或者是否需要首先从数组对中过滤出nodata对

多谢各位

Answer 1

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

我终于成功了。修改后的示例解决了原问题中提到的两个问题：

>>> 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.slopes不需要使用整形，但它使两个阵列的过滤更加容易

现在，我仍然不知道为什么要完成stats.slopes而不出错就需要进行过滤（尽管数据中有所有的零，结果无论如何都是错误的）。可能是，过滤掉零只会使集合小到足以容纳内存，但我不这么认为。我认为大多数零值使计算x，y点对的中值斜率变得不可能，因为如果大多数点具有相同的x，y值，那么它们的大多数对也具有未定义的斜率，中值约为大多数。但它可能完全是另一回事

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

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