基于python的Hurst计算预测未来发展趋势（长时序栅格影像）-编程知识

基于python的Hurst计算预测未来发展趋势（长时序栅格影像）

1.Hurst指数反映了时间序列长期记忆性的程度，即过去的信息对未来的影响程度。Hurst指数的取值范围为0到1之间，当Hurst指数等于0.5时，时间序列被认为是一种随机漫步，即具有随机性；当Hurst指数大于0.5时，时间序列显示出长期正相关性，表明趋势在未来可能持续；当Hurst指数小于0.5时，时间序列显示出长期负相关性，表明趋势在未来可能反转。

2.下面是一个例子，计算下面22幅影像的Hurst。
在这里插入图片描述
计算结果如下：

3.理论不再进一步讲解，直接上代码，有基础的可以自己调试，小白可以使用本人编好的exe程序，链接在下面。

import numpy as np
from osgeo import gdal
import osdef write(file_name, image, projection,geotransform,x_size,y_size):dtype = gdal.GDT_Float32# 数据格式driver = gdal.GetDriverByName('GTIFF')# 创建数据,设置文件路径及名称new_ds = driver.Create(file_name, x_size, y_size, 1, dtype)# 设置投影信息及6参数new_ds.SetGeoTransform(geotransform)new_ds.SetProjection(projection)# 将值写入new_ds中new_ds.GetRasterBand(1).WriteArray(image)# 把缓存数据写入磁盘new_ds.FlushCache()del new_dsdef Hurst(x):# x为numpy数组n = x.shape[0]t = np.zeros(n - 1)  # t为时间序列的差分for i in range(n - 1):t[i] = x[i + 1] - x[i]mt = np.zeros(n - 1)  # mt为均值序列,i为索引,i+1表示序列从1开始for i in range(n - 1):mt[i] = np.sum(t[0:i + 1]) / (i + 1)# Step3累积离差和极差,r为极差r = []for i in np.arange(1, n):  # i为taocha = []for j in np.arange(1, i + 1):if i == 1:cha.append(t[j - 1] - mt[i - 1])if i > 1:if j == 1:cha.append(t[j - 1] - mt[i - 1])if j > 1:cha.append(cha[j - 2] + t[j - 1] - mt[i - 1])r.append(np.max(cha) - np.min(cha))s = []for i in np.arange(1, n):ss = []for j in np.arange(1, i + 1):ss.append((t[j - 1] - mt[i - 1]) ** 2)s.append(np.sqrt(np.sum(ss) / i))r = np.array(r)s = np.array(s)xdata = np.log(np.arange(2, n))ydata = np.log(r[1:] / s[1:])h, b = np.polyfit(xdata, ydata, 1)return hdef main(path1,result_path):filepaths = []for file in os.listdir(path1):filepath1 = os.path.join(path1, file)filepaths.append(filepath1)# 获取影像数量num_images = len(filepaths)# 读取影像数据img1 = gdal.Open(filepaths[0])# 获取影像的投影，高度和宽度transform1 = img1.GetGeoTransform()proj = img1.GetProjection()height1 = img1.RasterYSizewidth1 = img1.RasterXSizearray1 = img1.ReadAsArray(0, 0, width1, height1)del img1# 读取所有影像for path1 in filepaths[1:]:if path1[-3:] == 'tif':img2 = gdal.Open(path1)array2 = img2.ReadAsArray(0, 0, width1, height1)array1 = np.vstack((array1, array2))del img2array1 = array1.reshape((num_images, height1, width1))# 输出矩阵，无值区用nan填充h_array = np.full([height1, width1], np.nan)# 只有有值的区域才进行计算c1 = np.isnan(array1)sum_array1 = np.sum(c1, axis=0)nan_positions = np.where(sum_array1 == num_images)positions = np.where(sum_array1<10)for i in range(len(positions[0])):x = positions[0][i]y = positions[1][i]hurst_list1 = array1[:, x, y]hurst_list1=hurst_list1[~np.isnan(hurst_list1)]h = Hurst(hurst_list1)h_array[x, y] = hwrite(result_path, h_array, proj, transform1, width1, height1)
if __name__ == "__main__":path1 = r"F:\rsei-2"result_path = r"F:\jieguo\h.tif"main(path1,result_path)