Matlab计算多幅遥感影像交集面积占比

Posted 2021-03-06 icydengyw

一个文件夹中有同一时间不同的遥感影像，有LST NDVI、反射率等，而想要知道这些像元数量相同，填充值和背景值均为nan的影像的交集部分占该影像研究区（已经经过研究区矢量裁剪）的比例，以及不同类别遥感影像的覆盖面积比例，主要目的是想以多种数据交集覆盖面积比例来判断改天是否适合数据反演

1 函数

函数基本结构

function[输出形参表: output1, ...,outptn] = 函数名(输入形参表: input1, ... , inputn)

注释说明部分
函数体代码部分

end

函数名：需要和m文件名保持一致

注1：函数中遇return语句时，将退出函数体，此函数调用结束；

注2：函数体里面也可以定义一个或几个函数，称为子函数；注意：子函数只能存在于主函数体内，不独立存在；子函数在主函数体内的位置可以任意，不影响使用；子函数只能被主函数以及其他位于同一主函数体下的子函数调用，但子函数“句柄”例外；

注3：在调用函数时，Matlab用两个永久变量nargin和nargout分别记录调用该函数时的输入实参和输出实参的个数。只要在函数文件中包含这两个变量，就可以准确地知道该函数文件被调用时的输入输出参数个数，从而决定函数如何进行处理。

内联函数

可以让用户编写简单的函数而不需要创建M文件。语法格式：

f=inline(‘函数表达式‘, ‘变量1‘, ‘变量2‘, ……)

调用方式：y=f(实参列表)

注意：实参列表顺序应与inline()定义中形参列表的顺序一致。

这种函数定义方式是将它作为一个内部函数调用。好处是，它是基于Matlab的数值运算内核的，所以它的运算速度较快，程序效率更高。缺点是，该方法只能对数值进行代入，不支持符号代入，且对定义后的函数不能进行求导等符号运算。

f=inline(‘x^2+y‘,‘x‘,‘y‘)

z=f(2,3)

 匿名函数　

 匿名函数具有内联的所有优点，并且效率比更高。匿名函数的主要功能是：

（1）可以代替“将函数编写为单独的m-文件”；

（2）可以实现符号函数的赋值运算；

（3）很方便地对含参变量函数进行操作。

基本格式：

f=@(参数1, 参数2, …) 函数表达式

其中，@为句柄操作符，可以定义指向函数的句柄，函数句柄，可以理解成函数的“代号”，适用于函数名比较长，用“代号”代替函数操作更方便。

f=@(x,y)x^2+y^2;
f(2,3)  %结果为13

a=1:5;
b=5:-1:1;
c=0.1:0.1:0.5;
g=@(x,y)x.^2+y.^2+c;
g(a,b) %结果为向量

函数参数可变

function re=add_multi(a,b,c,d)
% nargin 表示参数个数
if nargin==4
    re=a+b+c+d;
elseif nargin==3
    re=a+b+c;
elseif nargin==2
    re=a+b;
else
    error(‘输入有错‘);
end

2 多维数组

前两个维度就像一个矩阵，而第三个维度表示元素的页数或张数。

3 实现

原始文件：

clc
clear all

inpath=‘E:1SMAPDATAMOD11A1_testCategoryEqual‘;%输入文件夹，内含不同时间子文件夹
outpath=‘E:1SMAPDATAMOD11A1_testCategoryareaxls‘;%输出的面积比例xls表的位置
files = dir(inpath);
oldName = cell(length(files)-2,1);%含有所有子文件夹名字的元胞数组 前两个无效 
for i = 3:length(files)
    oldName{i-2} = files(i).name;%把文件夹名字提取出来重新排列
end
table={};%存放面积比例的元胞数组
table{1,1}=‘time‘;
for i = 1:length(oldName)%遍历每个子文件夹
    subFolder=[inpath ,‘‘,oldName{i}];
    fileFolder=fullfile(subFolder);%列出子文件夹中所有的tiff文件   
    perc=areaPerc(fileFolder,208300);%遍历计算出每个子文件夹中所有变量的面积覆盖比例 areaPerc函数返回某一时间面积比例的元胞数组   
    table{i+1,1}=oldName{i};%记录时间 
    for j=1:size(perc,2)  %写入面积比例数据             
        table{i+1,j+1}=perc{2,j};
    end
    if i==1 %第一次循环写入表头字段
        for j=1:size(perc,2)        
        table{1,j+1}=perc{1,j};        
        end
    end
    
end
writetable(cell2table(table),[outpath ,‘‘,‘areaPerc.xls‘]);%将元胞数组写入xls文件

areaPerc函数：

function perc = areaPerc( imagepath,jjjnum )
%areaPerc Summary of this function goes here
%imagepath是同一时间处的文件夹 里面有不同种类的影像
fileFolder=fullfile(imagepath);%%图像名称与路径 某PM2_5一文件夹下
dirOutput=dir(fullfile(fileFolder,‘*.tif‘)); %一文件夹下所有 *.tif 文件
fileNames={dirOutput.name}‘; %名字 为cell
table={};%数据种类字段

for i=1:length(fileNames)
filepath=[imagepath,‘‘,fileNames{i}];
[A,R]=geotiffread(filepath);%注意这里读取的数据矩阵A是 ‘uint16‘格式,这个格式不能为负数或者是nan
mid=strsplit(fileNames{i},‘.‘);
table{1,i}=mid{2};%有待改进
table{2,i}=(numel(A)-numel(find(isnan(A))))/jjjnum;%某一个tiff文件中的有效值比例
A=double(A);
B(:,:,i) = A;%创建多维数组 
end
%对三维数组B进行遍历 
[m,n]=size(A);
valid=0;
for i=1:m
    for j=1:n
    c=B(i,j,:);     
    if numel(find(isnan(c)))==0 %每一维度的第m行第n列均没有无效值
    valid = valid + 1;
    end
    end
end
allprc=valid/jjjnum;
table{1,length(fileNames)+1}=‘intersect‘;
table{2,length(fileNames)+1}=allprc;
perc=table;
end

结果为：

 后续可以根据intersect的比例进行判别是否比较适合进行反演了。

参考：

https://zhuanlan.zhihu.com/p/25530879

https://ww2.mathworks.cn/help/matlab/ref/struct.html  结构体  https://blog.csdn.net/jxlijunhao/article/details/18048359

https://ww2.mathworks.cn/help/matlab/math/multidimensional-arrays.html 多维数组