MOD13Q1数据最大合成NDVI数据

MOD13Q1数据的具体介绍可以到官方网站上进行阅读，这里就不再赘述，这里就直接上代码了。

一、处理过程

1、根据日期（天数）判断是那个月份。
2、根据月份将数据按照季度、年度进行NDVI最大值合成。

二、代码实现（MATLAB）

daycal.m

function [day]=daycal(year,month)
%根据某年某月推算出该月份的天数
num=[1,3,5,7,8,10,12];
if mod(year,4)==0   %&& mod(year,100) %判断year是否为闰年，因为闰年的2月是29天
    if month==2
        day=29;
    elseif sum(ismember(num,month))==1
        day=31;
    else
        day=30;
    end
else%非闰年的情况
    if month==2
        day=28;
    elseif sum(ismember(num,month))==1
        day=31;
    else
        day=30;
    end
end 
end

DayToMonth.m

function Month = DayToMonth(date)
%DAYTOMONTH 此处显示有关此函数的摘要
%根据天数确定月数
year = floor(date / 1000);
days = mod(date,1000);
months = MonthdaysInYear(year);
for i=1:12
   if days <= months(1,i)
       Month = i;
       break;
   end
end
end

MonthdaysInYear.m

function  months = MonthdaysInYear(year)
% 此处显示有关此函数的摘要
% 根据年份获取该年份所有月的天数
    months = zeros(1,12);
    months(1,1) = 31; 
    for i=2:12
        months(1,i) = months(1,i-1)+daycal(year,i);
    end
end

NDVIMerge.m

function ndvi = NDVIMerge(img1,img2)
%NDVIMERGE 此处显示有关此函数的摘要
%对输入数据进行NDVI最大合成
shape = size(img1);
row = shape(1,1);
col = shape(1,2);
ndvi = img1;
for i=1:row
    for j=1:col
        if ndvi(i,j) < img2(i,j)
            ndvi(i,j) = img2(i,j);
        end
    end
end
end

main.m

clc
clear
close all

%获取数据文件
dirPath=uigetdir();   %选择要进行计算文件路径

if isempty(dirPath) || length(dirPath) == 1
    fprintf("未选择目录！\n");
    return;
end
outputPath = 'C:\Users\23547\Desktop\others\yingying\result';
outputPath = [outputPath,'\'];
files=dir(dirPath);
fileNames={files.name};
pathName=[dirPath,'\'];

quarter=["春","夏","秋","冬"];
quarterNum = 1;
fprintf("开始计算...\n");
fullPath = [pathName,fileNames{3}];
fprintf("读取%s！\n",fileNames{3});
[Data,R]=geotiffread(fullPath);   %加载数据
info = geotiffinfo(fullPath);  % 读取tif数据的地理信息，为后面导出为tif数据提供地理信息
strArr = strsplit(fileNames{3},".");
dateRecord = strArr{2};
isNumber = isstrprop(dateRecord,'digit');
date=str2num(dateRecord(isNumber));
preyear = floor(date / 1000);
quarterResult = Data;
allResult = Data;
for m=4:size(fileNames,2)
    fullPath = [pathName,fileNames{m}];
    fprintf("读取%s！\n",fileNames{m});
    [Data,R]=geotiffread(fullPath);   %加载数据
    strArr = strsplit(fileNames{m},".");
    dateRecord = strArr{2};
    isNumber = isstrprop(dateRecord,'digit');
    date=str2num(dateRecord(isNumber));
    year = floor(date / 1000);
    month = DayToMonth(date);
    if year-preyear>0
        filename = [outputPath,num2str(floor(preyear)),num2str(quarter(quarterNum)),...
                strArr{3},'.tif'];  %存储位置和名字
        quarterResult=double(quarterResult)/10000;
        geotiffwrite(filename,quarterResult, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);  %导出

        filename = [outputPath,num2str(floor(preyear)),...
                strArr{3},'all.tif'];  %存储位置和名字
        allResult=double(allResult)/10000;
       geotiffwrite(filename,allResult, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);  %导出
       quarterResult = Data;
       allResult = Data;
       quarterNum = 1;
       preyear = year;
       continue;
    elseif month > 3*quarterNum
       filename = [outputPath,num2str(year),num2str(quarter(quarterNum)),...
           strArr{3},'.tif'];  %存储位置和名字
       quarterResult=double(quarterResult)/10000;
       geotiffwrite(filename,quarterResult, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);  %导出
       quarterResult = Data;
       allResult=NDVIMerge(allResult,Data);
       quarterNum=quarterNum+1;
       continue;
    end
    quarterResult=NDVIMerge(quarterResult,Data);
    allResult=NDVIMerge(allResult,Data);
end
filename = [outputPath,num2str(floor(year)),num2str(quarter(quarterNum)),...
    strArr{3},'.tif'];  %存储位置和名字
quarterResult=double(quarterResult)/10000;
geotiffwrite(filename,quarterResult, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);  %导出

filename = [outputPath,num2str(floor(year)),...
    strArr{3},'all.tif'];  %存储位置和名字
allResult=double(allResult)/10000;
geotiffwrite(filename,allResult, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);  %导出
fprintf("计算成功！\n");

三、实现效果