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

发布时间:2024年01月12日

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 os

def 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_ds

def 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为tao
        cha = []
        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 h


def 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.RasterYSize
    width1 = img1.RasterXSize
    array1 = 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 img2

    array1 = 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] = h

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

代码中导入了GDAL库,使用时需要安装该库。有基础的可以使用上述代码,自己调试,exe程序的链接如下,可以下载使用。

我用夸克网盘分享了「Hurst持续性检测.exe」,点击链接即可保存。打开「夸克APP」,无需下载在线播放视频,畅享原画5倍速,支持电视投屏。
链接:https://pan.quark.cn/s/83f34c91bc2e
提取码:MJ1L

希望大家可以关注下方的微信公众号:RS GIS遥感 地信学习,里面有更多的技术分享及软件使用教程,本人在这里谢谢大家,软件有什么问题可以微信公众号留言。感谢!!!!
在这里插入图片描述

文章来源:https://blog.csdn.net/weixin_43887139/article/details/135550594
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。