MATLAB 确定到目标网格的路径上经过的所有像元点

任务描述：确定站点i到目标网格的路径上经过的所有像元点

参考链接如下：
A*寻路 – 更加真实 的路径（一）
algorithm（一）：直线经过的网格区域计算

clc
clear
close all

[DEM,Georef] = geotiffread('F:\Data of Flash Flood\PRE_interpolation_test\Catc_DEM_resa.tif');  %读取目标区域DEM数据

%% 确定区域内像元点边界的经纬度和中心点的经纬度
% 像元点的边界的经纬度
[x,y] = size(DEM);
lon_boundary = Georef.LongitudeLimits(1):(Georef.LongitudeLimits(2)-Georef.LongitudeLimits(1))/y:Georef.LongitudeLimits(2);
lat_boundary = Georef.LatitudeLimits(1):(Georef.LatitudeLimits(2)-Georef.LatitudeLimits(1))/x:Georef.LatitudeLimits(2);
% 像元点中心的经纬度
lon_center = lon_boundary(1:y)+0.05;
lat_center = lat_boundary(1:x)+0.05;
lat_boundary = fliplr(lat_boundary);
lat_center = fliplr(lat_center);           %倒序排列，这样矩阵（1,1）位置处的经纬度就是（lon(1),lat(1)）

%% 确定站点所在位置的行列号
station_lon = station.location(1);                  %读取站点经度
station_lat = station.location(2);                  %读取站点纬度

station_lon_num = 0;
station_lat_num = 0;
for m = 1:length(lat_boundary)-1    %判断这些交点都分别在哪一行上
    if station_lat<lat_boundary(m) && station_lat>=lat_boundary(m+1)
       station_lat_num = m;         %直线所经过的像元点的纬度方向序号确定了
    end
end
for m = 1:length(lon_boundary)-1    %判断这些交点都分别在哪一行上
    if station_lon>=lon_boundary(m) && station_lon<lon_boundary(m+1)
       station_lon_num = m;        %直线所经过的像元点的经度方向序号确定了
    end
end
station.locationnum = [station_lat_num,station_lon_num];

%% 确定路径，并计算平均降水高程梯度
double S;
S_bar = zeros(x,y);
for i = 1:x
    for j = 1:y
       %% 计算站点坐标到目标像元中心点连线的斜率
        k = (lat_center(i)-station_lat)/(lon_center(j)-station_lon);
        b = station_lat-k*station_lon;           %求截距            
        if isequal([i,j] , station.locationnum)  %如果(i,j)是站点所在像元，就跳过，像元值等于站点值 
            continue
        end
       %% 判断直线斜率，以确定用哪种形式的判断路径点方式
        across_lon_num = 0;                    %记录路径线所经过的像元的列序号
        across_lat_num = 0;                    %记录路径线所经过的像元的行序号
        if abs(k)<1                            %此时直线坡度缓，穿过更多经度线，以与竖轴（经线）的交点判断路径
            % 判断需要经过几条经度线
            max_lon = max(lon_center(j),station_lon);
            min_lon = min(lon_center(j),station_lon);
            across_lon = lon_boundary(lon_boundary<max_lon);
            across_lon = across_lon(across_lon>min_lon);  %所需要穿过的经度线  
            across_lon_y = k*across_lon+b;                %穿过纬度线时的交点的纬度坐标
            for m = 1:length(across_lon_y)                %判断这些交点都分别在哪一行上
                for n = 1:x
                    if across_lon_y(m)<lat_boundary(n) && across_lon_y(m)>=lat_boundary(n+1)
                       across_lat_num(m) = n;  %直线所经过的像元点的纬度方向序号确定了
                    end
                end
            end
            %根据穿过的经度线确定像元点的列数
            for m = 1: length(across_lon)
                across_lon_num(m) = find(lon_boundary == across_lon(m));   %确定穿过的纬度线的序号
            end
            across_num = [across_lat_num;across_lon_num];
            across_num = [across_num,[across_num(1,:);across_num(2,:)-1]];
            across_num = unique(across_num','rows');      %删除重复的点
            across_num = across_num';
            across_num(:,across_num(2,:) == 0) = [];
        else                                    %此时直线坡度陡，所以以与横轴的交点判断路径点
            % 判断需要经过几条纬度线
            max_lat = max(lat_center(i),station_lat);
            min_lat = min(lat_center(i),station_lat);
            across_lat = lat_boundary(lat_boundary<max_lat);
            across_lat = across_lat(across_lat>min_lat);  %所需要穿过的纬度线
            across_lat_x = (across_lat-b)/k;              %穿过纬度线时的交点的经度坐标
            %检查点据是否合理：scatter([across_lat_x,lon_center(j),station_lon],[across_lat,lat_center(i),station_lat])
            for m = 1:length(across_lat_x)      %判断这些交点都分别在哪一列上
                for n = 1:y
                    if across_lat_x(m)>=lon_boundary(n) && across_lat_x(m)<lon_boundary(n+1)
                       across_lon_num(m) = n;   %直线所经过的像元点的经度方向序号确定了
                    end
                end
            end
            %根据穿过的纬度线确定像元点的行数
            for m = 1: length(across_lat)
                across_lat_num(m) = find(lat_boundary == across_lat(m));   %确定穿过的纬度线的序号
            end
            across_num = [across_lat_num;across_lon_num];
            across_num = [across_num,[across_num(1,:)-1;across_num(2,:)]]; 
            across_num = unique(across_num','rows');       %删除重复的点
            across_num = across_num';
            across_num(:,across_num(1,:) == 0) = [];
        end
        %此时across_num中存储的是路径上所有像元点的行列号，然而却不是按顺序排下来的（只是以经度或者纬度大小排列）
        %给路径上的像元重排序，顺序从站点位置开始，保证像元是能连起来的
        across_num_reorder = station.locationnum';
        across_num(:,find(sum(across_num == station.locationnum')==2)) = [];
        len = length(across_num(1,:));
        for num1 = 1:len
            col = find(sum(abs(across_num_reorder(:,num1)-across_num))==1);%之所以要求相减为1就是保证能连起来（比如(3,3)和(2,2)就是断开的）
            across_num_reorder(:,num1+1) = across_num(:,col);
            across_num(:,col) = [];
        end
        across_num = across_num_reorder;
    end
end

此时的across_num中存储的就是路径上像元点的行列号了，比如：
站点在（28,35），目标像元在（34,40），则路径为
Finish!