add old arcgis basinmaker

basinmaker/arcgisguiwarpper/arcgis_basinmaker_old/Batch_differentAcc.py

@@ -0,0 +1,73 @@
+# Copyright 2018-2018 Ming Han - ming.han(at)uwaterloo.ca
+#
+# License
+# This file is part of Ming Han's personal code library.
+#
+# Ming Han's personal code library is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Ming Han's personal code library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+
+# You should have received a copy of the GNU Lesser General Public License
+# along with Ming Han's personal code library.  If not, see <http://www.gnu.org/licenses/>.
+#
+###########################################################################################
+#
+#This is a python script to use An automated ArcGIS toolbox for watershed delineation with lakes v1.tbx
+#This script show how to use this tool box to generate lake-river routing sturcture with non HydroSHEDS dataset.
+#
+###########################################################################################
+
+import numpy as np
+import pandas as pd
+import os
+
+####define required inputs
+#pythonfile = 'C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/Script for woring with HydroSHEDS DEM.py'
+#accfile = 'C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/accs.csv'
+#Workfolder = "C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/"
+
+
+pythonfile = sys.argv[1]
+accfile = sys.argv[2]
+Workfolder = sys.argv[3] + '/'
+
+os.chdir(Workfolder)####set working folder
+accs = pd.read_csv(accfile,sep=',')### read accs
+
+##### generate .bat file to run python script with arcgis py 
+#ofile2 = open(Workfolder+"runaccpys.bat", "w")
+#ofile2.write('cd' + '      ' + Workfolder +'\n')
+#ofile2.write(pythonexc + '      ' + 'acc_py_temp.py'+'\n')
+
+#ofile2.close()
+
+#####begin loop for each acc in the list
+
+### generate 
+for i in range(0,len(accs)):
+    ##### generate python script for ith acc value
+    os.chdir(Workfolder)
+    iacc = accs['Accs'][i]
+    ofile = open(Workfolder+"acc_py_temp.py", "w")
+    ifile = open(pythonfile, "r")
+    for line in ifile:
+        if 'accthres =' not in line:
+            ofile.write(line)
+        else:
+            ofile.write('accthres =' + "\""+str(iacc)+"\""+'\n')
+    ifile.close()
+    ofile.close()
+    arcpy.AddMessage("Current acc is    " + str(iacc) )
+    execfile('acc_py_temp.py') 
+    # run python script for ith acc value
+#    os.system('cd' + '      ' + Workfolder +'\n')
+#    os.system(pythonexc + '      ' + 'acc_py_temp.py'+'\n')    
+#    os.system("runaccpys.bat")
+
+#### extract subbasins for each point of interest of each acc value.

basinmaker/arcgisguiwarpper/arcgis_basinmaker_old/CalculateRivLoss.py

