import time, sys, string, os, arcgisscripting, shutil from ModflowASCII import * from BoundaryCleaner import * def MODFLOWtxt(cellsize,DEM_cellsize,Recharge_cellsize,dem_minus,working_folder,outDirec,names): start = time.clock() gp = arcgisscripting.create(9.3) #Turn creation of data on or off create_AquiferPolygon = True create_DEM = True create_SPYD = True create_K = True create_Bedrock = True create_Recharge = True create_InitialHead = True create_RiverHead = True create_RiverBottom = True create_Conductance = True create_ModelBoundary = True create_SlopingBaseBoundary = True create_Wells = True SPYD = working_folder + names[0] Hyd_K = working_folder + names[1] Bedrock = working_folder + names[2] Recharge = working_folder + names[3] PredevWaterLv = working_folder + names[4] HighPlainsRivers = working_folder + names[5] DEM = working_folder + names[6] AquiferCoverage = working_folder + names[7] ToolboxFolder = names[8] Spatial_Ref_Folder = names[9] #Enter Paths """ KS Ogallala Model ToolboxFolder = "C:/Program Files/arcgis/Desktop10.0/ArcToolbox/Toolboxes/" Spatial_Ref_Folder = "C:/Program Files/arcgis/Desktop10.0/Coordinate Systems/Projected Coordinate Systems/UTM/NAD 1983/" SPYD = working_folder + "SPYD.shp" Hyd_K = working_folder + "Hyd_K.shp" Bedrock = working_folder + "bedrock_contours.shp" Recharge = working_folder + "rech_in_yr_ss" PredevWaterLv = working_folder + "pred_wlv.shp" HighPlainsRivers = working_folder + "High_Plains_Rivers.shp" DEM = working_folder + "ksogal_dem100" AquiferCoverage = working_folder + "OgallalaExtents_Clip.shp" """ """ HighPlains Model ToolboxFolder = "C:/Program Files/arcgis/Desktop10.0/ArcToolbox/Toolboxes/" Spatial_Ref_Folder = "C:/Program Files/arcgis/Desktop10.0/Coordinate Systems/Projected Coordinate Systems/UTM/NAD 1983/" SPYD = working_folder + "SPYD.shp" Hyd_K = working_folder + "Hyd_K.shp" Bedrock = working_folder + "bedrock_contours.shp" Recharge = working_folder + "rech_in_yr_ss" PredevWaterLv = working_folder + "pred_wlv.shp" HighPlainsRivers = working_folder + "HighPlainsMajorRiversNHD.shp" DEM = working_folder + "dem_250m" AquiferCoverage = working_folder + "hp_extents.shp" """ """ 3pt HighPlains Models ToolboxFolder = "C:/Program Files/arcgis/Desktop10.0/ArcToolbox/Toolboxes/" Spatial_Ref_Folder = "C:/Program Files/arcgis/Desktop10.0/Coordinate Systems/Projected Coordinate Systems/UTM/NAD 1983/" SPYD = working_folder + "SPYD.shp" Hyd_K = working_folder + "Hyd_K.shp" Bedrock = working_folder + "bedrock_contours.shp" Recharge = working_folder + "rech_in_yr_ss" PredevWaterLv = working_folder + "pred_wlv.shp" HighPlainsRivers = working_folder + "HighPlainsMajorRiversNHD.shp" DEM = working_folder + "dem_250m" AquiferCoverage = working_folder + "north_hp_extents.shp" #AquiferCoverage = working_folder + "central_hp_extents.shp" #AquiferCoverage = working_folder + "south_hp_extents.shp" """ #Create folder to hold new model information dirname = outDirec iteration = 0 try: os.makedirs(dirname) except OSError: while os.path.exists(dirname): iteration += 1 dirname = outDirec + "_" + str(iteration) try: os.makedirs(dirname) except OSError: raise os.makedirs(dirname + "/GIS") os.makedirs(dirname + "/MODFLOW_txt") os.makedirs(dirname + "/GIS/Temp") os.makedirs(dirname + "/txt") print "Working Directory:",working_folder print "Output Directory:", dirname print "cell size =",cellsize #Create Aquifer Polygon #Extents desc = gp.Describe(AquiferCoverage) extent = desc.Extent top = float(extent.ymax) bottom = float(extent.ymin) left = float(extent.xmin) right = float(extent.xmax) if top < 0: top = int(top) - 1 else: top = int(top) + 1 if bottom < 0 : bottom = int(bottom) - 1 else: bottom = int(bottom) + 1 if left < 0: left = int(left) - 1 else: left = int(left) + 1 if right < 0: right = int(right) - 1 else: right = int(right) + 1 #Determine required size of polygon width = float(abs(left - right)) new_width = width while new_width/cellsize - int(new_width/cellsize) <> 0: new_width += 1.0 #Adjust left and right extents change = (new_width - width)/2 if left < 0: left = abs(left) + int(change) left = -left else: left = left + int(change) if right < 0: if change - int(change) <> 0: right = abs(right) + int(change) + 1 right = -right else: right = abs(right) + int(change) right = -right else: if change - int(change) <> 0: right = right + int(change) + 1 else: right = right + int(change) length = float(abs(top - bottom)) new_length = length while new_length/cellsize - int(new_length/cellsize) <> 0: new_length += 1.0 #Adjust top and bottom extents change = (new_length - length)/2 if top < 0: top = abs(top) + int(change) top = -top else: top = top + int(change) if bottom < 0: if change - int(change) <> 0: bottom = abs(bottom) - int(change) - 1 bottom = -bottom else: bottom = abs(bottom) - int(change) bottom = -bottom else: if change - int(change) <> 0: bottom = bottom - int(change) - 1 else: bottom = bottom - int(change) left = left - cellsize*10 right = right + cellsize*10 top = top + cellsize*10 bottom = bottom - cellsize*10 outFolder = dirname + "/GIS" gp.workspace = outFolder textFolders = dirname #Create AquiferPolygon shapefile if create_AquiferPolygon or create_ModelBoundary or create_Conductance: print 'Creating Aquifer Polygon files' try: outFile = "AquiferPolygon.shp" StartandEndXcoord = left StartandEndYcoord = top PointOneXcoord = right PointOneYcoord = top PointTwoXcoord = right PointTwoYcoord = bottom PointThreeXcoord = left PointThreeYcoord = bottom sr = gp.CreateSpatialReference(Spatial_Ref_Folder + "NAD 1983 UTM Zone 14N.prj", "#", "#", "#", "#", "#") gp.CreateFeatureclass_management(outFolder, outFile, "POLYGON") cur = gp.InsertCursor(outFile) row = cur.NewRow() PolygonArray = gp.CreateObject("Array") pnt = gp.CreateObject("Point") pnt.x = StartandEndXcoord pnt.y = StartandEndYcoord PolygonArray.add(pnt) pnt.x = PointOneXcoord pnt.y = PointOneYcoord PolygonArray.add(pnt) pnt.x = PointTwoXcoord pnt.y = PointTwoYcoord PolygonArray.add(pnt) pnt.x = PointThreeXcoord pnt.y = PointThreeYcoord PolygonArray.add(pnt) pnt.x = StartandEndXcoord pnt.y = StartandEndYcoord PolygonArray.add(pnt) row.shape = PolygonArray cur.InsertRow(row) del row, cur except: raise AquiferPolygon = outFolder + "/AquiferPolygon.shp" #Create DEM if DEM_cellsize <> cellsize: print 'Resampling DEM to match input cell size' try: if DEM_cellsize > cellsize: print "Dem cellsize > model cellsize" else: # Local variables... DEM_resample = outFolder + "/dem" # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Data Management Tools.tbx") # Process: Resample... gp.Resample_management(DEM, DEM_resample, str(cellsize), "NEAREST") DEM = DEM_resample except: raise else: print 'Copying DEM to output folder' try: # Local variables... DEM_resample = outFolder + "/dem" gp.CopyRaster_management(DEM, DEM_resample, "", "", "", "NONE", "NONE", "") except: raise gp.GetMessage(2) #Create Specific Yield files if create_SPYD: print 'Creating Specific Yield files' try: # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Conversion Tools.tbx") #Local Variables Raster_SPYD = outFolder + "/SPYD" gp.Extent = str(left) + " " + str(bottom) + " " + str(right) + " " + str(top) SPYD_ASCII_file = textFolders + "/txt/spyd.txt" # Process: Feature to Raster... gp.FeatureToRaster_conversion(SPYD, "AVG_SPYD", Raster_SPYD, str(cellsize)) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Raster_SPYD, SPYD_ASCII_file) # Edit text file MODFLOW_ascii(SPYD_ASCII_file, textFolders + "/MODFLOW_txt/SPYD_MODFLOW.txt") del Raster_SPYD,SPYD_ASCII_file except: raise #Create Hydraulic Conductivity files if create_K: print 'Creating Hydraulic Conductivity files' try: # Allow overwrite gp.OverwriteOutput = True gp.workspace = outFolder + "/Temp/" # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Conversion Tools.tbx") gp.AddToolbox(ToolboxFolder + "Spatial Analyst Tools.tbx") #Local Variables gp.Extent = str(left) + " " + str(bottom) + " " + str(right) + " " + str(top) Raster_K = outFolder + "/Temp/hyd_K" K1_md = outFolder + "\\Temp\\k1_md" K2_md = outFolder + "\\Temp\\k2_md" K3_md = outFolder + "\\k_md" K_ASCII_file = textFolders + "\\txt\\k.txt" # Process: Feature to Raster... gp.FeatureToRaster_conversion(Hyd_K, "AVG_K", Raster_K, str(cellsize)) #print 'times' # Process: Times... feet_to_meters = .3048 gp.Times_sa(Raster_K, feet_to_meters, K1_md) # Process: Single Output Map Algebra... gp.SingleOutputMapAlgebra_sa("con(isnull(k1_md),17,k1_md)", K2_md, "") # Process: Single Output Map Algebra... gp.SingleOutputMapAlgebra_sa("con(k2_md == 0,1,k2_md)", K3_md, "") # Process: Raster to ASCII... gp.RasterToASCII_conversion(K3_md, K_ASCII_file) # Edit text file MODFLOW_ascii(K_ASCII_file, textFolders + "\\MODFLOW_txt\\K_MODFLOW.txt") gp.workspace = outFolder del Raster_K, K1_md, K2_md, K3_md, K_ASCII_file except: raise #Create Bedrock files if create_Bedrock or create_ModelBoundary or create_RiverHead or create_Conductance: print 'Creating Bedrock files' try: # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Spatial Analyst Tools.tbx") # Local variables... bedrock_raster = outFolder + "\\Temp\\bedrock_elev" Output_stream_polyline_features = "" Output_remaining_sink_point_features = "" Output_diagnostic_file = "" Output_parameter_file = "" bedrock_raster_m = outFolder + "\\bedrock_m" Bedrock_ASCII_file = textFolders + "\\txt\\bedrock.txt" #Create bedrock raster from topo lines gp.TopoToRaster_sa(Bedrock + " ELEV Contour", bedrock_raster, str(cellsize), AquiferCoverage, "20", "", "", "ENFORCE", "CONTOUR", "40", "", "1", "0", "", "", Output_stream_polyline_features, Output_remaining_sink_point_features, Output_diagnostic_file, Output_parameter_file) # Process: Times... feet_to_meters = .3048 gp.Times_sa(bedrock_raster, feet_to_meters, bedrock_raster_m) # Process: Raster to ASCII... gp.RasterToASCII_conversion(bedrock_raster_m, Bedrock_ASCII_file) # Edit text file MODFLOW_ascii(Bedrock_ASCII_file, textFolders + "\\MODFLOW_txt\\Bedrock_MODFLOW.txt") del bedrock_raster, Output_stream_polyline_features, Output_remaining_sink_point_features, Output_diagnostic_file del Output_parameter_file, Bedrock_ASCII_file except: raise #Create Predevelopment Water Level files if create_InitialHead or create_ModelBoundary or create_RiverHead or create_Conductance: print 'Creating Initial Head files' try: # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Spatial Analyst Tools.tbx") # Local variables... thickness = outFolder + "\\Temp\\thickness" thickness_reclass = outFolder + "\\thickness_rec" pred_wlv_raster = outFolder + "\\Temp\\pred_wlv" Output_stream_polyline_features = "" Output_remaining_sink_point_features = "" Output_diagnostic_file = "" Output_parameter_file = "" pred_wlv_raster_m = outFolder + "\\pred_wlv_m" InitialHead_raster = outFolder + "\\initialhead_m" InitialHead_ASCII_file = textFolders + "\\txt\\initialhead.txt" PredWlv_ASCII_file = textFolders + "\\txt\\predwlv.txt" #Create Predeveolopment Water Level Raster for future comparison # Process: Topo to Raster... gp.TopoToRaster_sa(PredevWaterLv + " ELEV Contour", pred_wlv_raster, str(cellsize), AquiferCoverage, "20", "", "", "ENFORCE", "CONTOUR", "60", "", "0.75", "0", "", "", Output_stream_polyline_features, Output_remaining_sink_point_features, Output_diagnostic_file, Output_parameter_file) # Process: Times... feet_to_meters = .3048 gp.Times_sa(pred_wlv_raster, feet_to_meters, pred_wlv_raster_m) # Process: Minus... #gp.Minus_sa(InitialHead_raster, bedrock_raster_m, thickness) #Create Initial Head #gp.Minus_sa(DEM, str(dem_minus), InitialHead_raster) gp.Minus_sa(pred_wlv_raster_m, bedrock_raster_m, thickness) #for use with optimization scripts # Process: Reclassify... gp.Reclassify_sa(thickness, "Value", "-10000 5 0;5 10000 1", thickness_reclass, "DATA") # Process: Raster to ASCII... gp.RasterToASCII_conversion(thickness_reclass, InitialHead_ASCII_file) # Process: Raster to ASCII... gp.RasterToASCII_conversion(pred_wlv_raster_m, PredWlv_ASCII_file) # Edit text file MODFLOW_ascii(InitialHead_ASCII_file, textFolders + "\\MODFLOW_txt\\InitialHead_MODFLOW.txt") MODFLOW_ascii(PredWlv_ASCII_file, textFolders + "\\MODFLOW_txt\\PredWlv_MODFLOW.txt") del pred_wlv_raster, Output_stream_polyline_features, Output_remaining_sink_point_features, Output_diagnostic_file del Output_parameter_file, pred_wlv_raster_m, InitialHead_ASCII_file, PredWlv_ASCII_file except: raise #Create Recharge files if create_Recharge: if Recharge_cellsize <> cellsize: print 'Resampling Recharge Raster' try: if Recharge_cellsize > cellsize: print "Recharge cellsize > model cellsize" else: # Local variables... Recharge_clip = outFolder + "\\Temp\\recharge_clip" Recharge_resample = outFolder + "\\Temp\\recharge_res" Recharge_md = outFolder + "\\recharge_md" Recharge_ASCII_file = textFolders + "\\txt\\recharge.txt" # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Data Management Tools.tbx") gp.AddToolbox(ToolboxFolder + "Spatial Analyst Tools.tbx") # Process: Clip (2)... gp.Clip_management(Recharge, str(left) + " " + str(bottom) + " " + str(right) + " " + str(top), Recharge_clip, AquiferPolygon, "", "ClippingGeometry") # Process: Resample... gp.Resample_management(Recharge_clip, Recharge_resample, str(cellsize), "NEAREST") # Process: Times... inchesPERyear_to_metersPERday = 0.0254/365 gp.Times_sa(Recharge_resample, inchesPERyear_to_metersPERday, Recharge_md) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Recharge_md, Recharge_ASCII_file) # Edit text file MODFLOW_ascii(Recharge_ASCII_file, textFolders + "\\MODFLOW_txt\\Recharge_MODFLOW.txt") del Recharge_clip,Recharge_resample,Recharge_md,Recharge_ASCII_file except: raise else: print 'Creating Recharge files' try: Recharge_md = outFolder + "\\recharge_md" Recharge_ASCII_file = textFolders + "\\txt\\recharge.txt" # Process: Times... inchesPERyear_to_metersPERday = 0.0254/365 gp.Times_sa(Recharge, inchesPERyear_to_metersPERday, Recharge_md) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Recharge_md, Recharge_ASCII_file) # Edit text file MODFLOW_ascii(Recharge_ASCII_file, textFolders + "\\MODFLOW_txt\\Recharge_MODFLOW.txt") del Recharge_md,Recharge_ASCII_file except: raise gp.GetMessage(2) #Create Boundary Files if create_ModelBoundary or create_RiverHead or create_Conductance or create_SlopingBaseBoundary: print 'Creating Model Boundary files' try: # Allow overwrite gp.OverwriteOutput = True # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Spatial Analyst Tools.tbx") gp.AddToolbox(ToolboxFolder + "Conversion Tools.tbx") gp.AddToolbox(ToolboxFolder + "Data Management Tools.tbx") # Local variables... gp.Extent = str(left) + " " + str(bottom) + " " + str(right) + " " + str(top) AquiferCoverage_raster = outFolder + "\\Temp\\aquifer_ras" Extents_raster = outFolder + "\\Temp\\model_extents" Extents_raster_reclassify = outFolder + "\\Temp\\mod_ext_rec" Extents_Polygon = outFolder + "\\Temp\\extents_poly.shp" Single_Model_Polygon = outFolder + "\\Temp\\single_boundary_polygon.shp" Single_Model_Polygon_Agg = outFolder + "\\Temp\\single_boundary_polygon_agg.shp" Single_Model_Polygon_Agg_Buffer = outFolder + "\\Temp\\single_boundary_polygon_agg_buffer.shp" Spec_Head_Bound_Poly = outFolder + "\\Temp\\spec_head_bound.shp" Single_Model_Raster = outFolder + "\\Temp\\bound_raster1" Single_Model_Raster_Agg = outFolder + "\\Temp\\bound_raster2" Single_Model_Raster_Agg_Buffer = outFolder + "\\Temp\\bound_raster3" Single_Model_Raster_Reclass = outFolder + "\\Temp\\bound_ras1_rc" Single_Model_Raster_Agg_Reclass = outFolder + "\\Temp\\bound_ras2_rc" Single_Model_Raster_Agg_Buffer_Reclass = outFolder + "\\Temp\\bound_ras3_rc" Bound_Raster4 = outFolder + "\\Temp\\bound_raster4" Bound_Raster5 = outFolder + "\\Temp\\bound_raster5" Bound_Raster6 = outFolder + "\\Temp\\bound_raster6" Boundary_Raster = outFolder + "\\boundary" Boundary_Centroids = outFolder + "\\Temp\\boundary_centroids.shp" Bound_X = outFolder + "\\Temp\\bound_x" Bound_Y = outFolder + "\\Temp\\bound_y" X_Coords = textFolders + "\\txt\\X_coords.txt" Y_Coords = textFolders + "\\txt\\Y_coords.txt" Boundary_ASCII_file = textFolders + "\\txt\\boundary.txt" # Process: Feature to Raster... fields = gp.ListFields(AquiferCoverage) names = [] for field in fields: names.append(field.Name) if "Raster_val" in names: gp.FeatureToRaster_conversion(AquiferCoverage, "Raster_val", AquiferCoverage_raster, str(cellsize)) else: gp.AddField(AquiferCoverage, "Raster_val", "short") gp.CalculateField_management (AquiferCoverage, "Raster_val", "1", "PYTHON_9.3") gp.FeatureToRaster_conversion(AquiferCoverage, "Raster_val", AquiferCoverage_raster, str(cellsize)) # Process: Plus... gp.Plus_sa(AquiferCoverage_raster, thickness_reclass, Extents_raster) # Process: Reclassify... gp.Reclassify_sa(Extents_raster, "Value", "1 1;2 2;NODATA 1", Extents_raster_reclassify, "DATA") # Process: Extents_raster_reclassify to Polygon... gp.RasterToPolygon_conversion(Extents_raster_reclassify, Extents_Polygon, "NO_SIMPLIFY", "VALUE") # Add Area Field gp.AddField_management (Extents_Polygon, "Area", "FLOAT", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "") # Calculate Area expression = "float(!SHAPE.AREA@SQUAREMILES!)" gp.CalculateField_management (Extents_Polygon, "Area", expression, "PYTHON_9.3") # Find, select, and save largest polygon within boundary to new shapefile: searchRows = gp.searchcursor(Extents_Polygon) searchRow = searchRows.next() Area = [] while searchRow: if searchRow.GRIDCODE == 2: Area.append(searchRow.Area) searchRow = searchRows.next() maxArea = max(Area) area_string = "\"Area\" = " + str(maxArea) area_string = "Area > " + str(int(maxArea) - 1) + " AND Area < " + str(int(maxArea) + 1) gp.MakeFeatureLayer(Extents_Polygon,"bound_lyr") gp.SelectLayerByAttribute_management("bound_lyr", "NEW_SELECTION", area_string) gp.CopyFeatures("bound_lyr", Single_Model_Polygon) # Process: Single_Model_Polygon Polygons... gp.AggregatePolygons_management(Single_Model_Polygon, Single_Model_Polygon_Agg, str(int(cellsize)) + " Unknown", "0 Unknown", "1E+20 SquareMiles", "NON_ORTHOGONAL") # Process: Buffer Single_Model_Polygon... gp.Buffer_analysis(Single_Model_Polygon_Agg, Single_Model_Polygon_Agg_Buffer, str(-int(cellsize)) + " Meters", "FULL", "FLAT", "NONE", "") # Add Raster_val Field to Single_Model_Polygon gp.AddField_management (Single_Model_Polygon, "Raster_val", "SHORT", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "") # Add Raster_val Field to Single_Model_Polygon_Agg gp.AddField_management (Single_Model_Polygon_Agg, "Raster_val", "SHORT", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "") # Add Raster_val Field to Single_Model_Polygon_Agg_Buffer gp.AddField_management (Single_Model_Polygon_Agg_Buffer, "Raster_val", "SHORT", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "") # Process: Calculate Field for Single_Model_Polygon... gp.CalculateField_management(Single_Model_Polygon, "Raster_val", "1", "PYTHON_9.3", "") # Process: Calculate Field for Single_Model_Polygon_Agg... gp.CalculateField_management(Single_Model_Polygon_Agg, "Raster_val", "1", "PYTHON_9.3", "") # Process: Calculate Field for Single_Model_Polygon_Agg_Buffer... gp.CalculateField_management(Single_Model_Polygon_Agg_Buffer, "Raster_val", "1", "PYTHON_9.3", "") # Process: PolygonToRaster... gp.PolygonToRaster_conversion(Single_Model_Polygon, "Raster_val", Single_Model_Raster, "CELL_CENTER", "NONE", str(cellsize)) # Process: Reclassify... gp.Reclassify_sa(Single_Model_Raster, "VALUE", "1 1;NODATA 0", Single_Model_Raster_Reclass, "DATA") # Process: PolygonToRaster... gp.PolygonToRaster_conversion(Single_Model_Polygon_Agg, "Raster_val", Single_Model_Raster_Agg, "CELL_CENTER", "NONE", str(cellsize)) # Process: Reclassify... gp.Reclassify_sa(Single_Model_Raster_Agg, "VALUE", "1 1;NODATA 0", Single_Model_Raster_Agg_Reclass, "DATA") # Process: PolygonToRaster... gp.PolygonToRaster_conversion(Single_Model_Polygon_Agg_Buffer, "Raster_val", Single_Model_Raster_Agg_Buffer, "CELL_CENTER", "NONE", str(cellsize)) # Process: Reclassify... gp.Reclassify_sa(Single_Model_Raster_Agg_Buffer, "VALUE", "1 1;NODATA 0", Single_Model_Raster_Agg_Buffer_Reclass, "DATA") # Add Single_Model_Polygon_Agg to Single_Model_Polygon_Agg_Buffer gp.Plus_sa(Single_Model_Raster_Agg_Reclass, Single_Model_Raster_Agg_Buffer_Reclass, Bound_Raster4) # Process: Reclassify... gp.Reclassify_sa(Bound_Raster4, "VALUE", "0 0;1 10;2 2", Bound_Raster5, "DATA") # Add Single_Model_Polygon to Bound_Raster5 gp.Plus_sa(Single_Model_Raster_Reclass, Bound_Raster5, Bound_Raster6) # Process: Reclassify... gp.Reclassify_sa(Bound_Raster6, "VALUE", "0 0;1 0;10 -1;11 -1;2 0;3 1;NODATA 0", Boundary_Raster, "DATA") # Process: Raster to Point... gp.RasterToPoint_conversion(Boundary_Raster, Boundary_Centroids, "VALUE") # Process: Add XY Coordinates... gp.AddXY_management(Boundary_Centroids) # Process: Point to Raster... gp.PointToRaster_conversion(Boundary_Centroids, "POINT_X", Bound_X, "MOST_FREQUENT", "NONE", str(cellsize)) # Process: Point to Raster... gp.PointToRaster_conversion(Boundary_Centroids, "POINT_Y", Bound_Y, "MOST_FREQUENT", "NONE", str(cellsize)) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Boundary_Raster, Boundary_ASCII_file) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Bound_X, X_Coords) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Bound_Y, Y_Coords) # Edit X_COORDS text file MODFLOW_ascii(X_Coords, textFolders + "\\MODFLOW_txt\\X_Coords.txt") # Edit Y_Coords text file MODFLOW_ascii(Y_Coords, textFolders + "\\MODFLOW_txt\\Y_Coords.txt") del thickness,thickness_reclass,Extents_raster,Extents_raster_reclassify,Extents_Polygon del Single_Model_Polygon_Agg,Single_Model_Polygon_Agg_Buffer,Spec_Head_Bound_Poly,Single_Model_Raster del Single_Model_Raster_Agg,Single_Model_Raster_Agg_Buffer,Single_Model_Raster_Reclass del Single_Model_Raster_Agg_Reclass,Single_Model_Raster_Agg_Buffer_Reclass,Bound_Raster4,Bound_Raster5 del Bound_Raster6,Boundary_Raster, Boundary_Centroids,Bound_X,Bound_Y #Create Boundary File for Sloping base if create_SlopingBaseBoundary: gp.AddToolbox(ToolboxFolder + "Conversion Tools.tbx") gp.AddToolbox(ToolboxFolder + "Data Management Tools.tbx") # Local variables... AquiferCoveragePoly = outFolder + "\\Temp\\aquifer_coverage_poly.shp" SlopingBaseExtentsPoly = outFolder + "\\Temp\\sloping_base_poly.shp" SlopingBaseExtentsRaster = outFolder + "\\slp_bse_bound" SlopingBaseExtents_ASCII_file = textFolders + "\\txt\\slopingbasebound.txt" DEM_ASCII_file = textFolders + "\\txt\\dem.txt" # Process: Raster to Polygon... gp.RasterToPolygon_conversion(AquiferCoverage_raster, AquiferCoveragePoly, "NO_SIMPLIFY", "VALUE") # Add Area Field gp.AddField_management (AquiferCoveragePoly, "Area", "FLOAT", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "") # Calculate Area expression = "float(!SHAPE.AREA@SQUAREMILES!)" gp.CalculateField_management (AquiferCoveragePoly, "Area", expression, "PYTHON_9.3") # Find, select, and save largest polygon within boundary to new shapefile: searchRows = gp.searchcursor(AquiferCoveragePoly) searchRow = searchRows.next() Area = [] while searchRow: if searchRow.GRIDCODE == 1: Area.append(searchRow.Area) searchRow = searchRows.next() maxArea = max(Area) area_string = "\"Area\" = " + str(maxArea) area_string = "Area > " + str(int(maxArea) - 1) + " AND Area < " + str(int(maxArea) + 1) gp.MakeFeatureLayer(AquiferCoveragePoly,"bound_lyr") gp.SelectLayerByAttribute_management("bound_lyr", "NEW_SELECTION", area_string) gp.CopyFeatures("bound_lyr", SlopingBaseExtentsPoly) # Process: PolygonToRaster... gp.PolygonToRaster_conversion(SlopingBaseExtentsPoly, "GRIDCODE", SlopingBaseExtentsRaster, "CELL_CENTER", "NONE", str(cellsize)) #Project Raster gp.DefineProjection_management(SlopingBaseExtentsRaster, "PROJCS['NAD_1983_UTM_Zone_14N',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-99.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]") # Process: Raster to ASCII... gp.RasterToASCII_conversion(SlopingBaseExtentsRaster, SlopingBaseExtents_ASCII_file) gp.RasterToASCII_conversion(DEM, DEM_ASCII_file) # Edit text file MODFLOW_ascii(SlopingBaseExtents_ASCII_file, textFolders + "\\MODFLOW_txt\\SlopingBaseBound_MODFLOW.txt") MODFLOW_ascii(DEM_ASCII_file, textFolders + "\\MODFLOW_txt\\DEM_MODFLOW.txt") del AquiferCoveragePoly,SlopingBaseExtentsPoly,SlopingBaseExtentsRaster,SlopingBaseExtents_ASCII_file,AquiferCoverage_raster #Create River Head files if create_RiverHead or create_RiverBottom or create_ModelBoundary or create_Conductance: print 'Creating River files' try: # Check out any necessary licenses gp.CheckOutExtension("spatial") # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Spatial Analyst Tools.tbx") gp.AddToolbox(ToolboxFolder + "Conversion Tools.tbx") gp.AddToolbox(ToolboxFolder + "Data Management Tools.tbx") gp.AddToolbox(ToolboxFolder + "Analysis Tools.tbx") # Local variables... gp.Extent = str(left) + " " + str(bottom) + " " + str(right) + " " + str(top) HighPlainsRivers_Clip = outFolder + "\\HighPlainsRivers_Clip.shp" River_Raster = outFolder + "\\Temp\\River_Raster1" River_Raster_Reclass = outFolder + "\\Temp\\River_Raster2" DEM_Resample = outFolder + "\\Temp\\DEM_Resample" RiverHead_Raster = outFolder + "\\river_head" RiverBottom_Raster = outFolder + "\\Temp\\riverbot_temp" RiverBottom_Raster_Reclass = outFolder + "\\river_bottom" RiverHead_ASCII_file = textFolders + "\\txt\\riverhead.txt" RiverBottom_ASCII_file = textFolders + "\\txt\\riverbottom.txt" # Process: Clip... gp.Clip_analysis(HighPlainsRivers, Single_Model_Polygon, HighPlainsRivers_Clip, "") # Process: Feature to Raster (2)... fields = gp.ListFields(HighPlainsRivers_Clip) names = [] for field in fields: names.append(field.Name) if "IsRiver" in names: gp.FeatureToRaster_conversion(HighPlainsRivers_Clip, "IsRiver", River_Raster, str(cellsize)) else: gp.AddField(HighPlainsRivers_Clip, "IsRiver", "short") gp.CalculateField_management (HighPlainsRivers_Clip, "IsRiver", "1", "PYTHON_9.3") gp.FeatureToRaster_conversion(HighPlainsRivers_Clip, "IsRiver", River_Raster, str(cellsize)) # Process: Reclassify... gp.Reclassify_sa(River_Raster, "VALUE", "1 1;NODATA 0", River_Raster_Reclass, "DATA") # Process: Times... gp.Times_sa(River_Raster_Reclass, DEM, RiverHead_Raster) #Process: Minus... gp.Minus_sa(RiverHead_Raster, "3", RiverBottom_Raster) # Process: Reclassify... gp.Reclassify_sa(RiverBottom_Raster, "VALUE", "-3 0", RiverBottom_Raster_Reclass, "DATA") # Process: Raster to ASCII... gp.RasterToASCII_conversion(RiverHead_Raster, RiverHead_ASCII_file) # Process: Raster to ASCII... gp.RasterToASCII_conversion(RiverBottom_Raster_Reclass, RiverBottom_ASCII_file) del River_Raster,River_Raster_Reclass,DEM_Resample except: raise #Create Conductance files if create_Conductance or create_ModelBoundary or create_RiverHead: print 'Creating Conductance files' try: # Load required toolboxes... gp.AddToolbox(ToolboxFolder + "Conversion Tools.tbx") #Local Variables gp.Extent = str(left) + " " + str(bottom) + " " + str(right) + " " + str(top) Conductance_Raster = outFolder + "\\conductance" Conductance_ASCII_file = textFolders + "\\txt\\conductance.txt" # Process: Feature to Raster... gp.FeatureToRaster_conversion(HighPlainsRivers_Clip, "Conduc", Conductance_Raster, str(cellsize)) # Process: Raster to ASCII... gp.RasterToASCII_conversion(Conductance_Raster, Conductance_ASCII_file) except: raise #Clean up River Head and Conductance text files CleanBoundaryTxt(textFolders + "\\txt\\",textFolders + "\\txt\\","boundary.txt") # Process: Conductance ASCII to Raster... gp.ASCIIToRaster_conversion(Conductance_ASCII_file, Conductance_Raster, "INTEGER") # Process: RiverHead ASCII to Raster... gp.ASCIIToRaster_conversion(RiverHead_ASCII_file, RiverHead_Raster, "INTEGER") # Process: RiverBottom ASCII to Raster... gp.ASCIIToRaster_conversion(RiverBottom_ASCII_file, RiverBottom_Raster_Reclass, "INTEGER") # Edit Boundary text file MODFLOW_ascii(Boundary_ASCII_file, textFolders + "\\MODFLOW_txt\\Boundary_MODFLOW.txt") # Edit Conductance text file MODFLOW_ascii(Conductance_ASCII_file, textFolders + "\\MODFLOW_txt\\Conductance_MODFLOW.txt") # Edit RiverHead text file MODFLOW_ascii(RiverHead_ASCII_file, textFolders + "\\MODFLOW_txt\\RiverHead_MODFLOW.txt") # Edit RiverBottom text file MODFLOW_ascii(RiverBottom_ASCII_file, textFolders + "\\MODFLOW_txt\\RiverBottom_MODFLOW.txt") del Conductance_Raster,Conductance_ASCII_file,RiverHead_Raster,RiverHead_ASCII_file,RiverBottom_Raster,RiverBottom_ASCII_file,HighPlainsRivers_Clip,RiverBottom_Raster_Reclass except: raise #Create Well Files if create_Wells: print 'Creating Well files' try: shutil.copy2(textFolders + "\\MODFLOW_txt\\RiverHead_MODFLOW.txt", textFolders + "\\MODFLOW_txt\\Wells_MODFLOW.txt") except: raise gp.GetMessage(2) finish = time.clock() print 'model run time:',finish-start,'seconds'