I've been developing my program with just one dataset and I got all the bugs worked out from that one dataset. Now, I am trying it on some other datasets to make sure everything is working as expected.
I am having a problem in the setting-up-the-environments stage of my program. The user inputs a polygon shapefile that contains a polygon representing the parcel outline. This polygon is used to create a mask for future raster calculations. The basic outline of this part of the program is as follows:
- Buffer the outline polygon by 50 ft
- Convert the buffered polygon into a raster
- Set the raster as the mask and the snapraster for the remainder of the program
An empty raster is created when converting the buffered polygon to a raster. I tried doing it manually in GIS and had the same problem. I added an integer field to my buffered shapefile and converted to a raster based on that integer field (rather than the ID or FID which were both 0.) That didn't fix it. After setting the processing extent to the same extent as my DEM, it worked. I made those changes to my code and it still is not producing the expected results. It's not giving me an error, but it tells me that the raster is empty and then keeps running through my code. I've exhausted my troubleshooting abilities. I also want to repeat: it works when I do it by hand in GIS but it isn't working via code for some reason.
My code is printing out the following:
Preparing necessary files...
Checking for and creating output folders for GIS files...
Checking for and creating output space for tables...
Setting cell size to DEM cell size: 5.0 ft...
Creating an environment mask from the site outline shapefile...
All cells in grid k:\gradwork\gis\colleg~1\cp\scratch\rastermask have NODATA value. VAT will not be built. (VatBldNew)
My code is below. Line 84 is where things are going wrong.
#### _______________IMPORT MODULES_________________###
print("Preparing necessary files...")
import os
import arcpy
import copy
### ______________INPUT FILES___________________###
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 J\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)
### ________________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 GIS files...")
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 output table folder if it does not exist and create project folder
print("Checking for and creating output space for tables...")
TabPath = r"C:\outputtable" "\\"
ProjFolder = TabPath + filename[:-4]
if not os.path.exists(TabPath):
os.makedirs(TabPath)
if not os.path.exists(ProjFolder):
os.makedirs(ProjFolder)
### _____________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
arcpy.env.extent = arcpy.Describe(DEM).extent
# 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",)
arcpy.AddField_management(maskshp, "RID", "SHORT")
with arcpy.da.UpdateCursor(maskshp, ["RID"]) as cursor:
for row in cursor:
row[0] = 8
cursor.updateRow(row)
# Convert buffer to raster
mask = arcpy.Raster(arcpy.PolygonToRaster_conversion(maskshp, "RID", 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 = maskEDIT: Changing line 69 to "arcpy.env.extent = DEM.extent" did not solve the problem.
EDIT 2: My files aren't in the same coordinate system!
