多年遥感TIF图像计算HURST指数和变异系数CV（matlab和python代码）

Hurst指数python代码：

import numpy as np
import rasterio
import os
from multiprocessing import Pool

eps = 1e-10

def read_geotiff(file_path):
    with rasterio.open(file_path) as src:
        data = src.read(1)
        meta = src.meta
    return data, meta

def calculate_hurst(ndvi):
    if np.min(ndvi) > 0:
        ndvi_diff = np.diff(ndvi)
        mean_diff = np.array([np.mean(ndvi_diff[:i + 1]) for i in range(len(ndvi_diff))])
        std_diff = np.array([np.std(ndvi_diff[:i + 1]) * np.sqrt(i / (i + 1)) for i in range(len(ndvi_diff))])
        rr = np.array([np.max(np.cumsum(ndvi_diff[:i + 1] - mean_diff[i])) - np.min(np.cumsum(ndvi_diff[:i + 1] - mean_diff[i])) for i in range(len(ndvi_diff))])
        rs = std_diff[1:] / (rr[1:] + eps)

        if np.any(np.isnan(rs)) or np.any(np.isinf(rs)):
            return np.nan

        lag = np.arange(2, len(ndvi_diff) + 1)
        valid_idx = (rs > 0) & (lag > 0)
        if np.sum(valid_idx) < 2:
            return np.nan

        H, _ = np.polyfit(np.log(lag[valid_idx]), np.log(rs[valid_idx]), 1)
        return H
    return np.nan

def main():
    base_path = r'G:\MOD16A3(2014-2024)\xn-4' #需要修改为你图像的目录
    initial_file = os.path.join(base_path, '2015.tif') #需要加载其中一张图像作为底图
    aa, meta = read_geotiff(initial_file)

    ndvisum = np.zeros((aa.size, 10))

    for year in range(2014, 2024):
        file_path = os.path.join(base_path, f'{year}.tif')
        with rasterio.open(file_path) as src:
            ndvi = src.read(1)
            ndvi = ndvi.flatten()
        ndvisum[:, year - 2014] = ndvi

    # 使用多进程
    with Pool() as pool:
        hsum = pool.map(calculate_hurst, ndvisum)

    hsum = np.array(hsum).reshape(aa.shape)

    output_file = os.path.join(base_path, 'Hurst.tif') #输出文件，可以添加需要的目录
    with rasterio.open(output_file, 'w', **meta) as dst:
        dst.write(hsum, 1)

if __name__ == '__main__':
    main()

变异系数python代码：

from osgeo import gdal
import numpy as np
import copy
import os


# 计算变异系数的函数
def coefficient_of_variation(data):
    mean = np.mean(data)
    std = np.std(data, ddof=0)
    cv = std / mean
    return cv


# 栅格图像组计算变异系数
def CV(images, outpath):
    images_pixels = []
    for image in images:
        open_tif = gdal.Open(image)
        band = open_tif.ReadAsArray()
        images_pixels.append(band)
    CV = copy.deepcopy(images_pixels[0])

    for i in range(len(CV)):
        for j in range(len(CV[i])):
            CV_data = [images_pixels[px][i][j] for px in range(len(images))]
            CV_value = coefficient_of_variation(CV_data)
            CV[i][j] = CV_value

    gtiff_driver = gdal.GetDriverByName('GTiff')
    out_tif = gtiff_driver.Create(outpath, CV.shape[1], CV.shape[0], 1, gdal.GDT_Float32)
    out_tif.SetProjection(open_tif.GetProjection())
    out_tif.SetGeoTransform(open_tif.GetGeoTransform())
    out_band = out_tif.GetRasterBand(1)
    out_band.WriteArray(CV)
    out_band.FlushCache()
    print('Raster CV computed')


# 自动读取文件夹中的TIF文件
def get_tif_files(folder):
    years = range(2014, 2024)  # 2014到2023，修改为你所对应的年份序列
    tif_files = [os.path.join(folder, f"{year}.tif") for year in years]
    return tif_files


