I wrote a code and it works as expected for one dataset. When expanding to other datasets, it usually crashes and always returns the following prompt:
Process finished with exit code -1073741819 (0xC0000005)
It doesn't always crash in the same place. I've isolated the part of my code where it crashes after numerous trial and errors--after line 365 at various numbers of iterations through both loops or right after line 378 but before the program terminates normally. Sometimes it crashes after the whole code has executed but before Python tells me it finished without errors. Sometimes it crashes when writing rows to a table (at different parts each time) and other times it crashes when converting those tables to Excel files. I am pretty sure it hasn't crashed when sorting the tables, only when writing (with InsertCursor) or converting them. There is no change in crashing frequency or exit code when I clear all the outputs before running the program again. A Windows notification usually pops up telling me that "Python.exe has stopped working" and then Python terminates and my code stops running and it prints that exit code.
I've tried to use the PDB debugger and post-mortem debugging stuff and I haven't been successful. I have Pycharm and have that debugging tool available to me as well. Honestly, this whole debugging thing is really foreign and I'm not sure what I'm supposed to be looking for. I may not have used it correctly, but the post-mortem debugging function didn't work to tell me about what happened in my code. Step-by-step analysis is super frustrating because it crashes at different points in the code and, again, I don't know what I'm really supposed to be looking for. After some investigation (explained below), I'm not sure a debugger will even help because I think my computer is causing the crash, not Python or my code.
I exchanged some of the files in the working dataset with larger files and it crashes every time now. Switching back to the original dataset does not result in a crash. So, I think it has something to do with the file size, even though those large files are not directly utilized more than once in the code and play no role in writing tables or conversion to Excel. I did try with a really small dataset (smaller than what I used to write the code) and it crashes less often than it completes, but it still crashes, so that's a bit perplexing and tells me that it may not really be related to file size at all.
I suspected it might have to do with a corruption of the output folder somehow, and when I changed my output folder name, the code completed using the large dataset for the first time! Then I ran it again and it crashed. I changed the output location of one of my tables within the output folder and it completed again! But then it crashed when I tried it the next time. Then I changed my output folder location from the C:\ drive to the K:\ drive (portable drive) and it crashed. I changed it back to a new folder I created in the C:\ drive and it crashed again.
Some of the suggested fixes have involved setting a recursion limit or setting a stack limit or going through stacks? I have no idea what that means. I don't even know what a stack is in this context, to be honest. I'm not a computer science person so the inner workings of programming are beyond my current level of understanding (but I am not opposed to having that change).
Googling the error suggests that this is an "Access Violation" that occurs when the computer tries to use memory it shouldn't be using. This computer has a lot of RAM, so that worries me a bit. Does deleting variables once I don't need them in the code anymore free up space? If I have a raster with 1 million pixels but I only need 20,000 pixels, would doing some arcpy operations and overwriting the variable name effectively remove the large dataset from the memory?
I think this problem may be able to be overcome by making my code more efficient, but that's a last-resort effort, honestly. I've streamlined my code as much as I am able and I think I may need some serious help getting my outputs the way I want them using different methods. I really don't feel like this code is particularly taxing on my machine. In fact, I just had the thought to check out what Task manager says about my computer's performance. Interestingly, Python only crashed about 50% of the time while I was looking at the analytics in the Task manager.
I'm running 64-bit Windows 8.1 with 32-bit Python 2.7.8 in Pycharm 4.0.6 and have ArcGIS 10.3 with the intermediate license, I think. I have 24 GB of RAM/Memory, a 2.13 GHz CPU, a 455 GB C:\ drive that's 65% full, and some other drives. I saw that Python (32-bit) was listed twice in the processes for some reason and ended one of them. Python never uses more than 14.5% of my CPU/processor during the whole code whether it crashes or not. My memory/RAM usage never goes above 4.1 GB (or 17%) during the whole code whether it crashes or not and my computer "idles" at 3.7 GB (15%). "Disk 0 (G: C:)" jumps up to 100% several times during the code. When it crashes, it seems to happen most often when Disk usage is 100%. However, it does not always crash when Disk usage is 100%; disk usage reaches 100% or into the 90's at multiple places during my code.
I understand that disk usage has something to do with the read/write speed of my drive. Is my hard drive not powerful enough? Not empty enough? I guess I'll be testing this program from another drive in the meantime to see what happens.
EDIT: I've been doing some more reading before playing around with moving my files to a different drive and I came across this solution for reading/writing problems in u-torrent: Disk Overload 100% solution - Speed Problems - µTorrent Community Forums. Does this sort of solution make sense in the context of Python and is it possible/safe to mess with those sort of settings within my code to get it to run?
EDIT 2: I found the "Resource Monitor" functionality and was looking at the read/write levels of Python.exe while I ran the program. I did not notice any patterns in the reading and writing B/sec usage by Python when it did or did not crash. It looks like Python stays open and active well after the program finishes both when it crashes and does not crash. This is in contrast to the "Processes" tab in the task manager; after crashing or finishing, Python.exe disappears. That's probably normal behavior though.
EDIT 3: I moved my input files and my output location to the D:\ drive, which has 1.6 TB of storage space free. There are no programs or processes set to run from that drive. So, the input and output locations were changed but my PyCharm and Python.exe are still located on the C:\ drive. It crashed the second time I ran the code and though the C:\ drive was busy, it wasn't at 100% and the D:\ drive was even less active.
EDIT 4: I updated PyCharm to 4.5.1 and then migrated both Python and PyCharm to the D:\ drive in addition to my inputs/outputs so a majority of the processes should be taking place from there. It still crashed.
EDIT 5: I thought about getting everything in 64-bit and upgrading to Python 3 but after a bit of digging it looks like ArcGIS 10.X *must* use 32-bit Python 2.7.X, so that's not an option. I upgraded to Python 2.7.9 (32-bit) and ran it again. Still crashed. Then I copied everything to my laptop. My laptop is running 64-bit Windows 8.1, 32-bit PyCharm 4.5.1, and has 32-bit Python 2.7.9. CPU/processor is 2.4 GHz, I have 11.9 GB Memory/RAM and my computer "idles" at 3.3 GB, the C:\ drive has 661 GB free out of 883 GB. I got the code running on my laptop and I actually got a different error and it crashed at a different section twice and the third time crashed at the end and gave the same exit code my desktop is giving. The different laptop exit code happened between lines 295 and 299 at the second iteration through the table both times. The exit code it gave on my laptop was:
Process finished with exit code -1073740940 (0xC0000374)
EDIT 6: My laptop seems to be alternating between crashing with these two exit codes. It may be worth noting that when the variable "totalArea" is used in line 354, PyCharm flags the variable and says "Name 'totalArea' can be not defined. This inspection warns about local variables referenced before assignment."
EDIT 7: It might be worth noting that my laptop has ArcGIS 10.2 while my desktop has ArcGIS 10.3. Also, I've confirmed that the program on my laptop does not crash 100% of the time. I've now managed to run it without incident twice in a row.
EDIT 8: I ended up "solving" the problem by rewriting my code. I figured out that the problem originated with the cursors regardless of where the program crashed. After trying to track where and why that was happening, I realized a better use of my time would be to find some other way to do what I wanted to do. When possible, I converted tables into dictionaries or lists and carried out the InsertCursor and UpdateCursor functions with pure Python. This was actually a good thing because the code completes so much faster now.
My original code is below, as requested.
#### _______________IMPORT MODULES_________________### print("Preparing necessary files...") import os import arcpy import copy ### ______________INPUT FILES___________________### outline = r"D:\Python\Inputs\Lakeridge.gdb\LRoutline" DA = r"D:\Python\Inputs\Lakeridge.gdb\SubAreas" DAID = "DA" # field in DA shapefile where unique DA values exist soil = r"D:\Python\Inputs\VB_SoilsPRJ.shp" WTin = r"D:\Python\Inputs\wtdepth1" DEMin = r"D:\Python\Inputs\2013 5ft DEM.img" MapLoc = r"D:\Python\Inputs\LakeRidge.mxd" WT = arcpy.Raster(WTin) DEM = arcpy.Raster(DEMin) ### ________________SET ENVIRONMENTS__________________### # Check out extension and overwrite outputs arcpy.CheckOutExtension("spatial") arcpy.env.overwriteOutput = True # Set Map Document mxd = arcpy.mapping.MapDocument(MapLoc) # Create project folder and set workspace print("Checking for and creating output folders for spatial data...") WorkPath = MapLoc[:-4] if not os.path.exists(WorkPath): os.makedirs(WorkPath) arcpy.env.workspace = WorkPath # Create scratch workspace ScratchPath = str(WorkPath) + r"\scratch" if not os.path.exists(ScratchPath): os.makedirs(ScratchPath) arcpy.env.scratchWorkspace = ScratchPath # Create GDB path, filename = os.path.split(MapLoc) GDB = filename[:-4] + ".gdb" GDBpath = MapLoc[:-4] + ".gdb" if not os.path.exists(GDBpath): arcpy.CreateFileGDB_management(path, GDB) # Create main output table folder if it does not exist and create project folder print("Checking for and creating output space for Excel files...") TabPath = r"D:\Python\Results" "\\" ProjFolder = TabPath + filename[:-4] if not os.path.exists(TabPath): os.makedirs(TabPath) if not os.path.exists(ProjFolder): os.makedirs(ProjFolder) # Define location of constraint database and establish GIS table output location print("Checking for and creating output space for GIS tables...") CRIT = TabPath + "constraints.xlsx" BMPFold = ProjFolder + r"\GIS-Tables" if not os.path.exists(BMPFold): os.makedirs(BMPFold) ### _________________VERIFY INPUTS________________### # Check that all inputs have the same projection and update list with projected file path names print("Verifying that coordinate systems are the same...") InSHP = [outline, DA, soil] InRAS = [WT] # The base projected coordinate system (PCS) is the DEM's PCS DEMSR = arcpy.Describe(DEM).spatialReference.PCSCode for i, l in enumerate(InSHP): sr = arcpy.Describe(l).spatialReference.PCScode if sr != DEMSR and ".gdb" not in l: l = arcpy.Project_management(l, l[:-4] + "PRJ.shp", DEMSR) InSHP = l elif sr != DEMSR and ".gdb" in l: l = arcpy.Project_management(l, l + "PRJ", DEMSR) InSHP = l sr = arcpy.Describe(WT).spatialReference.PCScode if sr != DEMSR: WTPRJ = arcpy.Raster(arcpy.ProjectRaster_management(WT, "WTPRJ", DEMSR, "CUBIC")) WTPRJ.save(WorkPath + r"\WT_PRJ") WT = WTPRJ # Assign projected file paths to variable names outline = InSHP[0] DA = InSHP[1] soil = InSHP[2] ### _____________SET PROCESSING EXTENTS____________ ### # Set cell size description = arcpy.Describe(DEM) cellsize = description.children[0].meanCellHeight print("Setting cell size to DEM cell size: " + str(cellsize) + " ft...") # Replace ft with code to get units!!! arcpy.env.cellSize = cellsize # Create buffer around outline to use as mask # Buffer distance is in feet print("Creating an environment mask from the site outline shapefile...") maskshp = arcpy.Buffer_analysis(outline, ScratchPath + r"\outline_buff", "50 Feet", "", "", "ALL",) # Convert buffer to raster mask = arcpy.Raster(arcpy.PolygonToRaster_conversion(maskshp, "Id", ScratchPath + r"\rastermask")) mask.save(ScratchPath + r"\rastermask") # Set raster mask and snap raster print("Setting raster mask and snap raster for project...") arcpy.env.mask = mask arcpy.env.snapRaster = mask arcpy.env.extent = mask.extent ### _______________ASSIGN HSG________________### # Many soils in the coastal plain are dual group soils, A/D, B/D, or C/D. # First letter is the drained condition and second letter is the undrained # condition. Soil is considered drained when the depth to water table is # greater than two feet from the surface. # This looks at the HSG assigned to the soil polygon and compares it # to the depth to WT layer. If HSG is unknown or invalid, # HSG is assigned D soil type. # Convert soils shapefile to raster and assign integer values to HSG. # A=1, B=2, C=3, 4=D and dual groups A/D=14, B/D=24, C/D=34 # "---" is treated as a D soil print("Converting dual group soils to single groups...") SoilUnclass = arcpy.PolygonToRaster_conversion(soil, "HSG", ScratchPath + r"\SoilUnclass", "MAXIMUM_COMBINED_AREA") SoilClass = arcpy.sa.Reclassify(SoilUnclass, "HSG", arcpy.sa.RemapValue([["A", 1], ["B", 2], ["C", 3], ["D", 4], ["A/D", 14], ["B/D", 24], ["C/D", 34], ["---", 4]]), "NODATA") SoilClass.save(ScratchPath + r"\HSGraster") # Determine whether locations with dual groups should be considered drained # or undrained and assign a single HSG value to those locations EffHSG = arcpy.sa.Con(SoilClass > 4, arcpy.sa.Con(WT >= 2.0, (SoilClass - 4) / 10, 4), SoilClass) EffHSG.save(WorkPath + r"\EffectiveHSG") ### ______________SUMMARIZE DA PROPERTIES________________ ### # Initialize expression to calculate area of polygons exparea = "float(!shape.area@ACRES!)" # Summarize total area for each DA print("Summarizing DA characteristics...") DAFld = [f.name for f in arcpy.ListFields(DA)] if "Area" not in DAFld: arcpy.AddField_management(DA, "Area", "FLOAT", 6, 3) arcpy.CalculateField_management(DA, "Area", exparea, "PYTHON") stat_field = [["Area", "SUM"]] field_combo = [DAID] DA_area = arcpy.Statistics_analysis(DA, BMPFold + r"\DA_area", stat_field, field_combo) # Convert DA shapefile to raster DAras = arcpy.Raster(arcpy.PolygonToRaster_conversion(DA, DAID, ScratchPath + r"\DAras", "MAXIMUM_AREA")) # Calculate Slope from DEM for the area of interest, convert to integer # and find median slope in each DA slope = arcpy.sa.Slope(DEM, "PERCENT_RISE") slope.save(WorkPath + r"\slope") roundslope = (slope + 0.005) * 100.00 # preserve the last 2 decimal places and round for truncation slopeINT = arcpy.sa.Int(roundslope) # convert to integer by truncation med_slope100 = arcpy.sa.ZonalStatistics(DAras, "VALUE", slopeINT, "MEDIAN", "DATA") # find median (integer operation) med_slope100.save(ScratchPath + r"\intslope") med_slope = med_slope100 / 100.00 # convert back to true median value med_slope.save(WorkPath + r"\medslope") # Find the median depth to water table in each DA rounded to 2 decimal places roundWT = (WT + 0.005) * 100.00 # preserve the last 2 decimal places and round for truncation WTINT = arcpy.sa.Int(roundWT) # convert to integer by truncation med_WT100 = arcpy.sa.ZonalStatistics(DAras, "VALUE", WTINT, "MEDIAN", "DATA") # find median (integer operation) med_WT100.save(ScratchPath + r"\intWT") med_WT = med_WT100 / 100.00 # convert back to true median value med_WT.save(WorkPath + r"\medWT") # Combine rasters to give unique combinations combo = arcpy.sa.Combine([DAras, EffHSG, med_WT100, med_slope100]) combo.save(WorkPath + r"\combo") combo_table = arcpy.BuildRasterAttributeTable_management(combo) # Convert integers to usable format arcpy.AddField_management(combo_table, "HSG", "TEXT", "", "", 6) arcpy.AddField_management(combo_table, "MEDSLOPE", "FLOAT", 5, 2) arcpy.AddField_management(combo_table, "MEDWT", "FLOAT", 5, 2) with arcpy.da.UpdateCursor(combo_table, ["EFFECTIVEHSG", "HSG", "INTSLOPE", "MEDSLOPE", "INTWT", "MEDWT"]) as cursor: for row in cursor: if row[0] == 1: row[1] = "A" if row[0] == 2: row[1] = "B" if row[0] == 3: row[1] = "C" if row[0] == 4: row[1] = "D" row[3] = float(row[2]) / 100.00 row[5] = float(row[4]) / 100.00 cursor.updateRow(row) ### _____________COMPARE CRITERIA_____________________ ### print("Loading constraint database...") # Convert Excel constraint file to GIS table compare = arcpy.ExcelToTable_conversion(CRIT, BMPFold + r"\BMP-constraints") Fields = [f.name for f in arcpy.ListFields(compare)] # Create dictionary from criteria table # Code is the key, other values are stored as a list D = {r[1]:(r[2:]) for r in arcpy.da.SearchCursor(compare, Fields)} # Codes: # SDAB Simple Disconnection A&B # SDCD Simple Disconnection C&D # SDSA Simple Disconnection C&D with Soil Amendments # CAAB Sheet Flow Conservation Area A&B # CACD Sheet Flow Conservation Area C&D # VFA Sheet Flow Veg Filter A # VFSA Sheet Flow Veg Filter B,C&D with Soil Amendments # GCAB Grass Channel A&B # GCCD Grass Channel C&D # GCSA Grass Channel C&D with Soil Amendments # MI1 Micro Infiltration- Level 1 # SI1 Small Infiltration- Level 1 # CI1 Conventional Infiltration- Level 1 # MI2 Micro Infiltration- Level 2 # SI2 Small Infiltration- Level 2 # CI2 Conventional Infiltration- Level 2 # BRE1 Bioretention Basin- Level 1 # BRE2 Bioretention Basin- Level 2 # DS1 Dry Swale- Level 1 # DS2 Dry Swale- Level 2 # WS1 Wet Swale- Level 1 # WS2 Wet Swale- Level 2 # F1 Filter- Level 1 # F2 Filter- Level 2 # CW1 Constructed Wetland- Level 1 # CW2 Constructed Wetland- Level 2 # WP1 Wet Pond- Level 1 # WP2 Wet Pond- Level 2 # WPGW1 Wet Pond with GW- Level 1 # WPGW2 Wet Pond with GW- Level 2 # EDP1 ED Pond- Level 1 # EDP2 ED Pond- Level 2 # Reference: # 0 - BMP # 1 - RR # 2 - PR # 3 - TPR # 4 - NR # 5 - TNR # 6 - SOIL # 7 - MAX_SLOPE # 8 - MIN_CDA # 9 - MAX_CDA # 10 - WT_SEP # 11 - WT_RELAX (boolean) # 12 - COAST_SEP # 13 - MIN_DEPTH # 14 - DEPTH_RELAX (boolean) # 15 - COAST_MIN_DEPTH # 16 - PWOP_PREF # 17 - YEAR_COST # Create output table for BMPs lumped by DA and HSG with criteria table as template Lump = arcpy.CreateTable_management(BMPFold + "\\", "BMP-Allowable") drop = ["OBJECTID", "FIELD1"] arcpy.AddField_management(Lump, "CODE", "TEXT", "", "", 8) arcpy.AddField_management(Lump, "DA", "TEXT", "", "", 15) arcpy.AddField_management(Lump, "HSG", "TEXT", "", "", 6) arcpy.AddField_management(Lump, "BMP", "TEXT", "", "", 50) arcpy.AddField_management(Lump, "MOD", "TEXT", "", "", 25) arcpy.AddField_management(Lump, "RR", "SHORT") arcpy.AddField_management(Lump, "PR", "SHORT") arcpy.AddField_management(Lump, "TPR", "SHORT") arcpy.AddField_management(Lump, "NR", "SHORT") arcpy.AddField_management(Lump, "TNR", "SHORT") arcpy.AddField_management(Lump, "PWOP_PREF", "TEXT", "", "", 25) arcpy.AddField_management(Lump, "YEAR_COST", "TEXT", "", "", 30) arcpy.DeleteField_management(Lump, drop) Fields = [f.name for f in arcpy.ListFields(Lump)] # Create table to build "Rejected BMP" table Fail = arcpy.Copy_management(Lump, BMPFold + "\\" + r"\BMP-Rejected") arcpy.AddField_management(Fail, "RSN_FAILED", "TEXT", "", "", 50) drop = ["BMP", "MOD", "RR", "PR", "TPR", "NR", "TNR", "PWOP_PREF", "YEAR_COST"] arcpy.DeleteField_management(Fail, drop) FFields = [f.name for f in arcpy.ListFields(Fail)] i = 0 print("Comparing site values to constraints...") # Compare the lumped parameters to the constraint dictionary for row in arcpy.da.SearchCursor(combo_table, ["DARAS", "HSG", "MEDSLOPE", "MEDWT"]): i += 1 print i # Temporarily store the area of each DA for later comparison print("Compare total area") for r in arcpy.da.SearchCursor(DA_area, [DAID, "SUM_AREA"]): # NOTE: DAID *must* be an integer value to function properly. if r[0] == row[0]: totalArea = r[1] # Duplicate criteria dictionary that can be amended throughout the loop print("Copy BMP") BMP = copy.deepcopy(D) # Initialize empty dictionary to store BMPs that fail each test print("Initialize empty dictionaries") NoBMP = {} Mod = {} # Compare lumped values in each DA/HSG pair to those in the constraint table print("Begin dictionary loop") for k, v in D.items(): # Test if soil type is incorrect for each BMP and store reason for failure print("Soil test") if row[1] not in v[6]: NoBMP= "Soil type mismatch" # Compare median slope to maximum slope if row[2] > v[7]: if k not in NoBMP.keys(): NoBMP= "Slope too steep" else: NoBMP+= ", Slope too steep" # Compare WT depths print("WT depth test") if v[10] == 0: Mod= "---" elif v[13] + v[10] <= row[3]: Mod= "---" elif v[13] + v[10] > row[3]: # Check if coastal modification allows use of practice if v[11] == 1: coast_WT = v[12] else: coast_WT = v[10] if v[14] == 1: coast_depth = v[15] else: coast_depth = v[13] # Notate if coastal modification allows for practice use if coast_WT + coast_depth <= row[3]: if v[11] == 1 and v[14] == 1: Mod= "Separation and Depth" elif v[11] == 1: Mod= "WT Separation" elif v[14] == 1: Mod= "Practice Depth" else: Mod= "---" # Remove the practice if coastal modifications do not help if coast_WT + coast_depth > row[3]: if k not in NoBMP.keys(): NoBMP= "WT proximity" else: NoBMP+= ", WT proximity" # Compare allowable contributing drainage areas (in acres) # Maximum CDA neglected because this is lumped analysis print("Compare areas") if v[8] >= totalArea: if k not in NoBMP.keys(): NoBMP= "CDA too small" else: NoBMP+= ", CDA too small" # Compare keys in BMP and NoBMP dictionaries. Remove matching pairs from the BMP dictionary. print("Removing bad BMPs from the BMP dictionary") for key in BMP.keys(): if key in NoBMP.keys(): del BMP[key] # Write remaining BMPs to table print("Writing BMPs to output table") with arcpy.da.InsertCursor(Lump, Fields) as cursor: for k,v in BMP.items(): cursor.insertRow((0, k, row[0], row[1], v[0], Mod, v[1], v[2], v[3], v[4], v[5], v[16], v[17])) # Sort values in table, effectively ranking them print("Ranking BMPs...") LumpSort = arcpy.Sort_management(Lump, BMPFold + "\\LumpSort", [["DA", "ASCENDING"], ["HSG", "ASCENDING"], ["TPR", "DESCENDING"]]) arcpy.DeleteField_management(LumpSort, ["ROWID"]) # Convert tables to readable format outside of GIS (.xls) print("Converting good BMPs to Excel format...") arcpy.TableToExcel_conversion(LumpSort, ProjFolder + r"\Lumped-Result.xls")
Can you post the code here so we can review it?
Absolutely. I updated the main post with another thing I tried and also my code.
Have your considered partitioning your code into separate blocks? So far, I see no reason that the code has to run from start to finish. Is this a matter of convenience? If split and constructed to work in a toolbox, there are logical tools that could be created to subdivide your workflow. This of course would add a few mouse clicks, but defined defaults in the tool settings, it might just allow the subprocess to finish clearing up memory and identifying errors etc in blocks. You have indicated in your last update that it has worked a few times on the new machine, but I am not sure you have identified what changes you had to make to the program to get it to run.
Thank you for your response! I will actually be running this as a standalone program, not as a script from GIS. I don't know how to partition a code into blocks that can run separately.
I actually don't remember if I made any changes to the code to get it to run on my laptop or if it just stopped throwing that other error. I feel like I changed how "totalarea" was assigned but that version of the code is long gone, so I can't share it now.
I ended up "solving" the problem by rewriting my code. I figured out that the problem originated with the cursors regardless of where the program crashed. After trying to track where and why that was happening, I realized a better use of my time would be to find some other way to do what I wanted to do. When possible, I converted tables into dictionaries or lists and carried out the InsertCursor and UpdateCursor functions with pure Python. This was actually a good thing because the code completes so much faster now.
With regards to Edit6, the variable:
totalArea
Appears to come from the code here:
totalArea = 0
#### _______________IMPORT MODULES_________________### print("Preparing necessary files...") import os import arcpy import copy ### ______________INPUT FILES___________________### outline = r"D:\Python\Inputs\Lakeridge.gdb\LRoutline" DA = r"D:\Python\Inputs\Lakeridge.gdb\SubAreas" DAID = "DA" # field in DA shapefile where unique DA values exist soil = r"D:\Python\Inputs\VB_SoilsPRJ.shp" WTin = r"D:\Python\Inputs\wtdepth1" DEMin = r"D:\Python\Inputs\2013 5ft DEM.img" MapLoc = r"D:\Python\Inputs\LakeRidge.mxd" WT = arcpy.Raster(WTin) DEM = arcpy.Raster(DEMin) ### ________________SET ENVIRONMENTS__________________### # Check out extension and overwrite outputs arcpy.CheckOutExtension("spatial") arcpy.env.overwriteOutput = True # Set Map Document mxd = arcpy.mapping.MapDocument(MapLoc) # Create project folder and set workspace print("Checking for and creating output folders for spatial data...") WorkPath = MapLoc[:-4] if not os.path.exists(WorkPath): os.makedirs(WorkPath) arcpy.env.workspace = WorkPath # Create scratch workspace ScratchPath = str(WorkPath) + r"\scratch" if not os.path.exists(ScratchPath): os.makedirs(ScratchPath) arcpy.env.scratchWorkspace = ScratchPath # Create GDB path, filename = os.path.split(MapLoc) GDB = filename[:-4] + ".gdb" GDBpath = MapLoc[:-4] + ".gdb" if not os.path.exists(GDBpath): arcpy.CreateFileGDB_management(path, GDB) # Create main output table folder if it does not exist and create project folder print("Checking for and creating output space for Excel files...") TabPath = r"D:\Python\Results" "\\" ProjFolder = TabPath + filename[:-4] if not os.path.exists(TabPath): os.makedirs(TabPath) if not os.path.exists(ProjFolder): os.makedirs(ProjFolder) # Define location of constraint database and establish GIS table output location print("Checking for and creating output space for GIS tables...") CRIT = TabPath + "constraints.xlsx" BMPFold = ProjFolder + r"\GIS-Tables" if not os.path.exists(BMPFold): os.makedirs(BMPFold) ### _________________VERIFY INPUTS________________### # Check that all inputs have the same projection and update list with projected file path names print("Verifying that coordinate systems are the same...") InSHP = [outline, DA, soil] InRAS = [WT] # The base projected coordinate system (PCS) is the DEM's PCS DEMSR = arcpy.Describe(DEM).spatialReference.PCSCode for i, l in enumerate(InSHP): sr = arcpy.Describe(l).spatialReference.PCScode if sr != DEMSR and ".gdb" not in l: l = arcpy.Project_management(l, l[:-4] + "PRJ.shp", DEMSR) InSHP = l elif sr != DEMSR and ".gdb" in l: l = arcpy.Project_management(l, l + "PRJ", DEMSR) InSHP = l sr = arcpy.Describe(WT).spatialReference.PCScode if sr != DEMSR: WTPRJ = arcpy.Raster(arcpy.ProjectRaster_management(WT, "WTPRJ", DEMSR, "CUBIC")) WTPRJ.save(WorkPath + r"\WT_PRJ") WT = WTPRJ # Assign projected file paths to variable names outline = InSHP[0] DA = InSHP[1] soil = InSHP[2] ### _____________SET PROCESSING EXTENTS____________ ### # Set cell size description = arcpy.Describe(DEM) cellsize = description.children[0].meanCellHeight print("Setting cell size to DEM cell size: " + str(cellsize) + " ft...") # Replace ft with code to get units!!! arcpy.env.cellSize = cellsize # Create buffer around outline to use as mask # Buffer distance is in feet print("Creating an environment mask from the site outline shapefile...") maskshp = arcpy.Buffer_analysis(outline, ScratchPath + r"\outline_buff", "50 Feet", "", "", "ALL",) # Convert buffer to raster mask = arcpy.Raster(arcpy.PolygonToRaster_conversion(maskshp, "Id", ScratchPath + r"\rastermask")) mask.save(ScratchPath + r"\rastermask") # Set raster mask and snap raster print("Setting raster mask and snap raster for project...") arcpy.env.mask = mask arcpy.env.snapRaster = mask arcpy.env.extent = mask.extent ### _______________ASSIGN HSG________________### # Many soils in the coastal plain are dual group soils, A/D, B/D, or C/D. # First letter is the drained condition and second letter is the undrained # condition. Soil is considered drained when the depth to water table is # greater than two feet from the surface. # This looks at the HSG assigned to the soil polygon and compares it # to the depth to WT layer. If HSG is unknown or invalid, # HSG is assigned D soil type. # Convert soils shapefile to raster and assign integer values to HSG. # A=1, B=2, C=3, 4=D and dual groups A/D=14, B/D=24, C/D=34 # "---" is treated as a D soil print("Converting dual group soils to single groups...") SoilUnclass = arcpy.PolygonToRaster_conversion(soil, "HSG", ScratchPath + r"\SoilUnclass", "MAXIMUM_COMBINED_AREA") SoilClass = arcpy.sa.Reclassify(SoilUnclass, "HSG", arcpy.sa.RemapValue([["A", 1], ["B", 2], ["C", 3], ["D", 4], ["A/D", 14], ["B/D", 24], ["C/D", 34], ["---", 4]]), "NODATA") SoilClass.save(ScratchPath + r"\HSGraster") # Determine whether locations with dual groups should be considered drained # or undrained and assign a single HSG value to those locations EffHSG = arcpy.sa.Con(SoilClass > 4, arcpy.sa.Con(WT >= 2.0, (SoilClass - 4) / 10, 4), SoilClass) EffHSG.save(WorkPath + r"\EffectiveHSG") ### ______________SUMMARIZE DA PROPERTIES________________ ### # Initialize expression to calculate area of polygons exparea = "float(!shape.area@ACRES!)" # Summarize total area for each DA print("Summarizing DA characteristics...") DAFld = [f.name for f in arcpy.ListFields(DA)] if "Area" not in DAFld: arcpy.AddField_management(DA, "Area", "FLOAT", 6, 3) arcpy.CalculateField_management(DA, "Area", exparea, "PYTHON") stat_field = [["Area", "SUM"]] field_combo = [DAID] DA_area = arcpy.Statistics_analysis(DA, BMPFold + r"\DA_area", stat_field, field_combo) # Convert DA shapefile to raster DAras = arcpy.Raster(arcpy.PolygonToRaster_conversion(DA, DAID, ScratchPath + r"\DAras", "MAXIMUM_AREA")) # Calculate Slope from DEM for the area of interest, convert to integer # and find median slope in each DA slope = arcpy.sa.Slope(DEM, "PERCENT_RISE") slope.save(WorkPath + r"\slope") roundslope = (slope + 0.005) * 100.00 # preserve the last 2 decimal places and round for truncation slopeINT = arcpy.sa.Int(roundslope) # convert to integer by truncation med_slope100 = arcpy.sa.ZonalStatistics(DAras, "VALUE", slopeINT, "MEDIAN", "DATA") # find median (integer operation) med_slope100.save(ScratchPath + r"\intslope") med_slope = med_slope100 / 100.00 # convert back to true median value med_slope.save(WorkPath + r"\medslope") # Find the median depth to water table in each DA rounded to 2 decimal places roundWT = (WT + 0.005) * 100.00 # preserve the last 2 decimal places and round for truncation WTINT = arcpy.sa.Int(roundWT) # convert to integer by truncation med_WT100 = arcpy.sa.ZonalStatistics(DAras, "VALUE", WTINT, "MEDIAN", "DATA") # find median (integer operation) med_WT100.save(ScratchPath + r"\intWT") med_WT = med_WT100 / 100.00 # convert back to true median value med_WT.save(WorkPath + r"\medWT") # Combine rasters to give unique combinations combo = arcpy.sa.Combine([DAras, EffHSG, med_WT100, med_slope100]) combo.save(WorkPath + r"\combo") combo_table = arcpy.BuildRasterAttributeTable_management(combo) # Convert integers to usable format arcpy.AddField_management(combo_table, "HSG", "TEXT", "", "", 6) arcpy.AddField_management(combo_table, "MEDSLOPE", "FLOAT", 5, 2) arcpy.AddField_management(combo_table, "MEDWT", "FLOAT", 5, 2) with arcpy.da.UpdateCursor(combo_table, ["EFFECTIVEHSG", "HSG", "INTSLOPE", "MEDSLOPE", "INTWT", "MEDWT"]) as cursor: for row in cursor: if row[0] == 1: row[1] = "A" if row[0] == 2: row[1] = "B" if row[0] == 3: row[1] = "C" if row[0] == 4: row[1] = "D" row[3] = float(row[2]) / 100.00 row[5] = float(row[4]) / 100.00 cursor.updateRow(row) ### _____________COMPARE CRITERIA_____________________ ### print("Loading constraint database...") # Convert Excel constraint file to GIS table compare = arcpy.ExcelToTable_conversion(CRIT, BMPFold + r"\BMP-constraints") Fields = [f.name for f in arcpy.ListFields(compare)] # Create dictionary from criteria table # Code is the key, other values are stored as a list D = {r[1]:(r[2:]) for r in arcpy.da.SearchCursor(compare, Fields)} # Codes: # SDAB Simple Disconnection A&B # SDCD Simple Disconnection C&D # SDSA Simple Disconnection C&D with Soil Amendments # CAAB Sheet Flow Conservation Area A&B # CACD Sheet Flow Conservation Area C&D # VFA Sheet Flow Veg Filter A # VFSA Sheet Flow Veg Filter B,C&D with Soil Amendments # GCAB Grass Channel A&B # GCCD Grass Channel C&D # GCSA Grass Channel C&D with Soil Amendments # MI1 Micro Infiltration- Level 1 # SI1 Small Infiltration- Level 1 # CI1 Conventional Infiltration- Level 1 # MI2 Micro Infiltration- Level 2 # SI2 Small Infiltration- Level 2 # CI2 Conventional Infiltration- Level 2 # BRE1 Bioretention Basin- Level 1 # BRE2 Bioretention Basin- Level 2 # DS1 Dry Swale- Level 1 # DS2 Dry Swale- Level 2 # WS1 Wet Swale- Level 1 # WS2 Wet Swale- Level 2 # F1 Filter- Level 1 # F2 Filter- Level 2 # CW1 Constructed Wetland- Level 1 # CW2 Constructed Wetland- Level 2 # WP1 Wet Pond- Level 1 # WP2 Wet Pond- Level 2 # WPGW1 Wet Pond with GW- Level 1 # WPGW2 Wet Pond with GW- Level 2 # EDP1 ED Pond- Level 1 # EDP2 ED Pond- Level 2 # Reference: # 0 - BMP # 1 - RR # 2 - PR # 3 - TPR # 4 - NR # 5 - TNR # 6 - SOIL # 7 - MAX_SLOPE # 8 - MIN_CDA # 9 - MAX_CDA # 10 - WT_SEP # 11 - WT_RELAX (boolean) # 12 - COAST_SEP # 13 - MIN_DEPTH # 14 - DEPTH_RELAX (boolean) # 15 - COAST_MIN_DEPTH # 16 - PWOP_PREF # 17 - YEAR_COST # Create output table for BMPs lumped by DA and HSG with criteria table as template Lump = arcpy.CreateTable_management(BMPFold + "\\", "BMP-Allowable") drop = ["OBJECTID", "FIELD1"] arcpy.AddField_management(Lump, "CODE", "TEXT", "", "", 8) arcpy.AddField_management(Lump, "DA", "TEXT", "", "", 15) arcpy.AddField_management(Lump, "HSG", "TEXT", "", "", 6) arcpy.AddField_management(Lump, "BMP", "TEXT", "", "", 50) arcpy.AddField_management(Lump, "MOD", "TEXT", "", "", 25) arcpy.AddField_management(Lump, "RR", "SHORT") arcpy.AddField_management(Lump, "PR", "SHORT") arcpy.AddField_management(Lump, "TPR", "SHORT") arcpy.AddField_management(Lump, "NR", "SHORT") arcpy.AddField_management(Lump, "TNR", "SHORT") arcpy.AddField_management(Lump, "PWOP_PREF", "TEXT", "", "", 25) arcpy.AddField_management(Lump, "YEAR_COST", "TEXT", "", "", 30) arcpy.DeleteField_management(Lump, drop) Fields = [f.name for f in arcpy.ListFields(Lump)] # Create table to build "Rejected BMP" table Fail = arcpy.Copy_management(Lump, BMPFold + "\\" + r"\BMP-Rejected") arcpy.AddField_management(Fail, "RSN_FAILED", "TEXT", "", "", 50) drop = ["BMP", "MOD", "RR", "PR", "TPR", "NR", "TNR", "PWOP_PREF", "YEAR_COST"] arcpy.DeleteField_management(Fail, drop) FFields = [f.name for f in arcpy.ListFields(Fail)] i = 0 print("Comparing site values to constraints...") # Compare the lumped parameters to the constraint dictionary for row in arcpy.da.SearchCursor(combo_table, ["DARAS", "HSG", "MEDSLOPE", "MEDWT"]): i += 1 print i # Temporarily store the area of each DA for later comparison print("Compare total area") for r in arcpy.da.SearchCursor(DA_area, [DAID, "SUM_AREA"]): # NOTE: DAID *must* be an integer value to function properly. if r[0] == row[0]: totalArea = r[1] # Duplicate criteria dictionary that can be amended throughout the loop print("Copy BMP") BMP = copy.deepcopy(D) # Initialize empty dictionary to store BMPs that fail each test print("Initialize empty dictionaries") NoBMP = {} Mod = {} # Compare lumped values in each DA/HSG pair to those in the constraint table print("Begin dictionary loop") for k, v in D.items(): # Test if soil type is incorrect for each BMP and store reason for failure print("Soil test") if row[1] not in v[6]: NoBMP= "Soil type mismatch" # Compare median slope to maximum slope if row[2] > v[7]: if k not in NoBMP.keys(): NoBMP = "Slope too steep" else: NoBMP += ", Slope too steep" # Compare WT depths print("WT depth test") if v[10] == 0: Mod = "---" elif v[13] + v[10] <= row[3]: Mod = "---" elif v[13] + v[10] > row[3]: # Check if coastal modification allows use of practice if v[11] == 1: coast_WT = v[12] else: coast_WT = v[10] if v[14] == 1: coast_depth = v[15] else: coast_depth = v[13] # Notate if coastal modification allows for practice use if coast_WT + coast_depth <= row[3]: if v[11] == 1 and v[14] == 1: Mod = "Separation and Depth" elif v[11] == 1: Mod = "WT Separation" elif v[14] == 1: Mod = "Practice Depth" else: Mod = "---" # Remove the practice if coastal modifications do not help if coast_WT + coast_depth > row[3]: if k not in NoBMP.keys(): NoBMP = "WT proximity" else: NoBMP += ", WT proximity" # Compare allowable contributing drainage areas (in acres) # Maximum CDA neglected because this is lumped analysis print("Compare areas") if v[8] >= totalArea: if k not in NoBMP.keys(): NoBMP = "CDA too small" else: NoBMP += ", CDA too small" # Compare keys in BMP and NoBMP dictionaries. Remove matching pairs from the BMP dictionary. print("Removing bad BMPs from the BMP dictionary") for key in BMP.keys(): if key in NoBMP.keys(): del BMP[key] # Write remaining BMPs to table print("Writing BMPs to output table") with arcpy.da.InsertCursor(Lump, Fields) as cursor: for k,v in BMP.items(): cursor.insertRow((0, k, row[0], row[1], v[0], Mod , v[1], v[2], v[3], v[4], v[5], v[16], v[17])) # Sort values in table, effectively ranking them print("Ranking BMPs...") LumpSort = arcpy.Sort_management(Lump, BMPFold + "\\LumpSort", [["DA", "ASCENDING"], ["HSG", "ASCENDING"], ["TPR", "DESCENDING"]]) arcpy.DeleteField_management(LumpSort, ["ROWID"]) # Convert tables to readable format outside of GIS (.xls) print("Converting good BMPs to Excel format...") arcpy.TableToExcel_conversion(LumpSort, ProjFolder + r"\Lumped-Result.xls")
The local variable referenced before assigned, is flagged as this part of the code may never get to line 11 in my quote, due to the if statements.
This may never happen, but pycharm cant see your data so flags it as a possible error as when you do the check in:
if v[8] >= totalArea:
Depending on your data, totalArea may not exist yet!
You could avoid this being flagged by assigning a default value for this variable as I have sort of shown above.
Thank you! This was very helpful. It is so good to know why it was happening and why the changes I made solved the problem. I ended up changing how that variable was used, amongst many other things.