POST
|
I searched and there's one question similar to mine, but it didn't have any responses. So, I'm trying again. I am using ArcGIS 10.3. I have a parcel layer that doesn't include the roads. I need parcels bordering empty space to extend to fill in the gaps, sharing the area evenly with opposite parcels. I provided an example of what I need to accomplish. I did the work in the example by hand, but I have a whole city to do. Is there a built-in tool that does something similar? Or a series of steps that can approximate the desired result? Thanks for your time.
... View more
06-09-2017
09:27 AM
|
1
|
2
|
903
|
POST
|
I recently finished a stand alone Python program that uses Arcpy extensively. I wanted to package and distribute it as an executable so my users would have minimal problems trying to get the program working. My program uses third-party libraries and I needed to update numpy on all the computers I tried to use it on. The Tcl/Tk package was out of date on one of the computers and trying to update it was a nightmare. Even though people using my program already have Python because they have ArcGIS, I really needed to get my program into an executable so my users would not have to endure the headache of setting up their Python environment to run my program. Cue Py2exe. I set it up, ran my .exe, and got an error immediately: "arcpy module could not be found" or something to that effect. Well, crap. A bit of Internet searching did not turn up the results I wanted. There were some complicated workarounds, such as calling the arcpy module in a subprocess, creating and using a wrapper script, or rewrite the script entirely. Someone even said that it just can't be done. Well, it can, and the solution is actually very simple, elegant (in my opinion), and requires very little extra effort by the user and the developer. Arcpy isn't packaged with the rest of the files when the executable is made, and forcing it to package the module is sure to be a gross violation of the ArcGIS user agreement. Every user of the program will already have the arcpy module somewhere on their computer, so to get the executable to work, we just need to tell the program where to look for the arcpy module. Unfortunately, the location of the arcpy module is likely to be different on everyone's computer. Fortunately, the arcpy module isn't located with the rest of the Python modules anyway! Instead, a .PTH file tells the program where to find the module. We can make this work for us with Py2exe. 1. Exclude "arcpy" from the Py2exe setup.py script. I'm not 100% sure this is necessary, to be honest. There's no harm in doing it though. options = {"py2exe": {"excludes": ["arcpy"]}}
setup(console=['MyProgram.py'], options=options) 2. Add the executable directory + "\\site-packages" to the Python site directory in the main code of MyProgram.py. sys.executable gives the file path of the Python interpreter. In this case, the Python interpreter is the executable. The .PTH file needs to go into a folder called "site-packages" in order for the arcpy module to import correctly. (You will add the site-packages folder manually to the dist folder after Py2exe has completed its job.) If you're not using Windows, change the backslashes to a single forward slash, ie + "/site-packages". I also added the code to any script that would be running as a separate process, just in case. from site import addsitedir
from sys import executable
from os import path
interpreter = executable
sitepkg = path.dirname(interpreter) + "\\site-packages"
addsitedir(sitepkg) 3. Run Py2exe as usual. 4. Add a "site-packages" folder to the "dist" folder. 5. Navigate to the GIS-installed Python folder, usually C:/Python27/ArcGIS10.x/Lib/site-packages 6. Copy the "Desktop10.x" PTH file and paste it into the "site-packages" folder in the "dist" directory for your application. That's it. The executable should load the arcpy module now. Instruct your users to copy the "Desktop10.x" PTH file into the executable's "site-packages" folder prior to running the application for the first time and it should work for them, too.
... View more
08-19-2015
08:14 AM
|
2
|
9
|
12049
|
POST
|
I'm using 10.3 with the latest version of Python 2.
... View more
07-24-2015
03:19 PM
|
0
|
1
|
1062
|
POST
|
Oh no, I thought this would work but that module only works when the code is being used within GIS. Still, it's very helpful for anyone who isn't making a standalone script. I wish I could set your answer to "90% correct answer" because in most situations, it would be the way to solve my problem!
... View more
07-24-2015
12:39 PM
|
0
|
3
|
1062
|
POST
|
I'm building a GUI for users to input their data with Tkinter and TkFileDialog. When I use TkFileDialog, it gives me a list of the individual files that Windows sees when I am trying to pick my file (top image), but I really want it to merge those files together to form one cohesive file like ArcCatalog does automatically (bottom image). I'm pretty sure my users are going to be up in arms about what file to pick if I don't make an effort to change it. I know I can specify the extensions that show and that is a reasonable workaround. However, a user is completely unable to select files from a .GDB in this way. I can think of a workaround for this but I would rather not have to do that as it is very labor intensive and makes me add a lot of extra parts to my GUI. I found IGxObjectFilter but that appears to not be a Python thing. Is there a way to read the files cohesively within Python?
... View more
07-23-2015
09:41 AM
|
0
|
5
|
3811
|
POST
|
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.
... View more
07-17-2015
08:07 AM
|
0
|
0
|
23926
|
POST
|
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.
... View more
07-17-2015
08:01 AM
|
0
|
0
|
23926
|
POST
|
These test runs have all been running on a 15 acre site split into 7 subareas with 2 soil types; the outside dictionary loop for this set of data only has 14 iterations. Dictionary D has around 30 key:value pairs, so at maximum, the final table will have 30 entries for each of the 14 rows, or 420 entries. There's really only about 250 rows in reality for this particular site. However, three of my inputs (DEM, WT, and soil) have spatial information for a whole city. The DEM and WT are rasters with cell size of 5 and are 2.4 GB each. I can't seem to find how big the soil shapefile is but it's not any bigger than 2.4 GB. Yes, ultimately if I choose to go with a direct-to-CSV/Excel route, I'd be storing an intermediate Python object. I was trying to figure out how to store everything as a dictionary since the key values would duplicate for each row and overwrite themselves, but the list of lists solves that issue. So, thank you! EDIT: Although, maybe I can try storing the intermediate and then using the cursor outside of the initial dictionary loop. EDIT 2: I changed the end of my code and it still crashes after converting stuff to Excel. But I do think it's going a lot faster, at least. I may have to bit the bullet and bring in an outside Excel writer. Here is the end of my code now. # Compare keys in BMP and NoBMP dictionaries. Remove matching pairs from the BMP dictionary.
for k in BMP.keys():
if k in NoBMP.keys():
del BMP
# Add remaining BMPs to an intermediate list of lists
for k, v in BMP.items():
BMPs.append([0, k, value[0], value[1], v[0], Mod , v[1], v[2], v[3], v[4], v[5], v[16], v[17]])
# Write BMP list to GIS table
cursor = arcpy.da.InsertCursor(Lump, Fields)
for l in BMPs:
cursor.insertRow((l[0], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8], l[9], l[10], l[11], l[12]))
del cursor
# 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")
... View more
06-15-2015
10:16 AM
|
0
|
0
|
992
|
POST
|
That's interesting. Is the "with" command treated differently from saying "cursor = X" or is that a preference/conventino/speed thing? This is my new code for that section: # Write remaining BMPs to table
cursor = arcpy.da.InsertCursor(Lump, Fields)
for k, v in BMP.items():
cursor.insertRow((0, k, value[0], value[1], v[0], Mod , v[1], v[2], v[3], v[4], v[5], v[16], v[17]))
del cursor It still crashed, unfortunately. But this time it crashed between lines 374 and 379, which is the first time that's happened I think. Hmm!
... View more
06-15-2015
09:56 AM
|
0
|
0
|
992
|
POST
|
This is kind of an update to my last post: Python Crashing: Process finished with exit code -1073741819 (0xC0000005) I understand a bit more *what* is happening to cause my code to crash, but I'm fuzzy on the details of why because my CompSci-Fu is very weak. Apparently, there is some problem with my code forgetting to ask for read/write access, and this is happening due to some problem within the underlying C code of the arcpy functions, most likely. I have PyCharm debugger at my disposal but I am not really sure what to be looking for in the debugger to figure out what I need to change to make my code happy. I'd really like to *not* change my code but I'm not sure I have an option unless I have some professional intervention. My code is crashing at two places: sometime between lines 368 and 371 at different numbers of iterations through the outside loop OR it's crashing after line 381 when the table has been converted to Excel and saved, but before the code finishes by itself. It crashes with the following notification in the console: Process finished with exit code -1073741819 (0xC0000005) Do you have any suggestions for me? I don't NEED it to have an intermediate step as a GIS table; I can write it directly to a CSV or Excel, I guess. I am concerned because the values I am writing contain commas and I don't want that to screw up the CSV. I also am hesitant to introduce a third party library for code upkeep reasons. Basically, I'm looking for suggestions on how to improve my code to avoid the error I am getting, knowing that the error most likely originates in the arcpy.da.InsertCursor or arcpy.TableToExcel_conversion modules, with preference for built-in Python functions or libraries. #### _______________IMPORT MODULES_________________###
print("Preparing necessary files...")
import os
import arcpy
import copy
### ______________INPUT FILES___________________###
outline = r"D:\Python\Inputs\LakeRidge\Lakeridge.gdb\LRoutline"
DA = r"D:\Python\Inputs\LakeRidge\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\LRpythontest.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", 7, 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 area lookup table to dictionary
arcpy.AddField_management(DA_area, "T_AREA", "FLOAT", 7, 3)
with arcpy.da.UpdateCursor(DA_area, ["SUM_AREA", "T_AREA"]) as cursor:
for r in cursor:
r[0] *= 100.00
r[0] = int(r[0])
r[0] = float(r[0])/100.00
r[1] = r[0]
print r[1]
cursor.updateRow(r)
DA_area = {r[0]: (r[1]) for r in arcpy.da.SearchCursor(DA_area, [DAID, "T_AREA"])}
# 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)
arcpy.AddField_management(combo_table, "T_AREA", "FLOAT", 7, 3)
with arcpy.da.UpdateCursor(combo_table, ["EFFECTIVEHSG", "HSG", "INTSLOPE", "MEDSLOPE", "INTWT", "MEDWT", "DARAS", "T_AREA"]) 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
for k, v in DA_area.items():
if row[6] == k:
row[7] = v
cursor.updateRow(row)
# Create dictionary for the DA information
DA_summary = {r[0]: (r[1:]) for r in arcpy.da.SearchCursor(combo_table, ["Rowid", "DARAS", "HSG", "MEDSLOPE", "MEDWT", "T_AREA"])}
### _____________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)]
print("Comparing site values to constraints...")
# Compare the lumped parameters to the constraint dictionary
for key, value in DA_summary.items():
# Duplicate criteria dictionary that can be amended throughout the loop
BMP = copy.deepcopy(D)
# Initialize empty dictionary to store BMPs that fail each test
NoBMP = {}
Mod = {}
# Compare lumped values in each DA/HSG pair to those in the constraint table
for k, v in D.items():
# Test if soil type is incorrect for each BMP and store reason for failure
if value[1] not in v[6]:
NoBMP = "Soil type mismatch"
# Compare median slope to maximum slope
if value[2] > v[7]:
if k not in NoBMP.keys():
NoBMP = "Slope too steep"
else:
NoBMP += ", Slope too steep"
# Compare WT depths
if v[10] == 0:
Mod = "---"
elif v[13] + v[10] <= value[3]:
Mod = "---"
elif v[13] + v[10] > value[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 <= value[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 > value[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
if v[8] >= value[4]:
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.
for k in BMP.keys():
if k in NoBMP.keys():
del BMP
# Write remaining BMPs to table
with arcpy.da.InsertCursor(Lump, Fields) as cursor:
for k, v in BMP.items():
cursor.insertRow((0, k, value[0], value[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")
... View more
06-15-2015
09:21 AM
|
0
|
5
|
4365
|
POST
|
Of course! I can't try this out yet but I bet this is the solution!
... View more
06-12-2015
06:28 AM
|
0
|
0
|
2543
|
POST
|
Absolutely. I updated the main post with another thing I tried and also my code.
... View more
06-11-2015
02:19 PM
|
0
|
2
|
23926
|
POST
|
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")
... View more
06-11-2015
09:45 AM
|
0
|
6
|
54212
|
POST
|
These are my original inputs: outline = r"K:\GradWork\GIS\CollegePark\CP_outline.shp"
DA = r"K:\GradWork\GIS\CollegePark\CPDrainage.shp"
DAID = "Id" # field in DA shapefile where unique DA values exist
soil = r"C:\Users\Rachael Johnson\Documents\GradWork\Project Files\BMPGIS\soils\VB_Soils.shp"
WTin = r"K:\GradWork\GIS\CollegePark\CPGW_adj1"
DEMin = r"K:\GradWork\GIS\CollegePark\cpelev"
MapLoc = r"K:\GradWork\GIS\CollegePark\CP.mxd"
WT = arcpy.Raster(WTin)
DEM = arcpy.Raster(DEMin) My code uses these inputs in the context of other arcpy functions to create new variables or set environments. The outline variable, for example, is used here and then isn't used again: ### _____________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...")
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 Here's an example of using one of the inputs to create a new variable and how that new variable is later transformed. # 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") Project_management does return a variable that I can use in that way; with most of these arcpy functions, setting a variable to them and then calling the variable prints the file path of the resulting shapefile or grid. So, in the context of the above code, "print SoilClass" returns the path "K:\GradWork\GIS\CollegePark\CP\scratch\HSGraster".
... View more
06-10-2015
09:02 AM
|
1
|
1
|
2543
|
POST
|
Is that code you posted preferable to my current code using enumerate? Your code still requires me to update the variables manually outside of the loop, correct?
... View more
06-10-2015
08:43 AM
|
0
|
0
|
2543
|
Title | Kudos | Posted |
---|---|---|
1 | 06-09-2017 09:27 AM | |
1 | 06-10-2015 09:02 AM | |
1 | 04-21-2015 06:57 AM | |
2 | 08-19-2015 08:14 AM | |
2 | 06-09-2015 11:35 AM |
Online Status |
Offline
|
Date Last Visited |
11-11-2020
02:23 AM
|