本文主要提供一套處理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ù)的行列號是反著的了袁，此時就需要對行列號進行變換，采用rot90函數(shù),并進行查看
data1=rot90(precipitation)
imshow(data1)
結(jié)果表明已經(jīng)反轉(zhuǎn)過來了湿颅，但是這個時候還要考慮到是不是需要南北轉(zhuǎn)換一下载绿，從圖中是無法分辨出來的，因此需要對反轉(zhuǎn)和未反轉(zhuǎn)的結(jié)果分別用先驗知識進行判斷油航。在判斷前崭庸，需要做個樣例出來，定義好投影信息谊囚，以便輸出其他圖像怕享。樣例制作過程如下
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)換教程函筋，里面有詳細說明。
一個月的數(shù)據(jù)不容易進行判斷奠伪，因此需要一年的數(shù)據(jù)進行判斷跌帐，首先將反轉(zhuǎn)和未反轉(zhuǎn)的數(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
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ū)進行裁剪來進行判斷哪種是正確的降水。本文采用黃河流域來進行裁剪谨敛，結(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ù)集的提取過程了。
更多需求佣盒，請查看個人介紹
