I'm working with script code in python (in ArcGIS) for some processes from raster layers, my code takes the raster files, create arrays of data(NumPyArray), manipulate the data and then creates a rater output with the results. My problem is, this process is very slow. I'm working with basins and therefore a large amount of data in each array. I want to know if there is some other way to optimize my code or if this process is best to carry it out with the use of SQL tables. In this case, i dont know how to use the data in raster layers with SQL tables, please help
I'm working with the kinematic wave method to determine the runoff of a basin, I am using the explicit method for resolution (The Lax-Wendroff scheme)
This is my Code:
import arcpy , numpy , mathfrom arcpy import envfrom fractions import Fractionarcpy.CheckOutExtension("3D")env.workspace = "C:/data" #I give the program input values, a raster layer of Slope,#a raster layer of Aspect, a number of iterations determined by the user#and an address output to keep the output raster (I do not know what kind #of data is the address and I have to put it as "any value")Slope = arcpy.GetParameterAsText(0)Aspect = arcpy.GetParameterAsText(1)ii = arcpy.GetParameterAsText(2)Salida = arcpy.GetParameterAsText(3)#I create 9 arrays, each represents a variable of the method (which #varies for each cell according to variable)# Array of the Slope dataS_raster=arcpy.Raster(Slope)S_array=arcpy.RasterToNumPyArray(S_raster)(height, width)=S_array.shape# Array of the Aspect dataA_raster=arcpy.Raster(Aspect)A_array=arcpy.RasterToNumPyArray(A_raster)#Array created as a copy of "A_array" to modify their values laterAX_array=arcpy.RasterToNumPyArray(A_raster)#Array created as a copy of "A_array" to modify their values laterAY_array=arcpy.RasterToNumPyArray(A_raster)#Array created as a copy of "A_array" to modify their values laterN_array=arcpy.RasterToNumPyArray(A_raster)#n is the constant Manning roughness of the surface runoff in this case#will be constant, for more general future cases, this value should be#associated with the use of land in each celln=0.014#Array created as a copy of "A_array" to modify their values laterH_array=arcpy.RasterToNumPyArray(S_raster)#h0 is the height of water in the initial step, approximately zeroh0=0.00001#Array created as a copy of "A_array" to modify their values laterQ_array = arcpy.RasterToNumPyArray(S_raster) #Array created as a copy of "A_array" to modify their values laterDt_array = arcpy.RasterToNumPyArray(S_raster) CFL=0.5 # constant used in the method#Array created as a copy of "A_array" to modify their values laterR_array = arcpy.RasterToNumPyArray(S_raster) r=0.000001 # rainfall intensity used in the method (in mm/hr)#In this example, R and N are constant, but for more general cases, N is a# raster related to the use of land and R is a raster obtained by#interpolation of data obtained from weather stations. So they are actually #variables desc=arcpy.Describe(Slope)b=desc.meanCellHeight dx=desc.meanCellWidth g=9.81 #constant of gravityiterations = int(ii) Area=(height+(height*0.07))*(width+(width*0.07)) qmx=Fraction(qmx)expo=Fraction('3/5')qmxx= qmx ** expon=Fraction
expo=Fraction('3/5')nn= n ** expos=S_array.max()expo=Fraction('0.3')ss= (s/100) ** expoMax=(qmxx*nn)/ss #This is maximum flow possible by an approximation by rational method#This "for" modify the values in the arrangements created for usefor row in range(0,height): for col in range(0,width): #Insert data into Water Height Array H_array.itemset((row,col),h0) #Insert data into Roughness Array N_array.itemset((row,col),n) #Insert data into rain intensity Array R_array.itemset((row,col),r) #Insert data into X and Y Arrays of Aspect x=numpy.sin(math.radians(A_array.item(row,col))) AX_array.itemset((row,col),x) y=numpy.cos(math.radians(A_array.item(row,col))) AY_array.itemset((row,col),y) #It verified that the slope values are nonnegative s=S_array.item(row,col) if s= 0: if AY_array.item(row-1,col) < 0: q1 = AY_array.item(row-1,col)*Q_array.item(row-1,col) else: q1 = 0 else: q1=0 if col-1 >= 0 : if AX_array.item(row,col-1) > 0: q2 = AX_array.item(row,col-1)*Q_array.item(row,col-1) else: q2 = 0 else: q2 = 0 if col+1 < width : if AX_array.item(row,col+1) < 0: q3 = AX_array.item(row,col+1)*Q_array.item(row,col+1) else: q3 = 0 else: q3 = 0 if row+1 < height : if AY_array.item(row+1,col) > 0: q4 = AY_array.item(row+1,col)*Q_array.item(row+1,col) else: q4 = 0 else: q4 = 0 h=ha+minDt*(r-((q+(q1-q2+q3-q4))/dx)) if h>0: h=h else: h=0.000000000000000001 if h>Max: h=Max else : h=h H_array.itemset((row,col),h) for row in range(0,height): for col in range(0,width): #Calculation data of flow array n=N_array.item(row,col) h=H_array.item(row,col) s=S_array.item(row,col) h=Fraction(h) expo=Fraction('5/3') hh=h ** expo s=Fraction(s) expo=Fraction('1/2') ss=((s/100) ** expo) q=((b/n)*hh*ss)/b Q_array.itemset((row,col),q) #Calculation data of delta time array h=Fraction(h) expo=Fraction('1/2') hh=h ** expo g=Fraction(g) expo=Fraction('1/2') gg=g ** expo dt=CFL*(dx/(gg*hh)) Dt_array.itemset((row,col),dt) minDt = Dt_array.min()#Create the raster with final data flowC_Cau = arcpy.NumPyArrayToRaster(Q_array,None,b)C_Cau.save(Salida)
أكثر...
I'm working with the kinematic wave method to determine the runoff of a basin, I am using the explicit method for resolution (The Lax-Wendroff scheme)
This is my Code:
import arcpy , numpy , mathfrom arcpy import envfrom fractions import Fractionarcpy.CheckOutExtension("3D")env.workspace = "C:/data" #I give the program input values, a raster layer of Slope,#a raster layer of Aspect, a number of iterations determined by the user#and an address output to keep the output raster (I do not know what kind #of data is the address and I have to put it as "any value")Slope = arcpy.GetParameterAsText(0)Aspect = arcpy.GetParameterAsText(1)ii = arcpy.GetParameterAsText(2)Salida = arcpy.GetParameterAsText(3)#I create 9 arrays, each represents a variable of the method (which #varies for each cell according to variable)# Array of the Slope dataS_raster=arcpy.Raster(Slope)S_array=arcpy.RasterToNumPyArray(S_raster)(height, width)=S_array.shape# Array of the Aspect dataA_raster=arcpy.Raster(Aspect)A_array=arcpy.RasterToNumPyArray(A_raster)#Array created as a copy of "A_array" to modify their values laterAX_array=arcpy.RasterToNumPyArray(A_raster)#Array created as a copy of "A_array" to modify their values laterAY_array=arcpy.RasterToNumPyArray(A_raster)#Array created as a copy of "A_array" to modify their values laterN_array=arcpy.RasterToNumPyArray(A_raster)#n is the constant Manning roughness of the surface runoff in this case#will be constant, for more general future cases, this value should be#associated with the use of land in each celln=0.014#Array created as a copy of "A_array" to modify their values laterH_array=arcpy.RasterToNumPyArray(S_raster)#h0 is the height of water in the initial step, approximately zeroh0=0.00001#Array created as a copy of "A_array" to modify their values laterQ_array = arcpy.RasterToNumPyArray(S_raster) #Array created as a copy of "A_array" to modify their values laterDt_array = arcpy.RasterToNumPyArray(S_raster) CFL=0.5 # constant used in the method#Array created as a copy of "A_array" to modify their values laterR_array = arcpy.RasterToNumPyArray(S_raster) r=0.000001 # rainfall intensity used in the method (in mm/hr)#In this example, R and N are constant, but for more general cases, N is a# raster related to the use of land and R is a raster obtained by#interpolation of data obtained from weather stations. So they are actually #variables desc=arcpy.Describe(Slope)b=desc.meanCellHeight dx=desc.meanCellWidth g=9.81 #constant of gravityiterations = int(ii) Area=(height+(height*0.07))*(width+(width*0.07)) qmx=Fraction(qmx)expo=Fraction('3/5')qmxx= qmx ** expon=Fraction
expo=Fraction('3/5')nn= n ** expos=S_array.max()expo=Fraction('0.3')ss= (s/100) ** expoMax=(qmxx*nn)/ss #This is maximum flow possible by an approximation by rational method#This "for" modify the values in the arrangements created for usefor row in range(0,height): for col in range(0,width): #Insert data into Water Height Array H_array.itemset((row,col),h0) #Insert data into Roughness Array N_array.itemset((row,col),n) #Insert data into rain intensity Array R_array.itemset((row,col),r) #Insert data into X and Y Arrays of Aspect x=numpy.sin(math.radians(A_array.item(row,col))) AX_array.itemset((row,col),x) y=numpy.cos(math.radians(A_array.item(row,col))) AY_array.itemset((row,col),y) #It verified that the slope values are nonnegative s=S_array.item(row,col) if s= 0: if AY_array.item(row-1,col) < 0: q1 = AY_array.item(row-1,col)*Q_array.item(row-1,col) else: q1 = 0 else: q1=0 if col-1 >= 0 : if AX_array.item(row,col-1) > 0: q2 = AX_array.item(row,col-1)*Q_array.item(row,col-1) else: q2 = 0 else: q2 = 0 if col+1 < width : if AX_array.item(row,col+1) < 0: q3 = AX_array.item(row,col+1)*Q_array.item(row,col+1) else: q3 = 0 else: q3 = 0 if row+1 < height : if AY_array.item(row+1,col) > 0: q4 = AY_array.item(row+1,col)*Q_array.item(row+1,col) else: q4 = 0 else: q4 = 0 h=ha+minDt*(r-((q+(q1-q2+q3-q4))/dx)) if h>0: h=h else: h=0.000000000000000001 if h>Max: h=Max else : h=h H_array.itemset((row,col),h) for row in range(0,height): for col in range(0,width): #Calculation data of flow array n=N_array.item(row,col) h=H_array.item(row,col) s=S_array.item(row,col) h=Fraction(h) expo=Fraction('5/3') hh=h ** expo s=Fraction(s) expo=Fraction('1/2') ss=((s/100) ** expo) q=((b/n)*hh*ss)/b Q_array.itemset((row,col),q) #Calculation data of delta time array h=Fraction(h) expo=Fraction('1/2') hh=h ** expo g=Fraction(g) expo=Fraction('1/2') gg=g ** expo dt=CFL*(dx/(gg*hh)) Dt_array.itemset((row,col),dt) minDt = Dt_array.min()#Create the raster with final data flowC_Cau = arcpy.NumPyArrayToRaster(Q_array,None,b)C_Cau.save(Salida) أكثر...