18年NDVI斜率图生成

想自己用python实现多幅栅格图像多年变化率图的制作，但网上始终没有找到合适的代码。

基础太差了55555

贴老师写的多波段图像实现多年降水与气温的月NDVI相关系数图代码：

import numpy as np
from pyrsgis import raster
from scipy.stats.stats import pearsonr

RasterFile1 = 'D:/XXX/NDVIAPR.tif'
RasterFile2 = 'D:/XXX/tempAPR.tif'
ds1, tempArr1 = raster.read(RasterFile1)
ds2, tempArr2 = raster.read(RasterFile2)
nBands1, rows1, cols1 = ds1.RasterCount, ds1.RasterXSize, ds1.RasterYSize
features1 = np.empty((nBands1, cols1, rows1))
features2 = np.empty((nBands1, cols1, rows1))
features1[:, :, :] = tempArr1
features2[:, :, :] = tempArr2
features_result1 = np.empty((cols1, rows1))
features_result2 = np.empty((cols1, rows1))

for i in range(0, cols1):
    for j in range(0, rows1):
        list1 = []
        list2 = []
        for k in range(0, nBands1):
            list1.append(features1[k, i, j])
            list2.append(features2[k, i, j])
        r, p_value = pearsonr(list1, list2)
        features_result1[i, j] = r
        features_result2[i, j] = p_value
outFile1 = 'd:/TEMPcor4.tif'
outFile2 = 'd:/TEMPsig4.tif'
raster.export(features_result1, ds1, filename=outFile1, dtype='float', bands='all')
raster.export(features_result2, ds1, filename=outFile2, dtype='float', bands='all')

其中的np.empty我开始以为是np.zeros的用法，后来发现不是：

np.empty（）：它可以生成一个根据参数规定的数组；若不初始化则会生成随机数。

我觉得要改的话应该就是相同的思路：

第一个是用皮尔逊函数，求十八波段某月NDVI和十八波段某月降水/气温之间的相关性与显著性；

第二个可以看作是把NDVI-->年NDVI，pre\temp-->十八年这个时间序列。

目前为止我尝试过两种思路（对不起TAT是真的太菜了）

1.更改皮尔逊函数，找到可以是一个变量与常数18的线性回归函数（对不起对不起！我也不知道我怎么想的真的太蠢了！可能这就是没学概率论的下场吧。。缺乏常识）

2.更改参与计算的维度，即怎么表达十八年这个时间序列?

我尝试了修改features2[:, :, :] = tempArr1为features2[:,0,0] = tempArr2，list2.append(features2[k, i, j])修改为list2.append(features2[k, 0, 0])。

得到了如下结果（aa是错的，下面的是正确的）

接着去查如何正确用切片切出自己需要的维度：

a[:, :, :]表示，第一个冒号表示选择第几维，也就是第几个数组；第二个表示第几行；第三个表示第几列

详情可见：https://www.cnblogs.com/vincent-sh/p/12730509.html

于是将

features2[:, :, :] = tempArr1 --》 features2[0, :, :] = tempArr1 ，想表达只取第一个维度，也就是十八个波段的那一层；还是错误，得到的图像与上面那张相同

为什么操作不同，最后的结果还是一样呢？？？我太困惑了！！

分分钟破案：没有保存。。。打脸来得太快

#后来发现pycharm原来是不需要保存的哦！写一句存一句的那种

接着就报错了：ValueError: could not broadcast input array from shape (18,184,211) into shape (184,211)

查！“出现这个问题的主要原因是因为list中array的shape不一致造成的，所以发生这个问题的时候，条件允许的话先将list中数组的shape全部打印出来观察一下”（很明显当下不可以）

试了这篇博客https://blog.csdn.net/sinat_29957455/article/details/103487477的第二种用mask的方法，没有研究这个函数到底什么意思，这样

features1[:, :, :] = tempArr1 #三个纬度数组的全部内容

