本文主要提供一套處理TRMM3B43數(shù)據(jù)的思維過程與技術(shù)流程，力求能夠讓讀者在處理其他國際通用數(shù)據(jù)時也能夠采用類似的方法來解決。
首先我們來看一下TRMM3B43數(shù)據(jù)的數(shù)據(jù)格式，是hdf格式的，相當(dāng)于nc，tif文件還是較難讀入的，不同的hdf文件都有著各自的文件結(jié)構(gòu)，那么如何來讀取呢。首先我們利用matlab的導(dǎo)入數(shù)據(jù)功能來導(dǎo)入一個數(shù)據(jù)，并在導(dǎo)入過程中查看參數(shù)，可以看到單位是mm/hr，即毫米每小時。
選擇好其中的一個數(shù)據(jù)，加載可以得到如下窗口，從下圖的紅色區(qū)域可以得到導(dǎo)入該hdf數(shù)據(jù)的代碼格式為
hdfread('D:\3B43.19980101.7.HDF', '/Grid/precipitation', 'Index', {[1 1],[1 1],[1440 400]});其中'D:\3B43.19980101.7.HDF'便是文件的名稱，可以用作后續(xù)的循環(huán)迭代。
查看其它的參數(shù)選項，可以得到TRMM3B43數(shù)據(jù)的空間范圍是上下50度，左右180度，獲取范圍后可以有針對的設(shè)置投影范圍等。
導(dǎo)入并查看數(shù)據(jù)，查看數(shù)據(jù)一般用imshow函數(shù)
precipitation = hdfread('D:\3B43.19980101.7.HDF', '/Grid/precipitation', 'Index', {[1 1],[1 1],[1440 400]});
imshow(precipitation)
通過查看數(shù)據(jù)發(fā)現(xiàn)數(shù)據(jù)的行列號是反著的，此時就需要對行列號進(jìn)行變換，采用rot90函數(shù),并進(jìn)行查看
data1=rot90(precipitation)
imshow(data1)
結(jié)果表明已經(jīng)反轉(zhuǎn)過來了，但是這個時候還要考慮到是不是需要南北轉(zhuǎn)換一下，從圖中是無法分辨出來的，因此需要對反轉(zhuǎn)和未反轉(zhuǎn)的結(jié)果分別用先驗知識進(jìn)行判斷。在判斷前，需要做個樣例出來，定義好投影信息，以便輸出其他圖像。樣例制作過程如下
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文件添加投影信息，如下所示，后續(xù)的在arcgis中去定義這個投影信息的過程見本博客的PDSI的轉(zhuǎn)換教程，里面有詳細(xì)說明。
一個月的數(shù)據(jù)不容易進(jìn)行判斷，因此需要一年的數(shù)據(jù)進(jìn)行判斷，首先將反轉(zhuǎn)和未反轉(zhuǎn)的數(shù)據(jù)分別進(jìn)行輸出，代碼如下：
% @author yinlichang3064@163.com
[~,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\未反轉(zhuǎn)的precp.tif',datasum,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite('H:\Global\3b43\反轉(zhuǎn)的precp.tif',datasum1,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
end
采用熟悉的研究區(qū)進(jìn)行裁剪來進(jìn)行判斷哪種是正確的降水。本文采用黃河流域來進(jìn)行裁剪，結(jié)果分別如下所示,結(jié)果表明是未經(jīng)過上下反轉(zhuǎn)的。
因此可用以下代碼來提取全球TRMM3B43每年的月和年數(shù)據(jù)集
%@author yinlichang3064@163.com
[~,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
提取出來的結(jié)果如下圖所示
通過上述過程就能夠?qū)崿F(xiàn)一個完整的全球普通數(shù)據(jù)集的提取過程了。
更多需求，請查看個人介紹
