Shapely | Binary mask of Polygons’ intersection

multi-polygonnumpypolygonpython 3shapely

Goal: Create a binary 2d-array, that represents Points of intersection of a square Polygon with a MultiPolygon.

from typing import List, Tuple
import numpy.typing as npt
from shapely.geometry import Polygon, Point
from shapely.geometry.multipolygon import MultiPolygon

def binary_mask(tile_poly: Polygon, wkt_mp: MultiPolygon, tile_shape: Tuple[int, int] = (256, 256)) -> npt.NDArray[np.int_]:
    mask = np.zeros(tile_shape)
    for x in range(tile_shape[0]):
        for y in range(tile_shape[1]):
            if tile_poly.contains(Point(x,y)) or wkt_mp.contains(Point(x,y)):
                mask[x, y] = 1

    return mask

However, all arrays only return 0.

mask = binary_mask(tile_poly, wkt_mp, tile_shape[:2])
print(np.unique(mask))
>>> [0.]

Most tile_poly objects lay within the wkt_poly; so should return an array of only 1s.

Any tile_poly that is on the wkt_poly boundary should contain some 0s; where its Points do not intersect, i.e. go outside wkt_poly's boundary.

Data

crs='EPSG:4326'

qupath_poly_flip_y[0]
>>> [<shapely.geometry.polygon.Polygon at 0x7fe4fe9f6af0>,
     <shapely.geometry.polygon.Polygon at 0x7fe4fe9f6bb0>,
     ...]

qupath_poly_flip_y[0][0].area  # size(256, 256)
>>> 65536.0

qupath_poly_flip_y[0][0].bounds
>>> (19200.0, 11336.0, 19456.0, 11592.0)

wkt_multipoly[0].area
>>> 17431020.0

wkt_multipoly[0].bounds
>>> (8474.0, 10026.0, 21732.0, 17267.0)

Best Answer

Instead of referencing a static tile_shape, for x and y coords, I needed the actual coords from the Polygon itself, using .bounds attribute.

For use in range(), they had to be cast as int().

def binary_mask(tile_poly: Polygon, wkt_mp: MultiPolygon) -> npt.NDArray[np.int_]:
    min_x, min_y, max_x, max_y = tuple(map(int, tile_poly.bounds))
    mask = np.zeros((max_x - min_x, max_y - min_y))
    
    for x in range(min_x, max_x):
        for y in range(min_y, max_y):
            if wkt_mp.contains(Point(x,y)):
                mask[x - min_x, y - min_y] = 1

    return mask
Related Question