mask = np.zeros((18, 184, 211))
features2_mask = copy.deepcopy(mask)
features2_mask[:features2.shape[0], :features2.shape[1]] = features2
features1_mask = copy.deepcopy(mask)
features1_mask[:features1.shape[0], :features1.shape[1]] = features1

features2[0, :, :] = tempArr2

套了以后，报错：ValueError: could not broadcast input array from shape (18,184,211) into shape (184,211)，跟上一个错误一模一样。

发现本身定义两个features的时候，就把人家维度定的是一样的，干嘛要再弄一次。。然后这样改了一波，

features2 = np.empty((nBands1))
features1[:, :, :] = tempArr1 #三个纬度数组的全部内容

mask = np.zeros((18,184,211))
features2_mask = copy.deepcopy(mask)
features2_mask[:features2.shape[0],:features2.shape[1]] = features2
features1_mask = copy.deepcopy(mask)
features1_mask[:features1.shape[0],:features1.shape[1]] = features1

features2[:] = tempArr2

报了个这个错误：

IndexError: tuple index out of range，超出数组范围。好吧，翻译了也看不懂，或许是mask里的参数设置错了呢？

学习mask：

即掩膜，很熟悉啊！

copy和deepcopy的区别：

看了https://blog.csdn.net/sinat_29957455/article/details/103487477后就觉得arcgis里的copy全是shallowcopy。

另外，docker又是什么鬼？

到这里忽然知道前面的为什么错了，每个像元不只是跟十八年的时间序列有关啊！它们也需要在纵横坐标上表达！那么问题出在哪里呢？

重新审视了这段代码，明白了是如何对两张图象中每个像元具有的值进行相关性分析的。于是有了一个奇奇怪怪的想法：能不能构建一个全部像元值相同，且随波段递增的十八波段图像呢？拿出了终极问题：如何arcgis中批量修改像元值？这什么鬼问题啊喂，怎么觉得似曾相识？？

成功改变了一个shape file的属性值，那又如何？转成栅格了之后还是啥也没有了。。先睡一觉，明天一定可以的！

-----------------------------------------------------------------------------分割线--------------------------------------------------------------------------

看了一篇博客https://blog.csdn.net/weixin_39559804/article/details/110078530?ops_request_misc=%25257B%252522request%25255Fid%252522%25253A%252522160842401516780274069014%252522%25252C%252522scm%252522%25253A%25252220140713.130102334.pc%25255Fall.%252522%25257D&request_id=160842401516780274069014&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~first_rank_v2~rank_v29-3-110078530.nonecase&utm_term=arcgis%E4%BF%AE%E6%94%B9%E5%83%8F%E5%85%83%E5%80%BCpython忽然想起栅格并没有属性表呀。。

然后明白了：

“1.多波段栅格是不会有属性表的。
2.单波段栅格，如果是浮点型，也不会有属性表。
3.单波段栅格，如果是整型栅格且没有属性表，可以通过ArcToolbox---数据管理工具---栅格---栅格属性---构建栅格属性表创建。”

好好地查着怎么批量修改属性表，怎么又走到clone这条山路TAT

查查查，栅格计算器好像还不错，http://blog.sina.com.cn/s/blog_627956c90101bfvl.html，过程中发现帮助文档很清楚地说有一个功能可以创建常量栅格！

https://blog.csdn.net/dsac1/article/details/37876381比较有逼格的博主，应该不会再错了吧！果然是见识短浅，很多个擦边球关键词都没有查到常量栅格这个东西，导致花费了更多时间去学了其他东西；也算不后悔吧！

想批量创建常量栅格、批量按掩膜提取的，奈何太菜，都是觉得自己写的是对的，结果又都运行不了。。。还试了C++的代码，运行不了是因为没有c++代码需要的库，于是乎喜提VScode赠送的vcpkg...当然这次都是没有用上的

