Thursday, April 23, 2015

Python - Batch process multiple MXD's (Map Documents) to PDF's

Lets say I have multiple MXD's into one folder and I would like to get a PDF of the layout view from each individual map documents.


Here is the Python code to accomplish this task. 


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.




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

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)