MATLAB计算多景tif和h5文件的斜率以及给定方程的参数

有时候需要我们自己计算多景tif的斜率以及给定方程的参数，这时候用的是MATLAB进行处理，代码如下：

clc
close all
clear
%读取时间
new_time = xlsread('time_04_06.xls');
[row,col]=size(new_time);
time=zeros(row,1);
for i=1:row
    time(i)=new_time(i,1)-2004+(datenum(new_time(i,:))-datenum(new_time(i,1),1,1)+1)/(datenum(new_time(i,1),12,31)-datenum(new_time(i,1),1,1)+1);
end
%%选择目标文件夹,注意文件夹里的数据名要按顺序排列
%特别是进行时间序列分析时
path= 'D:\2004_2006'; %选择文件夹位置，获取路径
fileList=dir(fullfile(path,'*.tif')); %获取文件夹下所有tif文件
n=length(fileList); %获取文件个数

%%创建三维矩阵，并获取参考矩阵和投影信息
%h=waitbar(0,'Running','name','3DMAT'); %进度条1
for i=1:n
    filePath=fullfile(path,fileList(i).name);
    data(:,:,i)=geotiffread(filePath);
    info=geotiffinfo(filePath);
    R=info.RefMatrix; %参考矩阵
    proj=info.GeoTIFFTags.GeoKeyDirectoryTag; %投影信息
    str=['Completed:','',num2str(round(i*100/n)),'%'];
    %waitbar(i/n,h,str); %进度条2
end
%delete(h); %删除进度条

%%进行所需各项计算
[m,n,p]=size(data); %行(range)、列(azimuth)、页(tims)
data2 = reshape(data,13582908,12);
slope=zeros(13582908,1); %形变速率
amp=zeros(13582908,1); %周期幅度
time=[ones(size(time)),time];

for i=1:13582908
        x=time(:,2);
        y=data2(i,:)';
        %p=polyfit(x,y,1);
        %slope(i,1)=p(1);
        %fprintf('i=%d\n',i)
        %%%%%%%%%给定的方程%%%%%%%%%%%%%%%
        f = fittype('a * sin(6.28 * t) + b * cos(6.28 * t) + v * t + c', 'independent','t','coefficients',{'a','b','v','c'});
        cfun = fit(x, y, f);
        coef = coeffvalues(cfun);
        slope(i,1) = coef(3);
        amp(i,1) = sqrt(coef(1)^2+coef(2)^2);
        fprintf('i=%d\n',i)
end
slope2 = reshape(slope,3687,3684);
amp2 = reshape(amp,3687,3684);

%%将计算结果写入到新的TIFF影像数据中
geotiffwrite('slope.tif', slope2, R,'GeoKeyDirectoryTag', proj);
geotiffwrite('amp.tif', amp2, R,'GeoKeyDirectoryTag', proj); %使输出的投影与原始一致

同样的也可以利用MATLAB计算HDF5文件里面的斜率或给定的参数：

clc
close all
clear
new_time = xlsread('time_ALL.xls');
[row,col]=size(new_time);
time=zeros(row,1);
for i=1:row
    time(i)=new_time(i,1)-2014+(datenum(new_time(i,:))-datenum(new_time(i,1),1,1)+1)/(datenum(new_time(i,1),12,31)-datenum(new_time(i,1),1,1)+1);
end

def_dir='PARAMS.h5';
def_data=hdf5read(def_dir,'rawts');
[range,azimuth,tims]=size(def_data);
lat=[31.7194,35.5459];
lon=[90.9887,93.7422];

%生成经纬度文件 lower left
originX = 90.9887; 
originY = 31.7194; 
pixelSize = 0.00027777777777778;
longitude=linspace(originX,pixelSize*range+originX,range);
latitude=linspace(pixelSize*azimuth+originY,originY,azimuth);

data1 = permute(def_data,[2 1 3]);
data = reshape(data1,136571401,189);

slope=zeros(136571401,1);
time=[ones(size(time)),time];
%fprintf('time is %d\n',time)

for i=1:136571401
        y=data(i,:)';
        x=time(:,2);
        p=polyfit(x(2:end),y(2:end),1);
        slope(i,1)=p(1);
        fprintf('i=%d\n',i)
end
slo2pe = reshape(slope,13777,9913);
slo3pe=flipud(slo2pe);

Refference = georasterref('RasterSize',[azimuth range],'Latlim',lat,'Lonlim',lon);%给文件一个地理坐标系WGS1984
outname='velocity.tif';
geotiffwrite(outname,slo3pe,Refference);