Find Indices Of Raster Cells That Intersect With A Polygon

I want to get a list of indices (row,col) for all raster cells that fall within or are intersected by a polygon feature. Looking for a solution in python, ideally with gdal/ogr mod

Solution 1:

Since you don't provide a working example, it's bit unclear what your starting point is. I made a dataset with 1 polygon, if you have a dataset with multiple but only want to target a specific polygon you can add SQLStatement or where to the gdal.Rasterize call.

Sample polygon

geojson = """{"type":"FeatureCollection",
"name":"test",
"crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:OGC:1.3:CRS84"}},
"features":[
{"type":"Feature","properties":{},"geometry":{"type":"MultiPolygon","coordinates":[[[[-110.254,44.915],[-114.176,37.644],[-105.729,36.41],[-105.05,43.318],[-110.254,44.915]]]]}}
]}"""

Rasterizing

Rasterizing can be done with gdal.Rasterize. You need to specify the properties of the target grid. If there is no predefined grid these could be extracted from the polygon itself

ds = gdal.Rasterize('/vsimem/tmpfile', geojson, xRes=1, yRes=-1, allTouched=True,
                    outputBounds=[-120, 30, -100, 50], burnValues=1, 
                    outputType=gdal.GDT_Byte)
mask = ds.ReadAsArray()
ds = None
gdal.Unlink('/vsimem/tmpfile')

Converting to indices

Retrieving the indices from the rasterized polygon can be done with Numpy:

y_ind, x_ind = np.where(mask==1)

Solution 2:

Clearly Rutger's solution above is the way to go with this, however I will leave my solution up. I developed a script that accomplished what I needed with the following:

1. Get the bounding box for each vector feature I want to check
2. Use the bounding box to limit the computational window (determine what portion of the raster could potentially have intersections)
3. Iterate over the cells within this part of the raster and construct a polygon geometry for each cell
4. Use ogr.Geometry.Intersects() to check if the cell intersects with the polygon feature

Note that I have only defined the methods, but I think implementation should be pretty clear -- just call match_cells with the appropriate arguments (ogr.Geometry object and geotransform matrix). Code below:

from osgeo import ogr

# Convert projected coordinates to raster cell indices
def parse_coords(x,y,gt):
    row,col =None,None
    if x:
        col =int((x - gt[0]) // gt[1])
        # If only x coordinate is provided, returncolumn index
        if not y:
            return col
    if y:
        row=int((y - gt[3]) // gt[5])
        # If only x coordinate is provided, returncolumn index
        if not x:
            returnrowreturn (row,col)

# Construct polygon geometry from raster cell
def build_cell((row,col),gt):
    xres,yres = gt[1],gt[5]
    x_0,y_0 = gt[0],gt[3]
    top = (yres*row) + y_0
    bottom = (yres*(row+1)) + y_0
    right= (xres*col) + x_0
    left= (xres*(col+1)) + x_0
    # Create ring topology
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(left,bottom)
    ring.AddPoint(right,bottom)
    ring.AddPoint(right,top)
    ring.AddPoint(left,top)
    ring.AddPoint(left,bottom)
    # Create polygon
    box = ogr.Geometry(ogr.wkbPolygon)
    box.AddGeometry(ring)
    return box

# Iterate over feature geometries &checkforintersection
def match_cells(inputGeometry,gt):
    matched_cells = []
    for f,feature in enumerate(inputGeometry):
        geom = feature.GetGeometryRef()
        bbox = geom.GetEnvelope()
        xmin,xmax = [parse_coords(x,None,gt) for x in bbox[:2]]
        ymin,ymax = [parse_coords(None,y,gt) for y in bbox[2:]]
        for cell_row inrange(ymax,ymin+1):
            for cell_col inrange(xmin,xmax+1):
                cell_box = build_cell((cell_row,cell_col),gt)
                if cell_box.Intersects(geom):
                    matched_cells += [[(cell_row,cell_col)]]
    return matched_cells 

Solution 3:

if you want to do this manually you'll need to test each cell for: Square v Polygon intersection and Square v Line intersection.

If you treat each square as a 2d point this becomes easier - it's now a Point v Polygon problem. Check in Game Dev forums for collision algorithms.

Good luck!