实现自相关分析的总体步骤

1、在mtalab中将NDVI数据转化为mat数据

folder_path = 'D:\1KM NDVI数据\裁剪数据';
 
% 获取文件夹下所有的 .tif 文件
tif_files = dir(fullfile(folder_path, '*.tif'));
 
% 读取第一个文件的经纬度信息
info = geotiffinfo(fullfile(folder_path, tif_files(1).name));
pop1990 = imread(fullfile(folder_path, tif_files(1).name));
 
% 获取各格点经纬度
lat = zeros(size(pop1990, 1), size(pop1990, 2));
lon = zeros(size(pop1990, 1), size(pop1990, 2));
x = lat;
y = lon;
 
for i = 1 : size(pop1990, 1)
    for j = 1 : size(pop1990, 2)
        [x(i, j), y(i, j)] = pix2map(info.RefMatrix, i, j);
    end
end
 
[lat, lon] = projinv(info, x, y);
 
% 将经纬度转换为一列
lat_col = lat(:);
lon_col = lon(:);
 
% 初始化结果数组
result = [lat_col, lon_col];
 
% 遍历所有 .tif 文件并读取数据
for k = 1:length(tif_files)
    % 读取当前文件的图像数据
    tif_path = fullfile(folder_path, tif_files(k).name);
    tif_data = imread(tif_path);
    
    % 将不在 -1 到 1 范围内的数据设置为 0
    tif_data(tif_data < -1 | tif_data > 1) = 0;
    
    % 将图像数据转换为一列
    tif_data_col = tif_data(:);
    
    % 将数据添加到结果数组中
    result = [result, tif_data_col];
end
 
% 保存结果到文件
save('C:\Users\y8143\Desktop\rs投稿\result.mat', 'result');

2、NDVI数据转换为R语言中的data数据框格式

# 安装 R.matlab 包（如果尚未安装）
# install.packages("R.matlab")
 
# 加载 R.matlab 包
library(R.matlab)
 
# 读取 MAT 文件
mat_data <- readMat("C:/Users/y8143/Desktop/rs投稿/result.mat")
 
# 查看读取的数据结构
str(mat_data)
 
# 提取数据
result <- mat_data$result
 
# 确保 result 是一个矩阵或数据框，并且有 22 列
if (ncol(result) != 22) {
  stop("The result matrix does not have 22 columns as expected.")
}
 
# 创建列名向量
col_names <- c("lat", "lon", as.character(2000:2019))
 
# 将 result 转换为数据框并设置列名
df <- as.data.frame(result)
names(df) <- col_names
 
# 查看数据框的前几行
head(df)
 
# 保存数据框到文件（可选）
write.csv(df, file = "C:/Users/y8143/Desktop/rs投稿/result.csv", row.names = FALSE)

3、在r语言中生成结果，并将结果输出为CSV格式

library(fitAR)  # 确保你有 fitAR 包
 
# 假设 data 是你的数据框，前两列是经纬度，后面所有列是 NDVI 数据
# 例如：
# data <- read.csv("path_to_your_data.csv")
 
# 时间点
years <- 2000:2019
time <- 1:length(years)
 
# 创建一个新的数据框来存储结果
results <- data.frame(
  lon = data$lon,
  lat = data$lat,
  AR_coefficient = NA,
  AR_pvalue = NA,
  Time_trend_coefficient = NA,
  Time_trend_pvalue = NA
)
 
# 对每一行进行自回归分析
for (i in 1:nrow(data)) {
  # 提取当前行的 NDVI 数据
  ndvi <- as.numeric(data[i, 3:ncol(data)])  # 提取 NDVI 列数据
  
  # 创建数据框用于拟合
  df <- data.frame(ndvi = ndvi, time = time)
  
  # 拟合 AR 模型，使用时间作为协变量
  AR_model <- fitAR(ndvi ~ time, data = df)
  
  # 提取结果
  results$AR_coefficient[i] <- AR_model$coefficients["time"]
  results$AR_pvalue[i] <- AR_model$pval["time"]
  results$Time_trend_coefficient[i] <- AR_model$coefficients["(Intercept)"]
  results$Time_trend_pvalue[i] <- AR_model$pval["(Intercept)"]
}
 
# 查看结果
print(results)

4、在matlab中将结果再转化为TIF

>> % 2. 提取经纬度和图像数据列
lat_col = result(:, 1);
lon_col = result(:, 2);
data_col = result(:, 3);

% 3. 确定经纬度范围和间隔（基于先前信息）
lat_min = min(lat_col);
lat_max = max(lat_col);
lon_min = min(lon_col);
lon_max = max(lon_col);

% 假设纬度和经度是等间距的
num_rows = 701;  % 从原始矩阵的纬度维度
num_cols = 928;  % 从原始矩阵的经度维度

% 4. 将列向量重新构建为原始的二维矩阵
lat_grid = reshape(lat_col, num_rows, num_cols);
lon_grid = reshape(lon_col, num_rows, num_cols);
data_grid = reshape(data_col, num_rows, num_cols);

% 5. 定义输出文件名
filename = 'C:\Users\y8143\Desktop\rs投稿\HUAIHE.tif';
imwrite(data_grid,filename)

上述结果还是没有经纬度坐标，在matlab补充经纬度坐标

% 1. 读取 2001.TIF 文件的地理信息
filename_source = 'D:\关于毕业论文\1KM NDVI数据\裁剪数据\2001.TIF';
[data_source, R] = readgeoraster(filename_source);
info=geotiffinfo('D:\关于毕业论文\1KM NDVI数据\裁剪数据\2001.TIF');

% 2. 读取要保存的数据矩阵
filename_source2='C:\Users\y8143\Desktop\rs投稿\HUAIHE.tif'
data_grid = readgeoraster(filename_source2)

% 3. 定义 HUAIHE.tif 的文件路径
filename_dest = 'C:\Users\y8143\Desktop\rs投稿\HUAIHE3.tif';

% 4. 使用 geotiffwrite 保存包含从 2001.TIF 复制的地理信息的TIFF文件
 geotiffwrite(filename_dest,data_grid,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

其实总结还原为tif的过程，必须有原始参考文件。