Skip to content

rect_grid

Overlay a square grid to the given areas of a gpd.GeoDataFrame.

Note

This function fits a rectangular grid with user defined resolution and optional overlap. The spatial predicates can be used to filter the grid cells that intersect, or are contained strictly within the given input GeoDataFrame.

Note

Returns None if the gdf is empty.

Parameters:

Name Type Description Default
gdf GeoDataFrame

GeoDataFrame to fit the grid to. Uses the bounding box of the GeoDataFrame to fit the grid.

 required
resolution Tuple[int, int]

Patch size/resolution of the grid (in pixels).

 (256, 256)
overlap int

overlap of the cells in the grid (in percentages).

 0
predicate str

Predicate to use for the spatial join, by default "intersects". Allowed values are "intersects", "within", "contains", "contains_properly".

 'intersects'

Returns:

Type Description
GeoDataFrame

gpd.GeoDataFrame: GeoDataFrame with the grid fitted to the given GeoDataFrame.

Raises:

Type Description
ValueError

If predicate is not one of "intersects" or "within".

Examples:

>>> from histolytics.spatial_ops.rect_grid import rect_grid
>>> from histolytics.data import cervix_tissue
>>>
>>> # get the stromal tissue
>>> tis = cervix_tissue()
>>> stroma = tis[tis["class_name"] == "stroma"]
>>>
>>> # fit a rectangular grid strictly within the stromal tissue
>>> grid = rect_grid(stroma, resolution=(256, 256), overlap=0, predicate="contains")
>>> print(grid.head(3))
                                            geometry
0  POLYGON ((5443 626, 5699 626, 5699 882, 5443 8...
1  POLYGON ((4419 882, 4675 882, 4675 1138, 4419 ...
2  POLYGON ((4675 882, 4931 882, 4931 1138, 4675 ...
>>> ax = tis.plot(column="class_name", figsize=(5, 5), aspect=1, alpha=0.5)
>>> grid.plot(ax=ax, edgecolor="black", facecolor="none", lw=1)
>>> ax.set_axis_off()

out

Source code in src/histolytics/spatial_ops/rect_grid.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def rect_grid(
    gdf: gpd.GeoDataFrame,
    resolution: Tuple[int, int] = (256, 256),
    overlap: int = 0,
    predicate: str = "intersects",
) -> gpd.GeoDataFrame:
    """Overlay a square grid to the given areas of a `gpd.GeoDataFrame`.

    Note:
        This function fits a rectangular grid with user defined resolution and optional
        overlap. The spatial predicates can be used to filter the grid cells that
        intersect, or are contained strictly within the given input GeoDataFrame.

    Note:
        Returns None if the gdf is empty.

    Parameters:
        gdf (gpd.GeoDataFrame):
            GeoDataFrame to fit the grid to. Uses the bounding box of the GeoDataFrame
            to fit the grid.
        resolution (Tuple[int, int]):
            Patch size/resolution of the grid (in pixels).
        overlap (int):
            overlap of the cells in the grid (in percentages).
        predicate (str):
            Predicate to use for the spatial join, by default "intersects".
            Allowed values are "intersects", "within", "contains", "contains_properly".

    Returns:
        gpd.GeoDataFrame:
            GeoDataFrame with the grid fitted to the given GeoDataFrame.

    Raises:
        ValueError: If predicate is not one of "intersects" or "within".

    Examples:
        >>> from histolytics.spatial_ops.rect_grid import rect_grid
        >>> from histolytics.data import cervix_tissue
        >>>
        >>> # get the stromal tissue
        >>> tis = cervix_tissue()
        >>> stroma = tis[tis["class_name"] == "stroma"]
        >>>
        >>> # fit a rectangular grid strictly within the stromal tissue
        >>> grid = rect_grid(stroma, resolution=(256, 256), overlap=0, predicate="contains")
        >>> print(grid.head(3))
                                                    geometry
        0  POLYGON ((5443 626, 5699 626, 5699 882, 5443 8...
        1  POLYGON ((4419 882, 4675 882, 4675 1138, 4419 ...
        2  POLYGON ((4675 882, 4931 882, 4931 1138, 4675 ...
        >>> ax = tis.plot(column="class_name", figsize=(5, 5), aspect=1, alpha=0.5)
        >>> grid.plot(ax=ax, edgecolor="black", facecolor="none", lw=1)
        >>> ax.set_axis_off()
    ![out](../../img/rect_grid.png)
    """
    if gdf.empty or gdf is None:
        return

    allowed = ["intersects", "within", "contains", "contains_properly"]
    if predicate not in allowed:
        raise ValueError(f"predicate must be one of {allowed}. Got {predicate}")

    if not (0 <= overlap < 100):
        raise ValueError("overlap must be in the range [0, 100)")

    stride = (
        int(resolution[0] * (1 - overlap / 100)),
        int(resolution[1] * (1 - overlap / 100)),
    )

    grid = _full_rect_grid(gdf, resolution, stride, pad=20)
    grid = grid.set_crs(gdf.crs, allow_override=True)
    _, grid_inds = grid.sindex.query(gdf.geometry, predicate=predicate)
    grid = grid.iloc[np.unique(grid_inds)]

    return grid.drop_duplicates("geometry").reset_index(drop=True)