How to check for overlap: Raster and polygon

[BennoPerchermeier]

I have two input files:

- a feature class with a number of polygons 

- a raster

I want to do an Extract by Mask for each of the polygons in my feature class. Due to my data structure, I sometimes do not have an overlap between one of the polygons and the raster. Polygon 3 in the picture below does not have an overlap with the raster. 

I already have a procedure that gives me the desired results, but I still have to do ExtractByMask for all of the polygons. All I do is not saving the temporary raster if the result is empty. Is there a way or a function that can check if a polygon and a raster intersect before I do the extract by mask procedure?

Below is my - probably way too complicated - code. It does the job, but it is very slow.

(The way I check if the raster is empty is by checking if the mean value and the standard Deviation is equal to 0)

import arcpy

arcpy.env.workspace = "my_workspace"
raster = "my_raster"
fc_polygons = "my_featureClass"

# loop over polygons in the feature class fc_polygons
flds = ["OBJECTID", "SHAPE@"]
with arcpy.da.SearchCursor(fc_polygons, flds) as curs:
    for row in curs:
        objectId = row[0]
        arcpy.AddMessage("Extracting by mask for polygon: {0}".format(objectId))
        polygon = row[flds.index("SHAPE@")]
        # extracting for the polygon
        new_raster = arcpy.sa.ExtractByMask(raster, polygon)
        # if new_raster is empty discard the temporary raster, otherwise save it
        if new_raster.mean == 0 and new_raster.standardDeviation == 0:
            arcpy.Delete_management(new_raster)
            arcpy.AddMessage("resulting raster is empty and is, therefore, discarded.")
        else:
            outname = "ras_buf_{0}".format(objectId)
            new_raster.save(outname)
            arcpy.AddMessage(" - saved raster as {0}".format(outname))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Edit: Curtis suggested using the extent of the raster, but my raster is not continuous. Below is another picture of how I imagine an intermediary step to look like. Only the polygons that overlap with the raster are selected:

[curtvprice]

I suggest using the extent of the raster as a polygon to select those that overlap before you start processing the polys. Then use the layer (with just the polygons that overlap) for your cursor loop.

lyrPoly = arcpy.MakeFeatureLayer_management("my_featureClass", "lyrPoly")
# make a polygon from extent
ext = arcpy.Describe("my_raster").extent
array = arcpy.Array()
array.add(arcpy.Point(ext.XMin, ext.YMin))
array.add(arcpy.Point(ext.XMin, ext.YMax))
array.add(arcpy.Point(ext.XMax, ext.YMax))
array.add(arcpy.Point(ext.XMas, ext.YMin))
array.add(arcpy.Point(ext.XMin, ext.YMin))
extPoly = arcpy.Polygon(array)
# Select polygons that intersect the extent of the raster
arcpy.SelectLayerByLocation(lyrPoly, "INTERSECT", extPoly)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

http://desktop.arcgis.com/en/arcmap/latest/analyze/python/using-geometry-objects-with-geoprocessing-...

[BennoPerchermeier]

Thank you for your answer Curtis! I like your idea because it seems much faster than what I am doing. Though I did not explain the issue properly. My raster is actually not continuous and there are empty polygons in the middle.

Running the script above selects all of the polygons:

I will update my question accordingly. Do you have another idea, how to do this?

[curtvprice · Accepted solution]

I would run the aggregate tool on your raster with maximum to make a coarse (ie small size) raster, convert this to a polygon. Select By Location against with polygon to select the polygons that overlap with your raster data cells. Note this workflow does require a Spatial Analyst license.

(UPDATE fixed a typo in the code the guaranteed writable temp folder path is env.scratchFolder)

from arcpy import env
from arcpy.sa import *

lyrPoly = arcpy.MakeFeatureLayer_management("my_featureClass", "lyrPoly")

