讲完了geotiff格式数据的读取和保存，本文讲下怎么用python处理一系列的栅格数据（本文以时间序列为例）。
假设我们有某地区每一年的降水序列，一共几十年，现在想要得到每个像元上年降水的变化趋势以及趋势的显著性检验（得到P值），怎么做呢？
思路
对于一个栅格数据，其包括元信息 数据。我们求每个像元上年降水的变化趋势以及对应的P值，实际上只是对数据进行处理，元信息基本上是不变的。
在处理的过程中，我们是求每个像元在时间维度上的变化趋势，类似下图：
引用自arcgis网站
也就是说我们对上图中的每一个条柱时间序列求趋势即可。有了思路，就非常简单了，我们直接上代码。
数据创建
这里说的数据创建就是把我们的多个栅格序列组成上面类似的时空立方体（这里仅指时空栅格）。
之前我们讲了怎么读取单张栅格，读取完之后是一个numpy的ndarray，那么只要进行相应的矩阵拼接即可：
- 导入包
import rasterio
import scipy.stats as ss
import numpy as np
import glob
from rasterio.plot import show
- 读取第一张图的元数据，方便最后写出结果
with rasterio.open(rs_files[0]) as src:
show(src)
profile = src.profile
data = src.read()
- 组成时空立方体
这里组成栅格立方体也非常简单，用numpy的concatenate函数拼接一下即可：
代码语言：javascript复制ds = [rasterio.open(f).read() for f in rs_files]
ds = np.concatenate(ds,axis=0)
ds.shape
>> (36, 133, 110)
show(ds[1,:,:])
趋势和p值计算
前面说过只要对每个条柱时间序列进行趋势计算即可，那么如何同时对所有的条柱时间序列进行计算呢？大神们已经为我们准备好了相应的工具，这就是numpy的apply_along_axis函数，具体的见参考链接【2】。
简单说就是这个函数可以沿着某一个维度应用我们定义的函数。所以我们只需要定义好相应的函数，沿着时间维度应用即可：
- 定义函数
def rs_slope(y,method='linear'):
"""
Compute the slope of a series of raster data along a given axis(always be time).
"""
l = len(y)
ys = range(l)
try:
if method == 'linear':
slope, intercept, r_value, p_value, std_err = ss.linregress(ys,y)
except:
slope = np.nan
return slope, r_value**2, p_value
- 应用函数
slope_rs,r2,pv = np.apply_along_axis(rs_slope, 0,arr=ds)
- 保存结果
profile['nodata'] = 0
with rasterio.open('slope.tif','w',**profile) as dest:
dest.write(slope_rs,1)
- 查看结果
with rasterio.open('./slope.tif') as src:
show(src)
到这里就完成了每个像元的线性趋势计算，不过上面的代码只保存了趋势值，并没有保存R方和p值，读者根据代码改一下即可。
小补充
为什么上面没有写年份？ 因为在计算趋势的时候，如果你不关心截距，那么年份是从0-35还是1980-2015，你算出来的趋势值（也就是下面公式中的a，x是年份）都是一样的，那么就不必要多浪费那点算力了：
总结
- 处理栅格序列的时候，元信息一般不变，所以可以利用某一个原始数据的元信息作为模版，方便保存处理后的结果；
- 对于栅格数据的值，就是一个数组而已，巧用numpy的函数可以实现很多我们需要的功能；
- 在能简化算法的时候，尽量简化。
参考
【1】
https://pro.arcgis.com/zh-cn/pro-app/latest/tool-reference/big-data-analytics/create-space-time-cube.htm
【2】https://numpy.org/doc/stable/reference/generated/numpy.apply_along_axis.html
