利用MTLAB代码批量计算每一年中的某一月的总初级生产力（GPP）

最新推荐文章于 2024-03-03 08:00:00 发布

地理学码农

最新推荐文章于 2024-03-03 08:00:00 发布

阅读量360

点赞数 1

文章标签：前端数据库 matlab

本文链接：https://blog.csdn.net/qq_44996109/article/details/133634561

版权

一、实验背景

1、利用栅格数据计算平面上每一个像元的每一年中的某一月的总初级生产力（GPP）时，不能简单的求和，需要考虑其时间分辨率。以现有的8天时间分辨率的GPP数据为例（本文实验数据），每一期数据代表的是8天内GPP的平均值，因此1月的GPP应该是将一月的所有期数据乘以8后再加和。

2、现有的GPP数据文件命名大部分使用年份和儒略日的组合来表示，例如：2001225.tif表示该数据是2001年的第225天后的8天GPP平均值，因此如果需要计算某一月的总GPP，我们需要从这一年中选出这一月的所有期数据进行计算，如果手动挑出再计算会造成巨大的工作量，尤其是一般会有很多年需要计算。

二、设计实验

1、实验数据：2001~2020年中国东北部8天时间分辨率，0.05°空间分辨率的GPP数据

2、软件：MATLAB2016

3、实验思路

（1）逐年挑选出8月的所有期数据后放在三维矩阵中

（2）批量计算每一年8月的GPP

4、代码实现

话不多说，上代码!

clear;
clc;
close all;

[a,R]=geotiffread('D:\text\data-text3_new\GOSIF_GPP\clip\GOSIF_GPP_2001001_Mean.tif');    %定义投影
info = geotiffinfo('D:\text\data-text3_new\GOSIF_GPP\clip\GOSIF_GPP_2001001_Mean.tif');
[m,n]=size(a);
temp = zeros(m,n);
result=zeros(m,n);
inputfile='D:\text\data-text3_new\GOSIF_GPP\clip/';    
outputfile='D:\text\data-text3_new\GOSIF_GPP\AUG/';

for year=2001:2021    %每一年的影像期数
    i=0;        
    ID=0;
    month2=zeros(1,46);        %存放每一个儒略日对应的月份
    for j=1:8:361
        ID=ID+1;
       [year2,month,day] = jday2date(year,j);
       month2(1,ID)=month;
       NDVI=zeros(m,n,4);
       if month2(1,ID)==8            %找出8月的数据，只计算每一年8月的GPP，不符合if条件的不做运算
            i=i+1;       
            counter=(ID-1)*8+1;
            NDVI(:,:,i)=imread(strcat(inputfile,'GOSIF_GPP_',num2str(year),num2str(counter,'%03d'),'_Mean.tif')); 
            temp=NDVI(:,:,i);
            temp(temp<0) = 0;      %python处理后的数据NoData值在matlab中显示为负无穷大，需要去除黑边
            temp(temp>6000)=0;   %GOSIF数据的水体和无植被区域覆盖为大于6000的值
            temp=temp*0.001;  %GOSIF GPP需要乘以尺度因子
            result = temp*8 + result;    %以2001年为例，将8月的四幅矩阵乘以8天后累加放在result中
       end      
    end
     geotiffwrite(strcat(outputfile,num2str(year),'aug.tif'),result, R,'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
     result=zeros(m,n);        %格式化处理，清空矩阵后准备存放下一年的三维数据
end

5、注意事项

（1）本代码中将儒略日转换为月份的部分参考自儒略日(Julian Day)与日期相互转化的python/matlab函数_matlab 儒略日转换为年月日-CSDN博客y

因此，需要将所有函数文件和脚本文件放在同一工作目录下以保证19行不报错，如图：