# create an generalized extent polygon of your raster
tmpAgg = os.path.join(env.scratchFolder, "xagg.tif")
tmpPoly = os.path.join(env.scratchFolder, "xaggpoly.shp")
agg = SetNull(IsNull(Aggregate("my_raster", 5, "MAXIMUM", "EXPAND")), 1)
agg.save(tmpAgg)
arcpy.RasterToPolygon_conversion(tmpAgg, tmpPoly)
# Select polygons that intersect the extent of the raster data areas
arcpy.SelectLayerByLocation(lyrPoly, "INTERSECT", tmpPoly)
arcpy.Delete_management(tmpAgg)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍arcpy.Delete_management(tmpPoly)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[BennoPerchermeier]

Awesome thanks! This idea works. The aggregation is relatively slow but still faster than trying to do an Extract by mask on all the empty polygons.

I had to amend your script to make it work, so here it is:

import os
import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")

raster = "my_raster"
fc_polygons = "my_featureClass"
tempFolder = r'C:\Path\to\temporary\folder'

# create temporary layer
lyrPoly = arcpy.MakeFeatureLayer_management(fc_polygons, "lyrPoly")

# create an generalized extent polygon of your raster
tmpAgg = os.path.join(tempFolder, "xagg2.tif")
tmpPoly = os.path.join(tempFolder, "xaggpoly.shp")
# create aggregated raster with larger cell sizes
agg = SetNull(IsNull(Aggregate(raster, 5, "MAXIMUM", "EXPAND")), 1)
agg.save(tmpAgg)
# convert to polygon feature class 
arcpy.RasterToPolygon_conversion(tmpAgg, tmpPoly)
# Select polygons that intersect the extent of the raster data areas
arcpy.SelectLayerByLocation_management(lyrPoly, "INTERSECT", tmpPoly)
arcpy.Delete_management(tmpAgg)
arcpy.Delete_management(tmpPoly)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I had to specify a temporary folder where I could save the temporary raster and the temporary polygon. For one, because I could not use the command "env.scratchRaster". For the other, because I could not save the tiff file to the workspace, because can't be saved to GDBs as far as I know.

Anyways, thank you very much for your answer, Curtis!!

[curtvprice]

This was a typo on my part, env.scratchFolder is what I meant. Sorry about that. I will edit the code above since it is marked as the answer. So glad that worked for you!

If performance is a problem (ie you do this a lot) you could probably speed it up two ways: 1) Up the aggregate from 5 to 10 or 20 and 2) use the "NO_SIMPLIFY" option to the Raster To Polygon step to avoid the generalizing which for a spatially complicated raster may slow things down a bit. T hough there is no way to speed it up too much because you have to look at every cell to capture all the raster geometry.

I just thought of a second approach, which would be to run Zonal Statistics As Table of your raster with the test polygons as zones. This would create a table with one row for every polygon that had an overlap. You could then perform an attribute join from your polygons to the Zonal Statistics As Table output to select the polygons that overlap. This may be much faster as there would be no raster to polygon conversion (only polygon to raster, which should be faster). But of course if your problem is solved, the discussion is kind of academic... I'm into that these days.

[BennoPerchermeier]

Thanks for the speeding up advice

If performance is a problem (ie you do this a lot) you could probably speed it up two ways: 1) Up the aggregate from 5 to 10 or 20 and 2) use the "NO_SIMPLIFY" option to the Raster To Polygon step to avoid the generalizing which for a spatially complicated raster may slow things down a bit.

I accelerated it by using 10 instead of 5 in the Aggregate function and use the "NO_SIMPLIFY" option. I also didn't have to save the temporary aggregated raster "agg". This saves time, too.

My final script:

import os
from arcpy import env
from arcpy.sa import *

lyrPoly = arcpy.MakeFeatureLayer_management("my_featureClass", "lyrPoly")

# create an generalized extent polygon of your raster
tmpPoly = os.path.join(tempFolder, "xaggpoly.shp")
agg = SetNull(IsNull(Aggregate("my_raster", 10, "MAXIMUM", "EXPAND")), 1)
arcpy.RasterToPolygon_conversion(agg, tmpPoly, "NO_SIMPLIFY")

# Select polygons that intersect the extent of the raster data areas
arcpy.SelectLayerByLocation_management(lyrPoly, "INTERSECT", tmpPoly)
arcpy.Delete_management(tmpPoly)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

PS: If that doesn't work out I will come back to our academic discussion about using zonal statistics