-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathTRT_process.m
164 lines (144 loc) · 5.64 KB
/
TRT_process.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
% Yiwen Mei (ymei2@gmu.edu)
% CEIE, George Mason University
% Last update: 3/2/2018
%% Functionality
% This function is used to process the TMPA 3B42-RT satellite precipitation
% product (Huffuman et al. 2007). It includes several functionalities:
% 1)unzip the TMPA 3B42-RT record (3B42-RT.2000030100.7R2.bin.gz);
% 2)read and crop the record based on an given lat/lon box;
% 3)output the cropped record as .asc file (optional);
% 4)project the record in .asc to another projection and output the record
% in GTiff (optional).
%% Input
% infname: full name with path of the input TMPA 3B42RT file (e.g.,
% G:\TRMM\3B42RT\2000\3B42RT.2000030100.7R2.bin.gz);
% fld : 0 or 1 stands for the real-time and CCA field;
% workpth: path to store the unzipped record;
% xl/xr : left/right longitude of the boundary (xl/xr can have either 1 or 2 number(s)
% where the first one represents the longitude and must be in the range of
% [-180 180]; the second one is boundary in the projected coordinate unit);
% yb/yt : bottom/top latitude of the boundary (yb/yt can have either 1 or 2 number(s)
% where the first one represents the latitude and must be in the range of
% [60 -60]; the second one is boundary in the projected coordinate unit);
% outpth : path to store the .asc and, if have, the .tif files (set it to "[]" if
% no need to output record in .asc format);
% out_pj : output coordinate system (e.g., EPSG:102009; set it to "[]" if no
% reprojection is required);
% rs : x and y resolution of the projected image (set it to "[]" if no
% reprojection is required).
%% Output
% p / p1 : cropped precipitation map in original/new projection (the
% orientation follows the human reading convention);
% 3B42RTyyyymmddhh.asc: precipitaiton map in original projection and resolution
% outputted to outpth as .asc file;
% 3B42RTpyyyymmddhh.tif: precipitaiton map in new projection, resolution and extend
% outputted to outpth as .tif file;
%% Additional note
% 1)Please make sure to have GDAL installed and the outpth set if you want to
% reproject the data into other coordinate system.
% 2)If reprojection is not required (i.e., out_pj is "[]") but record outputted
% as .asc is wanted, out_path is required to set.
% 3)The scale factor and no-data value of TMPA 3B42RT are preseved in the .asc
% and .tif record.
% 4)No scale factor is preserved in p and p1 and no-data value is replaced by NaN.
function [p,p1]=TRT_process(infname,fld,workpth,xl,xr,yb,yt,outpth,out_pj,rs)
% Lat/lon grids and other info of TMPA 3B42RT
rs_lon=360/1440;
rs_lat=120/480;
Lon=0:rs_lon:360;
Lat=60:-rs_lat:-60;
ndv=-31999; % no-data value
scf=100; % Scale factor
% Index of interested domain
if xl(1)<0 % Convert longitude to the range of [0 360];
xl(1)=xl(1)+360;
end
if xr(1)<0
xr(1)=xr(1)+360;
end
cl=find(xl(1)-Lon>=0,1,'last'); % left column
cr=find(xr(1)-Lon<=0,1,'first')-1; % right column
rt=find(yt(1)-Lat<=0,1,'last'); % top row
rb=find(yb(1)-Lat>=0,1,'first')-1; % bottom row
nr=length(rt:rb); % number of row
nc=length(cl:cr); % number of column
xll=(cl-1)*rs_lon; % longitude of lower left corner
yll=60-rb*rs_lat; % latitude of lower left corner
% unzip and read the input
% system(sprintf('7z e "%s" -o"%s" * -r',in_fname,workpth));
gunzip(infname,workpth); % unzip
[~,nm,~]=fileparts(infname);
uz_fn=[workpth nm];
fid=fopen(uz_fn,'r'); % read
p=fread(fid,[1440 1681],'integer*2','b');
if fld==0 % real-time
p=p(:,1202:1681)'; % original upper-left corner is (N,W). Flip to (W,N).
else % CCA
p=p(:,2:481)'; % original upper-left corner is (N,W). Flip to (W,N).
end
fclose(fid);
p=p(rt:rb,cl:cr); % crop
delete(uz_fn);
p1=[];
if ~isempty(out_pj)
% Creat the asc files
ds=nm(8:9+8); % pay careful attention to this
if fld==0
name=[outpth '3B42RT' ds '.asc'];
else
name=[outpth '3B42CCA' ds '.asc'];
end
fid=fopen(name,'w');
fprintf(fid,'%s\n%s\n%s\n%s\n%s\n%s\n',['ncols ' num2str(nc)],['nrows '...
num2str(nr)],['xllcorner ' num2str(xll,8)],['yllcorner ' num2str(yll,8)],...
['cellsize ' num2str(rs_lat)],['NODATA_value ' num2str(ndv)]);
fclose(fid);
dlmwrite(name,p,'delimiter',' ','-append'); % output .asc
% Project to a new coordinate
fun='gdalwarp -overwrite -of GTiff -r bilinear '; % GDAL function
pr1='-s_srs wgs84 '; % Projection of original record
pr2=['-t_srs ' out_pj ' '];
pr3=[];
if ~isempty(rs)
pr3=sprintf('-tr %i %i ',rs(1),rs(2));
end
pr4=[];
if length(xl)==2
pr4=sprintf('-te %i %i %i %i ',xl(2),yb(2),xr(2),yt(2));
end
par=[pr1 pr2 pr3 pr4];
inv=['"' name '" '];
if fld==0
ouv=['"' outpth '3B42RTp' ds '.tif"'];
else
ouv=['"' outpth '3B42CCAp' ds '.tif"'];
end
system([fun par inv ouv]); % project
delete(name);
if fld==0
p1=double(imread([outpth '3B42RTp' ds '.tif']));
else
p1=double(imread([outpth '3B42CCAp' ds '.tif']));
end
p1(p1==ndv)=NaN;
p1=p1/scf;
else
if ~isempty(outpth)
% Creat the asc files
ds=nm(8:end-8);
if fld==0
name=[outpth '3B42RT' ds '.asc'];
else
name=[outpth '3B42CCA' ds '.asc'];
end
fid=fopen(name,'w');
fprintf(fid,'%s\n%s\n%s\n%s\n%s\n%s\n',['ncols ' num2str(nc)],['nrows '...
num2str(nr)],['xllcorner ' num2str(xll,8)],['yllcorner ' num2str(yll,8)],...
['cellsize ' num2str(rs_lat)],['NODATA_value ' num2str(ndv)]);
fclose(fid);
dlmwrite(name,p,'delimiter',' ','-append');
end
end
p(p==ndv)=NaN;
p=p/scf;
end