Thursday, April 23, 2015
Monday, April 20, 2015
Python - Converting multiple File geodatabase (.gdb) into multiple Personal geodatabase (.mdb)
Lets say that I have several file geodatabase in a folder and each of the file geodatabase has 1 feature dataset name "URS" and rest are feature classes. "URS" feature dataset has several feature class as well. Individually I can export the file geodatabase to XML workspace and then import that into the personnel geodatabase and all domains go over fine. But the task for me is to convert multiple file geodatabase into multiple personal geodatabase without including feature dataset,. If the geodatabase don't have "URS" feature dataset or any other feature classes then the python program will create empty personal geodatabase.
This is how we need to accomplish converting multiple file geodatabase (.gdb) into multiple personal geodatabase.
This is how we need to accomplish converting multiple file geodatabase (.gdb) into multiple personal geodatabase.
Here is the Pyhton code to accomplish this task.
If you don't want any Feature class from feature dataset, then this how u need to reflect the changes in the above mentioned code.
Thursday, April 16, 2015
Matlab - PERSIANN Climate Data Modeling
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
Monday, April 6, 2015
VB - Get the DTHETA value for each of the grid cell
This is a VB code inserted in a Manifold program. Here is the task to compute DTHETA value to assign each of the grid cell as per the soil condition. The weighted average DTHETA value for the each of the grid cell has been computed as well.
Function DTHETA
Sat = Record.Data("InitSat")
XKSAT = Record.DataText("AvgXKSAT")
If Sat = "Dry" Then
If XKSAT <= 0.15 Then
DTHETA = EXP(-0.2394+(0.3616*log(XKSAT)))
ElseIf XKSAT <= 0.25 Then
DTHETA = EXP(-1.4122-(0.2614*log(XKSAT)))
Else
DTHETA = 0.35
End If
ElseIf Sat = "Normal" Then
If XKSAT <= 0.02 Then
DTHETA = EXP(1.6094+log(XKSAT))
ElseIf XKSAT <= 0.04 Then
DTHETA = EXP(-0.0142+(0.5850*log(XKSAT)))
ElseIf XKSAT <= 0.1 Then
DTHETA = 0.15
ElseIf XKSAT <= 0.15 Then
DTHETA = EXP(1.0038+(1.2599*log(XKSAT)))
ElseIf XKSAT <= 0.4 Then
DTHETA = 0.25
Else
DTHETA = EXP(-1.2342+(0.1660*log(XKSAT)))
End If
Else
DTHETA = 0
End If
End Function
Function ADTHETA
D = Record.Data("DTHETA")
A = Record.Data("Area (I)")
ADTHETA = D*A
End Function
Function WDTHETA
W = Record.Data("ADTHETA")
A = Record.Data("Area (I)")
WDTHETA = W/A
End Function
Thursday, April 2, 2015
R - Transpose a Matrix
Let's say that we need to arrange a Hydrograph output in such a way that we can insert that results into HEC-1 input file.
The structure of hydrograph output file is 721 (Rows) by 1 (Columns)
Inorder to input the result into HEC-1, I need to arrange in n (Rows) by 10 (Columns)
Here is the R code..
workDir <- "C:/HEC-1 Installation/HEC1EXE/"
setwd(workDir)
inData <- read.csv("100yr24hr PLR 107cs.txt", header=FALSE, sep=",", as.is=TRUE)
#Mat1 <- matrix(inData$V1, 721, 10, byrow = T)
Mat1 <- matrix(inData$V1, ncol = 10, byrow = T)
write.table(Mat1, "100yr24hr PLR 107cs Transpose.txt", sep="\t", col.names=FALSE, row.names=FALSE)
The resultant file will be
Python (ArcPy) - Get rid of suffixes from multiple shapefiles in a folder
Let's say that I have several shapefiles in a folder and I would like to remove the suffixes from each shapefile and overwrite the previous files.
For Example:
Here is the Python (ArcPy) code..
import os, sys, arcpy
#InFolder = sys.argv[1] # make this a hard set path if you like
InFolder = "I:\\python\\GIS-data-Management\\Globalmapper" # make this a hard set path if you like
arcpy.env.workspace = InFolder
CapitalKeywordsToRemove = ["_AREAS","_LINES"]# MUST BE CAPITALS
DS_List = arcpy.ListFeatureClasses("*.shp","ALL") # Get a list of the feature classes (shapefiles)
for ThisDS in DS_List:
NewName = ThisDS # set the new name to the old name
# replace key words, searching is case sensitive but I'm trying not to change the case
# of the name so the new name is put together from the original name with searching in
# upper case, use as many key words as you like
for CapKeyWord in CapitalKeywordsToRemove:
if CapKeyWord in NewName.upper():
# remove the instance of CapKeyWord without changing the case
NewName = NewName[0:NewName.upper().find(CapKeyWord)] + NewName[NewName.upper().find(CapKeyWord) + len(CapKeyWord):]
if NewName != ThisDS:
if not arcpy.Exists(NewName):
arcpy.AddMessage("Renaming " + ThisDS + " to " + NewName)
arcpy.Rename_management(ThisDS , NewName)
else:
arcpy.AddWarning("Cannot rename, " + NewName + " already exists")
else:
arcpy.AddMessage("Retaining " + ThisDS)
Subscribe to:
Posts (Atom)