@@ -0,0 +1,117 @@
+
+import numpy as np
+import sys
+import os
+import csv
+import pandas as pd
+from simpledbf import Dbf5
+import matplotlib
+import matplotlib.pyplot as plt
+import copy as copy
+
+##### Readed inputs
+#POI_FIDS_file = "C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/poi.csv"
+#OutputF = "C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/Outputs1/"
+#RefFNam = 'RivLloss_ACC_1000'
+
+POI_FIDS_file = sys.argv[1]
+RefFNam = sys.argv[2]
+OutputF = sys.argv[3] + '/'
+
+
+pois =  pd.read_csv(POI_FIDS_file,sep=',')### read fid of pois
+pois_new =  pd.read_csv(POI_FIDS_file,sep=',')
+pois_new2 =  pd.read_csv(POI_FIDS_file,sep=',')
+
+pois['Total_DA'] = np.nan
+pois_new['Total_DA'] = np.nan
+pois_new2['Total_DA'] = np.nan
+
+#### check if need process with several delineation results
+if os.path.exists(OutputF + 'finalcat_info.shp'):
+    Allfolders = ['']
+else:
+    Allfolders = os.listdir(OutputF)
+
+for i in range(0,len(Allfolders)):
+    hyshdply = OutputF +Allfolders[i] + '/'+'finalcat_info.shp'
+    hyshdply.replace('//','/')
+
+    if not os.path.exists(OutputF + Allfolders[i]+'/'+ 'finalcat_info.shp'):
+        print (OutputF + Allfolders[i]+'/'+ 'finalcat_info.shp')
+        continue
+    ####read dbf file of the catchment of ith acc
+    tempinfo = Dbf5(hyshdply[:-3]+'dbf')#np.genfromtxt(hyinfocsv,delimiter=',')
+    hyshdinfo2 = tempinfo.to_dataframe()
+
+    pois[Allfolders[i]] = 0.00
+    pois_new[Allfolders[i]] = 0.00
+    pois_new2[Allfolders[i]] = 0.00
+
+#### loop for each poi
+    for j in range(0,len(pois)):
+        ipoi = pois['POI'][j]
+
+        tsubid = hyshdinfo2.loc[hyshdinfo2['IsObs'] == ipoi,]['SubId']
+
+        if len(tsubid) > 0:
+            ### define output folder for catchment of jth point of interest
+            ij_shp =  OutputF +Allfolders[i]+'/'+'POI_'+str(ipoi)+'/'+'finalcat_info.shp'
+            if not os.path.exists(ij_shp):
+                continue
+            ij_tempinfo = Dbf5(ij_shp[:-3]+'dbf')
+            ij_info = ij_tempinfo.to_dataframe()
+            if np.isnan(pois.loc[j,'Total_DA']):
+                pois.loc[j,'Total_DA'] = np.sum(ij_info['Area2'].values) ### calcuate DA for each point of interest
+                pois_new.loc[j,'Total_DA'] = np.sum(ij_info['Area2'].values)
+                pois_new2.loc[j,'Total_DA'] = np.sum(ij_info['Area2'].values)
+            for k in range(0,len(ij_info)):
+                catid = ij_info['SubId'].values[k]
+                upcatid = ij_info.loc[ij_info['DowSubId'] == catid,]['SubId'].values
+                if len(upcatid) > 0:
+                    pois.loc[j,Allfolders[i]] = pois.loc[j,Allfolders[i]] + max(ij_info['Rivlen'].values[k],0)*ij_info['Area2'].values[k]
+                    pois_new2.loc[j,Allfolders[i]] = pois_new2.loc[j,Allfolders[i]] + max(ij_info['Rivlen'].values[k],0)
+                else:
+                    pois.loc[j,Allfolders[i]] = pois.loc[j,Allfolders[i]] + 0.0 ##upstream cat
+                    pois_new2.loc[j,Allfolders[i]] = pois_new2.loc[j,Allfolders[i]] + 0.0 ##upstream cat
+############calculate new methods
+                ccatid = catid
+                downcatid = ij_info.loc[ij_info['SubId'] == ccatid,]['DowSubId'].values
+                while len(downcatid) > 0 and downcatid[0] > 0 and len(ij_info.loc[ij_info['SubId'] == downcatid[0],]['Rivlen'].values) > 0:
+#                    arcpy.AddMessage(max(ij_info.loc[ij_info['SubId'] == downcatid[0],]['Rivlen'].values,0))
+                    pois_new.loc[j,Allfolders[i]] = pois_new.loc[j,Allfolders[i]] + max(ij_info.loc[ij_info['SubId'] == downcatid[0],]['Rivlen'].values,0)*ij_info['Area2'].values[k]
+                    ccatid = downcatid[0]
+                    downcatid = ij_info.loc[ij_info['SubId'] == ccatid,]['DowSubId'].values
+############
+            pois.loc[j,Allfolders[i]] = pois.loc[j,Allfolders[i]]/pois.loc[j,'Total_DA']
+            pois_new.loc[j,Allfolders[i]] = pois_new.loc[j,Allfolders[i]]/pois_new.loc[j,'Total_DA']
+        else:
+            print('Point of interest is not inculded in subbasins:     ' + str(ipoi))
+
+pois.to_csv(OutputF+'rivweightlength_hongli.csv')
+pois_new.to_csv(OutputF+'rivweightlength_new.csv')
+pois_new2.to_csv(OutputF+'rivwlength_alone.csv')
+
+pois_new_relative = copy.deepcopy(pois_new)
+
+for i in range(0,len(Allfolders)):
+    if RefFNam == Allfolders[i] or Allfolders[i] not in pois.columns:
+        continue
+    else:
+        pois[Allfolders[i]] = pois[RefFNam].values - pois[Allfolders[i]].values
+        pois_new[Allfolders[i]] = (pois_new[RefFNam].values - pois_new[Allfolders[i]].values)
+        pois_new2[Allfolders[i]] = pois_new2[RefFNam].values - pois[Allfolders[i]].values
+        pois_new_relative[Allfolders[i]] = (pois_new_relative[RefFNam].values - pois_new_relative[Allfolders[i]].values)/pois_new_relative[RefFNam].values
+
+
+for i in range(0,len(Allfolders)):
+    if RefFNam == Allfolders[i] or Allfolders[i] not in pois.columns:
+        pois[Allfolders[i]] = 0.0
+        pois_new[Allfolders[i]] = 0.0
+        pois_new2[Allfolders[i]] = 0.0
+
+pois.to_csv(OutputF+'rivlenloss_hongli.csv')
+pois_new_relative.to_csv(OutputF+'rivlenloss_new_relative.csv')
+pois_new.to_csv(OutputF+'rivlenloss_new.csv')
+pois_new2.to_csv(OutputF+'rivlenloss_rivalone.csv')
+

