forked from flozabel/Teddy
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwfde5_tiles.m
100 lines (86 loc) · 3.35 KB
/
wfde5_tiles.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
function wfde5_all=wfde5_tiles(upperyy,loweryy,leftxx,rightxx,tilesize,years_wfde5,startyear_wfde5,wfde5dir,parallel_cpus,mask)
%open parallel pool
delete(gcp('nocreate'));
parpool('Processes',min(parallel_cpus,12));
res=0.5;
parameters_wfde5={'Tair','Qair','SWdown','LWdown','PSurf','Wind','Rainf'};
for uppery=upperyy:-tilesize:loweryy
lowery=uppery-tilesize;
[lat,lon]=calc_coordinates_global_land(uppery,lowery,leftxx,rightxx,res,mask);
disp(['processing tile ',num2str(uppery),'° - ',num2str(lowery),'°']);
%convert lat/lon to row,col
yrow=floor((90-lat)/res)+1;
xcol=floor((180+lon)/res)+1;
%read hourly wfde5 reanalysis data from netcdf file
wfde5_all=zeros(366*24,years_wfde5,length(xcol),length(parameters_wfde5),'single');
for varloop=1:length(parameters_wfde5)
parameter_wfde5=parameters_wfde5{varloop};
path=[wfde5dir,filesep,parameter_wfde5];
for yy=1:years_wfde5
year=yy+startyear_wfde5-1;
data_all=zeros(744,length(xcol),12,'single');
parfor mm=1:12
filelist=dir([path,filesep,parameter_wfde5,'*_',num2str(year),num2str(mm,'%.2i'),'*.nc']);
filename=[filelist.name];
cdir=[filelist.folder];
data_all(:,:,mm)=read_wfde5(cdir,filename,xcol,yrow,parameter_wfde5);
end%months
for mm=1:12
if(leap_year(year)==1)
dayspermonth = [31 29 31 30 31 30 31 31 30 31 30 31];
else
dayspermonth = [31 28 31 30 31 30 31 31 30 31 30 31];
end
b=sum(dayspermonth(1:mm))*24;
a=b-dayspermonth(mm)*24+1;
wfde5_all(a:b,yy,:,varloop)=data_all(1:dayspermonth(mm)*24,:,mm);
end
end%years
end%varloop
if(~exist([wfde5dir,filesep,'wfde5_tiles'],'dir'))
mkdir([wfde5dir,filesep,'wfde5_tiles'])
end
save([wfde5dir,filesep,'wfde5_tiles',filesep,'wfde5_all_',num2str(uppery),'_',num2str(lowery),'_',num2str(leftxx),'_',num2str(rightxx),'.mat'],'wfde5_all','-v7.3');
end
end
function wfde5=read_wfde5(cdir,filename,xcol,yrow,parameter_wfde5)
finfo = ncinfo([cdir,filesep,filename]);
varname=finfo.Variables(1,4).Name;
x_nc=finfo.Dimensions(1,find(strcmpi({finfo.Dimensions.Name},'lon')==1)).Length;
y_nc=finfo.Dimensions(1,find(strcmpi({finfo.Dimensions.Name},'lat')==1)).Length;
z_nc=finfo.Dimensions(1,find(strcmpi({finfo.Dimensions.Name},'time')==1)).Length;
data=zeros(y_nc,x_nc,z_nc);
data2d=zeros(length(xcol),z_nc);
wfde5=zeros(744,length(xcol));
disp(['read hourly WFDE5 file ',filename]);
data=ncread([cdir,filesep,filename],varname,[1 1 1],[720 360 z_nc]);
%data=rot90(data,1);
for xloop=1:length(xcol)
x=xcol(xloop);
y=yrow(xloop);
x_rot90=360-y+1;
y_rot90=x;
data2d(xloop,:)=squeeze(data(y_rot90,x_rot90,:));
end
wfde5(1:z_nc,:)=permute(data2d(:,1:z_nc),[2 1]);
%add snowfall to rainfall
if strcmpi(parameter_wfde5,'Rainf')==1
wfde5_snow=zeros(744,length(xcol));
filename=strrep(filename,'Rainf','Snowf');
cdir=strrep(cdir,'Rainf','Snowf');
finfo = ncinfo([cdir,filesep,filename]);
varname=finfo.Variables(1,4).Name;
disp(['read hourly WFDE5 file ',filename]);
data=ncread([cdir,filesep,filename],varname,[1 1 1],[720 360 z_nc]);
%data=rot90(data,1);
for xloop=1:length(xcol)
x=xcol(xloop);
y=yrow(xloop);
x_rot90=360-y+1;
y_rot90=x;
data2d(xloop,:)=data(y_rot90,x_rot90,:);
end
wfde5_snow(1:z_nc,:)=permute(data2d(:,1:z_nc),[2 1]);
wfde5(:,:)=wfde5(:,:)+wfde5_snow(:,:);
end
end