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)
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