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:
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!
Solved! Go to Solution.
Is it possible that your outline shapefile is in a different projection than your DEM?
Is it possible that your outline shapefile is in a different projection than your DEM?
: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?
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
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?
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.
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!
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.
All I had to do was set an output coordinate system and it worked! Yay!
Of course, now I'm getting a different error but hopefully it's unrelated to this one.