利用MOD13Q1 NDVI数据生成月最大值

基于matlab的程序

%求每一年12月的最大值

clc

clear

[a,R]=geotiffread('E:\NDVI_WINTER\10-20CLIPNDVI\TPMOD13Q1.A2000337.250m_16_days_NDVI.tif');%定义投影，后续图像输出

info=geotiffinfo('E:\NDVI_WINTER\10-20CLIPNDVI\TPMOD13Q1.A2000337.250m_16_days_NDVI.tif');

[m,n]=size(a);

begin_year=2000;

end_year=2019;

qs=2;

%总共2期数据

for year=begin_year:end_year

    year

    ndvisum=zeros(m*n,qs);

    k1=1;

    for k=337:16:353

        filename=['E:\NDVI_WINTER\10-20CLIPNDVI\TPMOD13Q1.A',int2str(year),int2str(k),'.250m_16_days_NDVI.tif']; 

        ndvi=double(importdata(filename));

        ndvi=reshape(ndvi,m*n,1);

        ndvisum(:,k1)=ndvi;

        k1=k1+1;

        ndvisum(ndvisum==-3000 | ndvisum==-32768)=NaN;%为空值赋值

    end

    max=zeros(m*n,1);

    a337=ndvisum(:,1);

    a353=ndvisum(:,2);

    %length(find(max==0))

    logt= a337 >=a353;

    max = a337 .*logt + a353 .*(1 - logt);

    max=reshape(max,m,n);

    filename=['E:\NDVI_WINTER\月最大值\空值\TP',int2str(year),int2str(12),'.month_maxNDVI2.tif'];

    geotiffwrite(filename,max,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

end