基于Matlab的TRMM3B43数据处理的思维过程与技术流程

本文主要提供一套处理TRMM3B43数据的思维过程与技术流程，力求能够让读者在处理其他国际通用数据时也能够采用类似的方法来解决。
首先我们来看一下TRMM3B43数据的数据格式，是hdf格式的，相当于nc，tif文件还是较难读入的，不同的hdf文件都有着各自的文件结构，那么如何来读取呢。首先我们利用matlab的导入数据功能来导入一个数据，并在导入过程中查看参数，可以看到单位是mm/hr，即毫米每小时。

导入数据.jpg
单位.jpg

选择好其中的一个数据，加载可以得到如下窗口，从下图的红色区域可以得到导入该hdf数据的代码格式为
hdfread('D:\3B43.19980101.7.HDF', '/Grid/precipitation', 'Index', {[1 1],[1 1],[1440 400]});其中'D:\3B43.19980101.7.HDF'便是文件的名称，可以用作后续的循环迭代。

trmm.jpg

查看其它的参数选项，可以得到TRMM3B43数据的空间范围是上下50度，左右180度，获取范围后可以有针对的设置投影范围等。

trmm范围.jpg

导入并查看数据，查看数据一般用imshow函数

precipitation = hdfread('D:\3B43.19980101.7.HDF', '/Grid/precipitation', 'Index', {[1  1],[1  1],[1440   400]});
imshow(precipitation)

查看数据.jpg

通过查看数据发现数据的行列号是反着的，此时就需要对行列号进行变换，采用rot90函数,并进行查看

data1=rot90(precipitation)
imshow(data1)

结果表明已经反转过来了，但是这个时候还要考虑到是不是需要南北转换一下，从图中是无法分辨出来的，因此需要对反转和未反转的结果分别用先验知识进行判断。在判断前，需要做个样例出来，定义好投影信息，以便输出其他图像。样例制作过程如下

precipitation = hdfread('D:\3B43.19980101.7.HDF', '/Grid/precipitation', 'Index', {[1  1],[1  1],[1440   400]});
data1=rot90(precipitation);
dlmwrite('H:\Global\3b43\example.txt',data1,'\t',1,1);

打开example.txt文件添加投影信息，如下所示，后续的在arcgis中去定义这个投影信息的过程见本博客的PDSI的转换教程，里面有详细说明。

案例.jpg

一个月的数据不容易进行判断，因此需要一年的数据进行判断，首先将反转和未反转的数据分别进行输出，代码如下：

[~,R]=geotiffread('H:\Global\3b43\example.tif');
info=geotiffinfo('H:\Global\3b43\example.tif');
mon_length_pingnian=[31 28 31 30 31 30 31 31 30 31 30 31];
mon_length_runnian=[31 29 31 30 31 30 31 31 30 31 30 31];
for year=1998
    datasum1=0;
    datasum=0;
    if mod(year,4)==0;
        cd=mon_length_runnian;
    else
       cd= mon_length_pingnian;
    end
    for mon=1:12
        if mon<10
            filename=['H:\Global\3b43\3B43.',int2str(year),'0',int2str(mon),'01.7.HDF'];
        else
            filename=['H:\Global\3b43\3B43.',int2str(year),int2str(mon),'01.7.HDF'];
        end
        data=hdfread(filename, '/Grid/precipitation', 'Index', {[1  1],[1  1],[1440   400]});
        data=rot90(data);
        data=data.*24*cd(mon);
        datasum=datasum+data;
        data1=flipud(data);
        data1=data1.*24.*cd(mon);
        datasum1=datasum1+data1;
    end
    geotiffwrite('H:\Global\3b43\未反转的precp.tif',datasum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
    geotiffwrite('H:\Global\3b43\反转的precp.tif',datasum1,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end

采用熟悉的研究区进行裁剪来进行判断哪种是正确的降水。本文采用黄河流域来进行裁剪，结果分别如下所示,结果表明是未经过上下反转的。

判断结果.jpg

因此可用以下代码来提取全球TRMM3B43每年的月和年数据集

[~,R]=geotiffread('H:\Global\3b43\example.tif');
info=geotiffinfo('H:\Global\3b43\example.tif');
mon_length_pingnian=[31 28 31 30 31 30 31 31 30 31 30 31];
mon_length_runnian=[31 29 31 30 31 30 31 31 30 31 30 31];
for year=1998:2017
    datasum=0;
    if mod(year,4)==0;
        cd=mon_length_runnian;
    else
       cd= mon_length_pingnian;
    end
    for mon=1:12
        if mon<10
            filename=['H:\Global\3b43\3B43.',int2str(year),'0',int2str(mon),'01.7.HDF'];
        else
            filename=['H:\Global\3b43\3B43.',int2str(year),int2str(mon),'01.7.HDF'];
        end
        data=hdfread(filename, '/Grid/precipitation', 'Index', {[1  1],[1  1],[1440   400]});
        data=rot90(data);
        data=data.*24*cd(mon);
        filename_mon=['H:\Global\3b43\TRMM3B43_precp_',int2str(year),'_',int2str(mon),'月降水.tif'];
        geotiffwrite(filename_mon,data,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
        datasum=datasum+data;
    end
    filename=['H:\Global\3b43\TRMM3B43_year_precp_',int2str(year),'年降水.tif'];
    geotiffwrite(filename,datasum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end

提取出来的结果如下图所示

结果文件.jpg

通过上述过程就能够实现一个完整的全球普通数据集的提取过程了。