basinmaker/arcgisguiwarpper/arcgis_basinmaker_old/Extract_Frompan_caproduct.py

@@ -0,0 +1,114 @@
+def dbftocsv(filename,outname):
+    if filename.endswith('.dbf'):
+#        print "Converting %s to csv" % filename
+        csv_fn = outname
+        with open(csv_fn,'wb') as csvfile:
+            in_db = dbf.Dbf(filename)
+            out_csv = csv.writer(csvfile)
+            names = []
+            for field in in_db.header.fields:
+                names.append(field.name)
+            out_csv.writerow(names)
+            for rec in in_db:
+                out_csv.writerow(rec.fieldData)
+            in_db.close()
+#            print "Done..."
+    else:
+        print "Filename does not end with .dbf"
+
+###################################################################3
+def Defcat(out,outletid):
+    otsheds = np.full((1,1),outletid)
+    Shedid = np.full((10000000,1),-99999999999999999)
+    psid = 0
+    rout = copy.copy(out)
+    while len(otsheds) > 0:
+        noutshd = np.full((10000000,1),-999999999999999)
+        poshdid = 0
+        for i in range(0,len(otsheds)):
+            Shedid[psid] = otsheds[i]
+            psid = psid + 1
+            irow = np.argwhere(rout[:,1]==otsheds[i]).astype(int)
+            for j in range(0,len(irow)):
+                noutshd[poshdid] = rout[irow[j],0]
+                poshdid = poshdid + 1
+        noutshd = np.unique(noutshd)
+        otsheds = noutshd[noutshd>=0]
+    Shedid = np.unique(Shedid)
+    Shedid = Shedid[Shedid>=0]
+    return Shedid
+
+def writeraster(w_filname,nraster,dataset):
+    orvh = open(w_filname,"w")
+    ncols = arcpy.GetRasterProperties_management(dataset, "COLUMNCOUNT")
+    nrows = arcpy.GetRasterProperties_management(dataset, "ROWCOUNT")
+    xllcorner = arcpy.GetRasterProperties_management(dataset, "LEFT")
+    yllcorner = arcpy.GetRasterProperties_management(dataset, "BOTTOM")
+    orvh.write("ncols      "+str(ncols) + "\n")
+    orvh.write("nrows      "+ str(nrows) + "\n")
+    orvh.write("xllcorner    "+str(xllcorner) + "\n")
+    orvh.write("yllcorner    "+str(yllcorner) + "\n")
+    orvh.write("cellsize     "+str(cellSize) + "\n")
+    orvh.write("NODATA_value  -9999" + "\n")
+    orvh.close()
+    f_handle = open(w_filname, 'a')
+    np.savetxt(f_handle,nraster,fmt='%i')
+    f_handle.close()
+
+
+import numpy as np
+from scipy.optimize import curve_fit
+import arcpy
+from arcpy import env
+from arcpy.sa import *
+import copy
+import sys
+import shutil
+import os
+import csv
+import csv
+from dbfpy import dbf
+import pandas as pd
+from shutil import copyfile
+from simpledbf import Dbf5
+arcpy.env.overwriteOutput = True
+arcpy.CheckOutExtension("Spatial")
+##### Readed inputs
+OutHyID = int(sys.argv[1])
+hyshdply = sys.argv[2]
+OutputFolder = sys.argv[3] + "/"
+OutHyID2 = int(sys.argv[4])
+tempinfo = Dbf5(hyshdply[:-3]+'dbf')#np.genfromtxt(hyinfocsv,delimiter=',')
+hyshdinfo2 = tempinfo.to_dataframe()
+hyshdinfo = hyshdinfo2[['SubId','DowSubId']].values
+
+SptailRef = arcpy.Describe(hyshdply).spatialReference
+
+arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(int(SptailRef.factoryCode)) ### WGS84
+
+if not os.path.exists(OutputFolder):
+    os.makedirs(OutputFolder)
+HydroBasins1 = Defcat(hyshdinfo,OutHyID)
+if OutHyID2 > 0:
+    HydroBasins2 = Defcat(hyshdinfo,OutHyID2)
+#    arcpy.AddMessage(HydroBasins1)
+#    arcpy.AddMessage(HydroBasins2)
+    for i in range(len(HydroBasins2)):
+        rows =np.argwhere(HydroBasins1 == HydroBasins2[i])
+#        arcpy.AddMessage(rows)
+        HydroBasins1 = np.delete(HydroBasins1, rows)
+#    arcpy.AddMessage(HydroBasins1)
+    HydroBasins = HydroBasins1
+else:
+    HydroBasins = HydroBasins1
+#arcpy.AddMessage(HydroBasins)
+out_feature_class = OutputFolder +"finalcat_info.shp"
+where_clause = '"SubId" IN'+ " ("
+for i in range(0,len(HydroBasins)):
+    if i == 0:
+        where_clause = where_clause + str(HydroBasins[i])
+    else:
+        where_clause = where_clause + "," + str(HydroBasins[i])
+where_clause = where_clause + ")"
+arcpy.Select_analysis(hyshdply, out_feature_class, where_clause)
+##################