# 主程序入口
folder_path = r'G:\MOD16A3(2014-2024)\xn-4' #修改为你图像的目录
tif_files = get_tif_files(folder_path)
output_path = os.path.join(folder_path, 'MOD16-cv.tif')
CV(tif_files, output_path)

注意：需要添加相关的运行库才能运行。

HURST指数MATLAB代码：

clear
[aa,R]=geotiffread('G:/2014.tif');%先投影信息
info=geotiffinfo('G:/2014.tif');
ndvisum=zeros(size(aa,1)*size(aa,2),10);%n期数据
for year=2014:2023
    filename=strcat('G:/',int2str(year),'.tif');
    ndvi=importdata(filename);
    ndvi=reshape(ndvi,size(ndvi,1)*size(ndvi,2),1);
    ndvisum(:,year-2013)=ndvi;%year=2014时应该减2013
end
hsum=zeros(size(aa,1),size(aa,2))+NaN;
for kk=1:size(ndvisum,1);
    ndvi=ndvisum(kk,:);
    if min(ndvi)>0
        ndvi_cf=[];
        for i=1:length(ndvi)-1
            ndvi_cf1=ndvi(i+1)-ndvi(i);
            ndvi_cf=[ndvi_cf,ndvi_cf1];
        end
        M=[];
        for i=1:size(ndvi_cf,2)
            M1=mean(ndvi_cf(1:i));
            M=[M,M1];
        end
        S=[];

        for i=1:size(ndvi_cf,2)
            S1=std(ndvi_cf(1:i))*sqrt((i-1)/i);
            S=[S,S1];
        end

        for i=1:size(ndvi_cf,2)
            for j=1:i
                der(j)=ndvi_cf(1,j)-M(1,i);
                cum=cumsum(der);
                RR(i)=max(cum)-min(cum);
            end
        end

        RS=S(2:size(ndvi_cf,2)).\RR(2:size(ndvi_cf,2));
        T=[];
        for i=1:size(ndvi_cf,2)
            T1=i;
            T=[T,T1];
        end
        lag=T(2:size(ndvi_cf,2));                  
        g=polyfit(log(lag/2),log(RS),1);                
        H=g(1);
        hsum(kk)=H;
        clear der
    end
end
geotiffwrite('G:/MOD16A3(2014-2024)/xn-4/WGS84/Hurst.tif',hsum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

变异系数MATLAB代码：

% 导入第一年的图像以获取空间参考和信息
[a, R] = geotiffread('G:\MOD16A3(2014-2024)\xn-4\2014.tif');
info = geotiffinfo('G:\MOD16A3(2014-2024)\xn-4\2014.tif');

% 初始化参数
[m, n] = size(a);
years = 10; % 从2014到2023年，如果你是2000-2020，则这个值为21，根据你的年份长度而定
data = zeros(m*n, years);
k = 1;

% 读取每年的数据并转化为合适的形式
for year = 2014:2023
    file = fullfile('G:\MOD16A3(2014-2024)\xn-4', sprintf('%d.tif', year));
    if exist(file, 'file')
        bz = geotiffread(file);
        bz = reshape(bz, m*n, 1);
        data(:, k) = bz;
    else
        disp(['File does not exist: ', file]);
    end
    k = k + 1;
end

% 计算变异系数
bya = zeros(m, n);

for i = 1:m*n
    bz = data(i, :);
    if min(bz) > 0 % 根据数据实际最小值范围更改
        A = mean(bz); % 求平均值
        S = std(bz, 0, 2); % 求标准差，第二个参数0表示样本标准差
        V = S / A; % 变异系数
        bya(i) = V;
    end
end

bya = reshape(bya, m, n); % 重塑bya为原始图像尺寸

% 保存结果为GeoTIFF文件
name1 = 'G:\MOD16A3(2014-2024)\xn-4\变异系数.tif';
geotiffwrite(name1, bya, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);