基于 matlab 计算 TPI（地形位置指数）

1. TPI 简介

        地形位置指数算法由 Weiss 提出，主要是根据局部地形高程对各类地貌单元提取。

        其基本原理为：在邻域分析方法的基础上，计算每个栅格的高程值和该栅格领域内所有栅格的平均高程之间的差值，正值表示该栅格点高于领域内栅格，而负值意味着它低于领域内栅格。当TPI值接近于零值时，该栅格点可能位于平面或斜面，如下图所示：

TPI值与地形的关系

2. TPI 计算

        TPI 计算公式如下：

        其中：为点处的高程值，代表分析窗口半径，代表分析窗口半径内所有栅格点数。

3. Arcgis 计算 TPI

        参考链接：利用ArcGIS计算TPI(地形位置指数)_哔哩哔哩_bilibili

4. Matlab 计算 TPI

%% 读取dem数据
% [dem_mat,dem_count]=geotiffread('cut_dem.tif');
[dem_mat,dem_count]=readgeoraster('cut_dem.tif');
info=geotiffinfo('cut_dem.tif');
dem_mat=double(dem_mat);
dem_mat(abs(dem_mat)>10000)=0;
% imagesc(dem_mat)

%% 计算TPI
% 初始化TPI矩阵
tpi_mat = zeros(size(dem_mat));
% 设置邻域窗口大小
win_size = 3;
% 计算窗口的半径
radius = floor(win_size / 2);
% 获取DEM数据的行数和列数
[rows, cols] = size(dem_mat);
% 遍历DEM数据中的每个点
for i = radius+1:rows-radius
    for j = radius+1:cols-radius
        % 提取3x3邻域
        window = dem_mat(i-radius:i+radius, j-radius:j+radius);
        % 计算邻域内高程的平均值
        mean_dem = mean(window(:));
        % 计算TPI
        tpi_mat(i, j) = dem_mat(i, j) - mean_dem;
    end
end
% 显示TPI矩阵
imagesc(tpi_mat);
colorbar;

% 保存tif影像
tpi_mat=int16(tpi_mat);
geotiffwrite('tpi.tif',tpi_mat,dem_count,'GeoKeyDirectoryTag',...
    info.GeoTIFFTags.GeoKeyDirectoryTag);