老老实实创建了十八年的常量栅格，成功合成多波段图了以后，分析出来的值和老师的不一样！我真的不知道该怎么办了！55555TAT

-----------------------------------------------------分割线----------------------------------------------------------------

找到了一个提取NDVI的视频！完整回顾一下

处理过程：

辐射定标

掩膜

提取

计算植被覆盖度

----------------------------------------------分割线--------------------------------------------------------------------------

圣诞节去找了老师，值不同但长相相似的真相大白：

我的用了皮尔逊函数，所以是两个变量之间的相关关系；而变化趋势，指的是一个变量的变化情况，用的是它的斜率！所以这俩本质上是不同的。

相关性：其实可以从值看出来的

变化率：

只说图的话，看起来相关性对比比较强烈(￣(工)￣)

------------------------------------------关于anaconda的小（大）插曲----------------------------------------------------------

平安夜那天去找老师的时候，老师说anaconda3.8太新了，不够稳定，换成3.7的，结果这一换就是俩小时；我的电脑莫名其妙就不支持其中的一系列库，GDAL、numpy等等报错，遇到了老师也解决不了的问题。。回到寝室后自己慢慢琢磨，原来有安的库和anaconda版本始终不匹配的问题在里面，

ERROR: GDAL-3.1.4-cp39-cp39-win_amd64.whl is not a supported wheel on this platform.

ERROR: GDAL-3.1.4-cp37-cp37m-win_amd64.whl is not a supported wheel on this platform.
#总之就是因为不会看名字中的玄机，也不知道有命令可以读出电脑中的python对应支持的库的版本，就镇傻不拉几地一个一个试。。

cpXX里的XX表示的就是对应的python版本；而看自己电脑上的python支持的库的版本，（之前）可以使用

python -m pip debug --verbose

来查看自己电脑上的库的，但是之后突然就不行了TAT

然后又去查其他可以看支持包版本的命令，找到了一个用来更新pip的命令：

python -m pip install --upgrade pip

把电脑上的19.1.1---->20.3.3，再执行上一条就又可以看了！出来的结果长这样

回到圣诞节，老师给了生成斜率的函数，原来是他自己写的哦

import numpy as np

def trendline(data):
    order=1
    index=[i for i in range(1,len(data)+1)]
    coeffs = np.polyfit(index, list(data), order)
    slope = coeffs[-2]
    return float(slope)

调用就van事。不够之前有查到一个回归序列函数，叫lingress，改天试试吧；看来还是编程素养不行，自己编写函数的意识不够

结果出来和一开始那张图不一样？？？！！（得知真相后笑死）

怎么回事呢？哪个是正确的呢？前面那个错哪了呢？

老师坚持要查个水落石出！（步骤有点忘了。。或许每次请教老师的第一件事就是打开录屏）

先把出来的斜率栅格图转成点图层，这样每个栅格对应了一个点，图层就有了属性表；

然后提取一个点要素到一个新图层，把斜率图的十八波段数据多值提取到这个点里

把之前的多波段放进去，再把这个点的属性表导出到excel，转置后就有了十八年这个点的NDVI序列。（应该就是这样的吧？X）

但是我的arcgis的多值提取到点功能出毛病了。。。用地理处理里的提取出来出了第一年，十七年的栅格值全是一样的，试了几次才发现源数据就TM是错的，用了工具箱里的才成功了。

........记不得了，步骤过于混乱，到最后甚至不知道老师在点什么QAQ还手输每个像元值算斜率，结果出来比较大小过后最后才反应过来开始那张是2001-2019的，现在这张是2001-2018的图，当然不一样了。。。

不断地遇到问题，不断地解决，还是这样学得更快！虽然过程中挫败感无比强烈。

这是一篇写了一周的博客，甚至逻辑还很混乱，因为水平垃圾，犯错太多，记不得自己犯了哪些，也不是每一个都能找到或者学会并讲出解决方法。daydayup,加油哇