Let's say you have downloaded PERSIANN climate data files from the Centre of Hydrometerology and Remote Sensing.
Executing below mentioned codes on the downloaded .bin files will give you below mentioned output.
% this program is for training of CHRS graduate created by Prashant Kuranjekar
% last edit: 03/12/2007
% Usage: follow the hint in the screen and enjoy!
% Notice: You NEED yyyymmdd2yydoy.m to run me!
clear all;
%close all;
% read data
NDIR = '/mnt/data01/goesdata/globccs/nesdis/';
yyyy = input('Please input the year (format: yyyy):');
mm = input('Please input the month (format: mm):');
ddsta = input('Please input the start day (format: dd):');
hrsta = input('Please input the start hour (format: hr):');
mnsta = input('Please input the start minute (format: mn):');
rdcnt = input('Please input the number of records (48 for a day):');
lgtlw = input('Please input the lower boundary of longitude degree ("-" for West):');
lgtup = input('Please input the upper boundary of longitude degree ("-" for West):');
lttlw = input('Please input the lower boundary of latitude degree ("-" for North):');
lttup = input('Please input the upper boundary of latitude degree ("-" for North):');
[col_l,col_u]=mapping(lgtlw,lgtup,-180,180,1,9000);
[row_l,row_u]=mapping(lttlw,lttup,-60,60,1,3000);
col=col_l:1:col_u;
row=row_l:1:row_u;
mask=nan(3000,9000);
mask(row,col)=1;
if lgtlw<0
strlgtlw=sprintf('W%.3f',0-lgtlw);
else
strlgtlw=sprintf('E%.3f',lgtlw);
end
if lgtup<0
strlgtup=sprintf('W%.3f',0-lgtup);
else
strlgtup=sprintf('E%.3f',lgtup);
end
if lttlw<0
strlttlw=sprintf('N%.3f',0-lttlw);
else
strlgtlw=sprintf('S%.3f',lttlw);
end
if lttup<0
strlttup=sprintf('W%.3f',0-lttup);
else
strlttup=sprintf('E%.3f',lttup);
end
disp(['The selected boundary is ',strlgtlw,'~',strlgtup,', ',strlttlw,'~',strlttup]);
disp(['Relative zone in the matrix is row:',num2str(row_l),'~',num2str(row_u),', column:',num2str(col_l),'~',num2str(col_u)]);
%initialize
totp = zeros(3000,9000);
cnt = 0;
dd = ddsta;
hr = hrsta;
mn = mnsta;
scrsz = get(0,'ScreenSize');
figure1 = figure('Position',[1 scrsz(4)/2 (scrsz(3)-2)/2 scrsz(4)/2],'Name','[I] Changing Process of average data','NumberTitle','off');
figure2 = figure('Position',[scrsz(3)/2 scrsz(4)/2 (scrsz(3)-2)/2 scrsz(4)/2],'Name','[II] Changing Process of actural data','NumberTitle','off');
nperseries = zeros(3000:9000:rdcnt);
% data read and process
for rdnum = 1:1:rdcnt
yydoy = yyyymmdd2yydoy(yyyy,mm,dd);
shr = sprintf('%02d',hr);
smn = sprintf('%02d',mn);
gzfilename = [NDIR, 'rgccs', yydoy, shr, smn, '.bin.gz'];
filename = ['rgccs', yydoy, shr, smn, '.bin'];
disp(['Reading Data in ', filename,' ......']);
unix(['gunzip -c ', gzfilename, ' > ', filename]);
fid = fopen(filename, 'r', 'b');
npers = fread(fid, [9000 3000], 'int16', 'ieee-le');
fclose(fid);
npers = (npers/100)';
npers2 = npers;
npers3 = ones(3000,9000);
npers(find(npers<0)) = nan;
npers2(find(npers2<0)) = 0;
npers3(isnan(npers)) = 0;
nperseries(:,:,rdnum)=npers;
totp = totp+npers2;
cnt = cnt+npers3;
disp(['Read Data in ', filename, ' sucessfully!']);
unix(['rm ',filename]);
if mn == 15
mn = 45;
else
mn = 15;
hr = hr+1;
if hr == 24
hr = 00;
dd = dd+1;
end
end
% mean
avep=totp./cnt;
avep_selected=avep.*mask;
npers_selected=npers.*mask;
% plot daily data & mean
figure(figure1);
imagesc(avep_selected);
colorbar;
figure(figure2);
imagesc(npers_selected);
colorbar;
end
IsRepeat = input('Do you want to see the simulation of areal precipitaion changing process (Y/N):','s');
if IsRepeat == 'Y' || IsRepeat == 'y'
disp('The animation is shown in figure [II], press "ctrl+c" to exit!')
rptcnt = 0;
while 1
rdnum = rem(rptcnt,rdcnt)+1;
figure(figure2);
imagesc(nperseries(:,:,rdnum));
colorbar;
pause(0.5);
rptcnt = rptcnt+1;
end
else
close all;
disp('===================End of persiann simulation @Prashant Kuranjekar. 2007===================');
end
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
clear all;
RDIR = ['/mnt/data01/goesdata/globccs/nesdis/'];
%RDIR = ['/mnt/data01/goesdata/globccs/rapids/'];
%RDIR = ['/mnt/dans/goesdata/globccs/merge/out_post/'];
RES = 0.04;
%nesdis
hthr=20;
%dthr=8;
dthr=7;
%Rapids
%hthr=40;
%dthr=10;
%Ask for the input
% file = Erin;
% file3 = Ptot_hrly;
% file1 = Pcnt_dly;
% file2 = Pcnt;
%
% year = 2007;
% dstart= 221;
% dend = 232;
% hstart= 0;
% hend = 23;
%________________________________________________________________________
% lower_left_long = -110;
% lower_right_long = -85;
% lower_left_lat = 25;
% upper_left_lat = 38;
%Ask for the input
file = input('Enter the file\n','s');
file3 = input('Enter the filename for Ptot_hrly\n','s');
file1 = input('Enter the filename for Pcnt_dly\n','s');
file2 = input('Enter the filename for Pcnt\n','s');
year = input('Enter the year\n');
dstart= input('Enter the start day\n');
dend = input('Enter the end day\n');
hstart= input('Enter the start hour\n');
hend = input('Enter the end hour\n');
% %________________________________________________________________________
lower_left_long = input('Enter the lower left longitude\n');
lower_right_long = input('Enter the lower right longitude\n');
lower_left_lat = input('Enter the lower left latitude\n');
upper_left_lat = input('Enter the upper left latitude\n');
%________________________________________________________________________
if year<2000
year = num2str(year-1900);
yr=[year];
end
if year==2000
year = strcat(['0',num2str(year-2000)]);
yr=[year];
end
if 2000<year<2010
year = strcat(['0',num2str(year-2000)]);
yr=[year];
end
if year>=2010
year = num2str(year-2000);
yr=[year];
end
cellsize = 0.04;
ROWS = 3000;
COLS = 9000;
nodata_value = -99;
Xref = 0.02;
Yref = 59.8;
X = 0;
Y = 0;
Xref_adj = round(((Xref-X)/cellsize)+0.5)+4500;
Yref_adj = 1500-round((Yref-Y)/cellsize);
%______________________________________________________________________
%PART 1
if (lower_left_long>=-180 & lower_left_long<=180)
left_col= round(((lower_left_long-Xref)/cellsize)+0.5)+Xref_adj;
end
if (lower_right_long>=-180 & lower_right_long<=180)
Right_col= round(((lower_right_long-Xref)/cellsize)-0.5)+Xref_adj;
end
ncols=(Right_col-left_col)+1;
%PART 2
if (lower_left_lat>=-90 & lower_left_lat<=90)
lower_lat= Yref_adj+round((Yref-lower_left_lat)/cellsize);
end
if (upper_left_lat>=-90 & upper_left_lat<=90)
upper_lat= Yref_adj+1+round((Yref-upper_left_lat)/cellsize);
end
nrows = (lower_lat-upper_lat)+1;
Ptot=zeros(nrows,ncols);
Pcnt=zeros(nrows,ncols);
%___________________________________________________________________
for day = dstart:dend
if (day<10)
day = strcat(['0','0',num2str(day)]);
elseif (day<100)
day = strcat(['0',num2str(day)]);
else day = num2str(day);
day =[day];
end
%_____________________________________________________________________
Ptot_dly=zeros(nrows,ncols);
Pcnt_dly=zeros(nrows,ncols);
for hour=hstart:hend
if (hour<10)
hh = strcat(['0', num2str(hour)]);
else hh = num2str(hour);
end
%_____________________________________________________________________
Ptot_hrly=zeros(nrows,ncols);
for minutes=15:30:45
%for minutes=mstart:mend
if (minutes<10)
mm = strcat(['0', num2str(minutes)]);
else mm = num2str(minutes);
end
%_______________________________________________________________________
%_______________________________________________________________________
%_______________________________________________________________________
yrdayhhmm = strcat(yr,day,hh,mm);
%________________________________________________________________________
%str = [ '!gunzip -c ' RDIR 'rgccs' yrdayhhmm '.bin.gz > /tmp/rgccs' yrdayhhmm '.bin' ];
%eval(str);
str = [ '!gunzip -c ' RDIR 'rgccs' yrdayhhmm '.bin.gz > rgccs' yrdayhhmm '.bin' ];
eval(str);
%_________________________________________________________________________
fid=fopen(['rgccs' yrdayhhmm '.bin'],'r','b');
%if fid>=0
pers_mat1=fread(fid,[9000 3000],'short','ieee-le');
%pers_mat1=fread(fid,[9000 3000],'short','ieee-be');
pers_mat = pers_mat1/100;
fclose(fid);
str = [ '!rm rgccs' yrdayhhmm '.bin' ];
eval(str);
Ind = find(pers_mat<0);
pers_mat(Ind)=nan;
p = (pers_mat)';
a = [p(:,4501:9000) p(:,1:4500)];
a;
%end
%_____________________________________________________________________
m= [a(upper_lat:lower_lat,left_col:Right_col)];
m = 0.5*m;
%save m m -ascii;
%save(strcat([file,yr,day,hh,mm,'.asc']),'m','-ascii');
%_________________________________________________________________________
nanid=~isnan(m);
m(find(isnan(m)==1))=0;
%m=nanid.*m;
Ptot_hrly=Ptot_hrly+m;
Ptot_dly=Ptot_dly+m;
Pcnt_dly=Pcnt_dly+nanid;
%___________________________________________________________________________
end%minutes
save(strcat([file3,yr,day,hh,'.asc']),'Ptot_hrly','-ascii');
end%hrly
save(strcat([file1,yr,day,'.asc']),'Pcnt_dly','-ascii');
Ptot_dly=Ptot_dly./Pcnt_dly*(hend-hstart+1);
ind=find(Pcnt_dly<hthr);
Ptot_dly(ind)=NaN;
nanid=~isnan(Ptot_dly);
Ptot_dly(find(isnan(Ptot_dly)==1))=0;
save(strcat([file,yr,day,'.asc']),'Ptot_dly','-ascii');
%Ptot_dly=nanid.*Ptot_dly;
%Ptot=zeros(nrows,ncols);
%Pcnt=zeros(nrows,ncols);
Ptot=Ptot+Ptot_dly;
Pcnt=Pcnt+nanid;
end%day
Total_accumulated_days = dend-dstart+1;
Total_accumulated_days = num2str(Total_accumulated_days);
save(strcat([file2,yr,Total_accumulated_days,'.asc']),'Pcnt','-ascii');
Ptot=Ptot./Pcnt*(dend-dstart+1);
ind=find(Pcnt<dthr);
Ptot(ind)=NaN;
save(strcat([file,yr,Total_accumulated_days,'.asc']),'Ptot','-ascii');
imagesc(Ptot);
colorbar
% %read country boundary
% fid=fopen('country.bin','r','b');
% coun_bin=fread(fid,[2 478870/2],'float');
% country = (coun_bin)';
% Ind = find(country<-400);
% country(Ind)=nan;
% country;
% %plotting
% xx = (country(:,1)*25)+4500.5;
% yy = (country(:,2)*-25)+1500.5;
% xy = [xx yy];
% figure;
% imagesc(Ptot);
% hold on;
% plot(xy(:,1),xy(:,2),'w');
% xlabel('longitude');
% ylabel('Latitude');
% %title('Trend of snow water equivalent')
% vec1 = [1000 1000 1000 1000 1000 1000 1000 1000 1000];
% vec2 = [0 vec1];
% vec3 = cumsum(vec2);
% vectick=[' 0'; ' 40'; ' 80'; '120'; '160'; '200'; '240'; '280'; '320'; '360'];
% set(gca,'Xtick',vec3);
% set(gca,'XtickLabel',vectick);
% %_____________________________________________________________________
% vec4 = [500 500 500 500 500 500];
% vec5 = [0 vec4];
% vec6 = cumsum(vec5);
% vectick=['60'; '40'; '20'; ' 0'; '20'; '40'; '60'];
% set(gca,'Ytick',vec6);
% set(gca,'YtickLabel',vectick);
% %_________________________________________________________________________
% %k = strcat([file,yr,day,hh,'.jpeg']);
k = strcat([file,yr,Total_accumulated_days,'.bmp']);
%print('-djpeg',k);
print('-dbmp',k);
colorbar
Executing below mentioned codes on the downloaded .bin files will give you below mentioned output.
% this program is for training of CHRS graduate created by Prashant Kuranjekar
% last edit: 03/12/2007
% Usage: follow the hint in the screen and enjoy!
% Notice: You NEED yyyymmdd2yydoy.m to run me!
clear all;
%close all;
% read data
NDIR = '/mnt/data01/goesdata/globccs/nesdis/';
yyyy = input('Please input the year (format: yyyy):');
mm = input('Please input the month (format: mm):');
ddsta = input('Please input the start day (format: dd):');
hrsta = input('Please input the start hour (format: hr):');
mnsta = input('Please input the start minute (format: mn):');
rdcnt = input('Please input the number of records (48 for a day):');
lgtlw = input('Please input the lower boundary of longitude degree ("-" for West):');
lgtup = input('Please input the upper boundary of longitude degree ("-" for West):');
lttlw = input('Please input the lower boundary of latitude degree ("-" for North):');
lttup = input('Please input the upper boundary of latitude degree ("-" for North):');
[col_l,col_u]=mapping(lgtlw,lgtup,-180,180,1,9000);
[row_l,row_u]=mapping(lttlw,lttup,-60,60,1,3000);
col=col_l:1:col_u;
row=row_l:1:row_u;
mask=nan(3000,9000);
mask(row,col)=1;
if lgtlw<0
strlgtlw=sprintf('W%.3f',0-lgtlw);
else
strlgtlw=sprintf('E%.3f',lgtlw);
end
if lgtup<0
strlgtup=sprintf('W%.3f',0-lgtup);
else
strlgtup=sprintf('E%.3f',lgtup);
end
if lttlw<0
strlttlw=sprintf('N%.3f',0-lttlw);
else
strlgtlw=sprintf('S%.3f',lttlw);
end
if lttup<0
strlttup=sprintf('W%.3f',0-lttup);
else
strlttup=sprintf('E%.3f',lttup);
end
disp(['The selected boundary is ',strlgtlw,'~',strlgtup,', ',strlttlw,'~',strlttup]);
disp(['Relative zone in the matrix is row:',num2str(row_l),'~',num2str(row_u),', column:',num2str(col_l),'~',num2str(col_u)]);
%initialize
totp = zeros(3000,9000);
cnt = 0;
dd = ddsta;
hr = hrsta;
mn = mnsta;
scrsz = get(0,'ScreenSize');
figure1 = figure('Position',[1 scrsz(4)/2 (scrsz(3)-2)/2 scrsz(4)/2],'Name','[I] Changing Process of average data','NumberTitle','off');
figure2 = figure('Position',[scrsz(3)/2 scrsz(4)/2 (scrsz(3)-2)/2 scrsz(4)/2],'Name','[II] Changing Process of actural data','NumberTitle','off');
nperseries = zeros(3000:9000:rdcnt);
% data read and process
for rdnum = 1:1:rdcnt
yydoy = yyyymmdd2yydoy(yyyy,mm,dd);
shr = sprintf('%02d',hr);
smn = sprintf('%02d',mn);
gzfilename = [NDIR, 'rgccs', yydoy, shr, smn, '.bin.gz'];
filename = ['rgccs', yydoy, shr, smn, '.bin'];
disp(['Reading Data in ', filename,' ......']);
unix(['gunzip -c ', gzfilename, ' > ', filename]);
fid = fopen(filename, 'r', 'b');
npers = fread(fid, [9000 3000], 'int16', 'ieee-le');
fclose(fid);
npers = (npers/100)';
npers2 = npers;
npers3 = ones(3000,9000);
npers(find(npers<0)) = nan;
npers2(find(npers2<0)) = 0;
npers3(isnan(npers)) = 0;
nperseries(:,:,rdnum)=npers;
totp = totp+npers2;
cnt = cnt+npers3;
disp(['Read Data in ', filename, ' sucessfully!']);
unix(['rm ',filename]);
if mn == 15
mn = 45;
else
mn = 15;
hr = hr+1;
if hr == 24
hr = 00;
dd = dd+1;
end
end
% mean
avep=totp./cnt;
avep_selected=avep.*mask;
npers_selected=npers.*mask;
% plot daily data & mean
figure(figure1);
imagesc(avep_selected);
colorbar;
figure(figure2);
imagesc(npers_selected);
colorbar;
end
IsRepeat = input('Do you want to see the simulation of areal precipitaion changing process (Y/N):','s');
if IsRepeat == 'Y' || IsRepeat == 'y'
disp('The animation is shown in figure [II], press "ctrl+c" to exit!')
rptcnt = 0;
while 1
rdnum = rem(rptcnt,rdcnt)+1;
figure(figure2);
imagesc(nperseries(:,:,rdnum));
colorbar;
pause(0.5);
rptcnt = rptcnt+1;
end
else
close all;
disp('===================End of persiann simulation @Prashant Kuranjekar. 2007===================');
end
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
clear all;
RDIR = ['/mnt/data01/goesdata/globccs/nesdis/'];
%RDIR = ['/mnt/data01/goesdata/globccs/rapids/'];
%RDIR = ['/mnt/dans/goesdata/globccs/merge/out_post/'];
RES = 0.04;
%nesdis
hthr=20;
%dthr=8;
dthr=7;
%Rapids
%hthr=40;
%dthr=10;
%Ask for the input
% file = Erin;
% file3 = Ptot_hrly;
% file1 = Pcnt_dly;
% file2 = Pcnt;
%
% year = 2007;
% dstart= 221;
% dend = 232;
% hstart= 0;
% hend = 23;
%________________________________________________________________________
% lower_left_long = -110;
% lower_right_long = -85;
% lower_left_lat = 25;
% upper_left_lat = 38;
%Ask for the input
file = input('Enter the file\n','s');
file3 = input('Enter the filename for Ptot_hrly\n','s');
file1 = input('Enter the filename for Pcnt_dly\n','s');
file2 = input('Enter the filename for Pcnt\n','s');
year = input('Enter the year\n');
dstart= input('Enter the start day\n');
dend = input('Enter the end day\n');
hstart= input('Enter the start hour\n');
hend = input('Enter the end hour\n');
% %________________________________________________________________________
lower_left_long = input('Enter the lower left longitude\n');
lower_right_long = input('Enter the lower right longitude\n');
lower_left_lat = input('Enter the lower left latitude\n');
upper_left_lat = input('Enter the upper left latitude\n');
%________________________________________________________________________
if year<2000
year = num2str(year-1900);
yr=[year];
end
if year==2000
year = strcat(['0',num2str(year-2000)]);
yr=[year];
end
if 2000<year<2010
year = strcat(['0',num2str(year-2000)]);
yr=[year];
end
if year>=2010
year = num2str(year-2000);
yr=[year];
end
cellsize = 0.04;
ROWS = 3000;
COLS = 9000;
nodata_value = -99;
Xref = 0.02;
Yref = 59.8;
X = 0;
Y = 0;
Xref_adj = round(((Xref-X)/cellsize)+0.5)+4500;
Yref_adj = 1500-round((Yref-Y)/cellsize);
%______________________________________________________________________
%PART 1
if (lower_left_long>=-180 & lower_left_long<=180)
left_col= round(((lower_left_long-Xref)/cellsize)+0.5)+Xref_adj;
end
if (lower_right_long>=-180 & lower_right_long<=180)
Right_col= round(((lower_right_long-Xref)/cellsize)-0.5)+Xref_adj;
end
ncols=(Right_col-left_col)+1;
%PART 2
if (lower_left_lat>=-90 & lower_left_lat<=90)
lower_lat= Yref_adj+round((Yref-lower_left_lat)/cellsize);
end
if (upper_left_lat>=-90 & upper_left_lat<=90)
upper_lat= Yref_adj+1+round((Yref-upper_left_lat)/cellsize);
end
nrows = (lower_lat-upper_lat)+1;
Ptot=zeros(nrows,ncols);
Pcnt=zeros(nrows,ncols);
%___________________________________________________________________
for day = dstart:dend
if (day<10)
day = strcat(['0','0',num2str(day)]);
elseif (day<100)
day = strcat(['0',num2str(day)]);
else day = num2str(day);
day =[day];
end
%_____________________________________________________________________
Ptot_dly=zeros(nrows,ncols);
Pcnt_dly=zeros(nrows,ncols);
for hour=hstart:hend
if (hour<10)
hh = strcat(['0', num2str(hour)]);
else hh = num2str(hour);
end
%_____________________________________________________________________
Ptot_hrly=zeros(nrows,ncols);
for minutes=15:30:45
%for minutes=mstart:mend
if (minutes<10)
mm = strcat(['0', num2str(minutes)]);
else mm = num2str(minutes);
end
%_______________________________________________________________________
%_______________________________________________________________________
%_______________________________________________________________________
yrdayhhmm = strcat(yr,day,hh,mm);
%________________________________________________________________________
%str = [ '!gunzip -c ' RDIR 'rgccs' yrdayhhmm '.bin.gz > /tmp/rgccs' yrdayhhmm '.bin' ];
%eval(str);
str = [ '!gunzip -c ' RDIR 'rgccs' yrdayhhmm '.bin.gz > rgccs' yrdayhhmm '.bin' ];
eval(str);
%_________________________________________________________________________
fid=fopen(['rgccs' yrdayhhmm '.bin'],'r','b');
%if fid>=0
pers_mat1=fread(fid,[9000 3000],'short','ieee-le');
%pers_mat1=fread(fid,[9000 3000],'short','ieee-be');
pers_mat = pers_mat1/100;
fclose(fid);
str = [ '!rm rgccs' yrdayhhmm '.bin' ];
eval(str);
Ind = find(pers_mat<0);
pers_mat(Ind)=nan;
p = (pers_mat)';
a = [p(:,4501:9000) p(:,1:4500)];
a;
%end
%_____________________________________________________________________
m= [a(upper_lat:lower_lat,left_col:Right_col)];
m = 0.5*m;
%save m m -ascii;
%save(strcat([file,yr,day,hh,mm,'.asc']),'m','-ascii');
%_________________________________________________________________________
nanid=~isnan(m);
m(find(isnan(m)==1))=0;
%m=nanid.*m;
Ptot_hrly=Ptot_hrly+m;
Ptot_dly=Ptot_dly+m;
Pcnt_dly=Pcnt_dly+nanid;
%___________________________________________________________________________
end%minutes
save(strcat([file3,yr,day,hh,'.asc']),'Ptot_hrly','-ascii');
end%hrly
save(strcat([file1,yr,day,'.asc']),'Pcnt_dly','-ascii');
Ptot_dly=Ptot_dly./Pcnt_dly*(hend-hstart+1);
ind=find(Pcnt_dly<hthr);
Ptot_dly(ind)=NaN;
nanid=~isnan(Ptot_dly);
Ptot_dly(find(isnan(Ptot_dly)==1))=0;
save(strcat([file,yr,day,'.asc']),'Ptot_dly','-ascii');
%Ptot_dly=nanid.*Ptot_dly;
%Ptot=zeros(nrows,ncols);
%Pcnt=zeros(nrows,ncols);
Ptot=Ptot+Ptot_dly;
Pcnt=Pcnt+nanid;
end%day
Total_accumulated_days = dend-dstart+1;
Total_accumulated_days = num2str(Total_accumulated_days);
save(strcat([file2,yr,Total_accumulated_days,'.asc']),'Pcnt','-ascii');
Ptot=Ptot./Pcnt*(dend-dstart+1);
ind=find(Pcnt<dthr);
Ptot(ind)=NaN;
save(strcat([file,yr,Total_accumulated_days,'.asc']),'Ptot','-ascii');
imagesc(Ptot);
colorbar
% %read country boundary
% fid=fopen('country.bin','r','b');
% coun_bin=fread(fid,[2 478870/2],'float');
% country = (coun_bin)';
% Ind = find(country<-400);
% country(Ind)=nan;
% country;
% %plotting
% xx = (country(:,1)*25)+4500.5;
% yy = (country(:,2)*-25)+1500.5;
% xy = [xx yy];
% figure;
% imagesc(Ptot);
% hold on;
% plot(xy(:,1),xy(:,2),'w');
% xlabel('longitude');
% ylabel('Latitude');
% %title('Trend of snow water equivalent')
% vec1 = [1000 1000 1000 1000 1000 1000 1000 1000 1000];
% vec2 = [0 vec1];
% vec3 = cumsum(vec2);
% vectick=[' 0'; ' 40'; ' 80'; '120'; '160'; '200'; '240'; '280'; '320'; '360'];
% set(gca,'Xtick',vec3);
% set(gca,'XtickLabel',vectick);
% %_____________________________________________________________________
% vec4 = [500 500 500 500 500 500];
% vec5 = [0 vec4];
% vec6 = cumsum(vec5);
% vectick=['60'; '40'; '20'; ' 0'; '20'; '40'; '60'];
% set(gca,'Ytick',vec6);
% set(gca,'YtickLabel',vectick);
% %_________________________________________________________________________
% %k = strcat([file,yr,day,hh,'.jpeg']);
k = strcat([file,yr,Total_accumulated_days,'.bmp']);
%print('-djpeg',k);
print('-dbmp',k);
colorbar
No comments:
Post a Comment