Troubleshooting a Polygon to Raster conversion

[RachaelJohnson]

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:

1. Buffer the outline polygon by 50 ft
2. Convert the buffered polygon into a raster
3. 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 = mask

EDIT:  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!

[SepheFox · Accepted solution]

Is it possible that your outline shapefile is in a different projection than your DEM?

[RachaelJohnson]

:U  OH NO!!

My DEM is a projected coordinate system (NAD 1983 HARN VA State Plane South) but my buffer is in a GCS, NAD 1983 HARN.  How do I resolve this discrepancy in Python?

[SepheFox]

Well, want your data to be in projected coordinates rather than geographic, so I would start by reprojecting the original outline shapefile, before buffering, using the Project tool.

# project data

  arcpy.Project_management(infc, outfc, outCS)

Help file here: ArcGIS Desktop

[RachaelJohnson]

Would it be just as good to set my output CS to the projected system my DEM is in and have GIS take care of the conversion automatically?

[SepheFox]

Do you mean the output from buffering? Honestly, you're playing with fire anytime you start with data that is not all in the same projected coordinate system. For instance, your buffers are being created from shapes that are in a different CS, so there could be unexpected results. It's best practice to get everything into the same PCS before you start any geoprocessing, and the Project tool is the only sure way to do that correctly.

[RachaelJohnson]

I meant setting my environments so the output CS is the same as my DEM CS, so that the output of any tool results in PCS.  But, it probably is better to write a code block that checks projections before the code actually starts.  Thanks for the tip! 

[SepheFox]

Yeah, I really would recommend checking and the Projecting anything that is different as the first step. Sometimes there can be subtle effects you may not realize.

[RachaelJohnson]

All I had to do was set an output coordinate system and it worked! Yay!

• arcpy.env.outputCoordinateSystem = arcpy.Describe(DEM).spatialReference

Of course, now I'm getting a different error but hopefully it's unrelated to this one.