basinmaker/arcgisguiwarpper/arcgis_basinmaker_old/ExtractbyPOI.py

@@ -0,0 +1,63 @@
+
+import numpy as np
+import arcpy
+from arcpy import env
+from arcpy.sa import *
+import sys
+import os
+import csv
+import pandas as pd
+from simpledbf import Dbf5
+arcpy.env.overwriteOutput = True
+arcpy.CheckOutExtension("Spatial")
+
+##### Readed inputs
+#POI_FIDS_file = "C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/poi.csv"
+#OutputF = "C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Samples/Examples/Outputs1/"
+#Pathoftoolbox = "C:/Users/dustm/Documents/ubuntu/share/OneDrive/OneDrive - University of Waterloo/Documents/RoutingTool/Code/Toolbox/An automated ArcGIS toolbox for watershed delineation with lakes.tbx"
+
+
+POI_FIDS_file = sys.argv[1]
+Pathoftoolbox = sys.argv[2]
+OutputF = sys.argv[3] + '/'
+
+arcpy.ImportToolbox(Pathoftoolbox)
+
+
+pois =  pd.read_csv(POI_FIDS_file,sep=',')### read fid of pois
+
+#### check if need process with several delineation results
+if os.path.exists(OutputF + 'finalcat_info.shp'):
+    Allfolders = ['']
+else:
+    Allfolders = os.listdir(OutputF)
+
+for i in range(0,len(Allfolders)):
+    hyshdply = OutputF +Allfolders[i] + '/'+'finalcat_info.shp'
+    hyshdply.replace('//','/')
+    arcpy.AddMessage(hyshdply[:-3]+'dbf')
+    if not os.path.exists(OutputF + Allfolders[i]+'/'+ 'finalcat_info.shp'):
+        print (OutputF + Allfolders[i]+'/'+ 'finalcat_info.shp')
+        continue
+    ####read dbf file of the catchment of ith acc
+    tempinfo = Dbf5(hyshdply[:-3]+'dbf')#np.genfromtxt(hyinfocsv,delimiter=',')
+    hyshdinfo2 = tempinfo.to_dataframe()
+
+#### loop for each poi
+    for j in range(0,len(pois)):
+        ipoi = pois['POI'][j]
+
+        tsubid = hyshdinfo2.loc[hyshdinfo2['IsObs'] == ipoi,]['SubId']
+
+        if len(tsubid) > 0:
+            ### define output folder for catchment of jth point of interest
+            OutputfoldercaseRes =  OutputF +Allfolders[i]+'/'+'POI_'+str(ipoi)+'/'
+            if not os.path.exists(OutputfoldercaseRes):
+                os.makedirs(OutputfoldercaseRes)
+#            arcpy.AddMessage(tsubid,ipoi,hyshdply)
+            arcpy.AutomatedLocalRoutingNetworkExtractionToolset(str(int(tsubid.values[0])),hyshdply,OutputfoldercaseRes,'-1')
+        else:
+            print('Point of interest is not inculded in subbasins:     ' + str(ipoi))
+
+
+