Skip to content

clip_polygon crop=True drops boundary pixels when all_touched=True #1197

@brendancol

Description

@brendancol

Describe the bug

clip_polygon with crop=True and all_touched=True silently drops boundary pixels whose cell area overlaps the polygon but whose center coordinate falls outside the polygon's bounding box.

_crop_to_bbox compares pixel center coordinates against the geometry bounding box with >=/<=. In normal mode (all_touched=False) this is fine because only pixels whose centers fall inside the polygon get rasterized. But all_touched=True tells rasterize to include any pixel whose cell intersects the polygon, and cells extend half a pixel past their center. The crop step has no knowledge of all_touched, so it removes those edge pixels before rasterize gets a chance to see them.

Reproducer

import numpy as np
from shapely.geometry import Polygon
from xrspatial.polygon_clip import clip_polygon
from xrspatial.tests.general_checks import create_test_raster

data = np.arange(48, dtype=np.float64).reshape(8, 6)
raster = create_test_raster(data)  # pixel spacing 0.5, x in [0, 2.5]

# Left edge at x=0.15 -- inside the cell of pixel x=0.0 (cell [-0.25, 0.25])
poly = Polygon([(0.15, 0.0), (0.15, 3.5), (2.5, 3.5), (2.5, 0.0)])

result_crop = clip_polygon(raster, poly, crop=True, all_touched=True)
result_nocrop = clip_polygon(raster, poly, crop=False, all_touched=True)

n_crop = np.count_nonzero(np.isfinite(result_crop.values))
n_nocrop = np.count_nonzero(np.isfinite(result_nocrop.values))
print(f"crop=True: {n_crop}, crop=False: {n_nocrop}")
# crop=True: 40, crop=False: 48 -- 8 pixels silently lost

Expected behavior

crop=True, all_touched=True should include exactly the same pixels as crop=False, all_touched=True, just on a smaller output array.

Root cause

_crop_to_bbox doesn't receive the all_touched flag. When all_touched is set, the bbox comparison needs to expand by half a pixel on each side so that pixels whose cells overlap the geometry boundary aren't thrown out during the crop.

Proposed fix

  1. Pass all_touched into _crop_to_bbox.
  2. Compute pixel spacing from the coordinate arrays.
  3. When all_touched=True, pad minx/miny/maxx/maxy by half a pixel before comparing against pixel centers.
  4. Add a test confirming crop=True, all_touched=True matches crop=False, all_touched=True for a polygon whose edges land between pixel centers.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions