diff --git a/docs/getting_started.rst b/docs/getting_started.rst index 1f5869f..55e171d 100644 --- a/docs/getting_started.rst +++ b/docs/getting_started.rst @@ -259,8 +259,11 @@ What's next Pandas DataFrame, or Earth Engine Feature Collection (but they *must* have point geometries). Or use :func:`~geebeam.grid_and_run_pipeline` to sample on a regular grid rather than randomly for when you need complete spatial - coverage without gaps. See :mod:`~geebeam.sampler` - for more info, including custom dataset splits (e.g. train/val/test). + coverage without gaps. By default a grid tile is kept only if its sample point falls + inside the region; pass ``tile_coverage="intersect"`` to instead keep every + tile whose footprint touches the region, for gap-free coverage right up to the + edges. See :mod:`~geebeam.sampler` for more info, including custom + dataset splits (e.g. train/val/test). - **Split processing of large patches** — Earth Engine limits single request sizes to 50 MB, but geebeam lets you split processing: :doc:`split_processing` - **API reference** — full parameter documentation for all three pipeline diff --git a/src/geebeam/pipeline.py b/src/geebeam/pipeline.py index e8fb41e..536aaf3 100644 --- a/src/geebeam/pipeline.py +++ b/src/geebeam/pipeline.py @@ -70,7 +70,6 @@ def _prepare_run_metadata(config, align_transform=None): return scale_x, scale_y - def _apply_position_offset(input_records, position, patch_size, scale_x, scale_y, align_transform=None): """Add x_topleft/y_topleft to each record; original x/y is preserved. @@ -79,23 +78,7 @@ def _apply_position_offset(input_records, position, patch_size, scale_x, scale_y target crs) is provided, each patch's top-left corner is snapped onto that transform's pixel grid so the extracted patch aligns exactly. """ - _valid_positions = {'center', 'top-left', 'top-right', 'bottom-left', 'bottom-right'} - if position == 'top-left': - dx, dy = 0, 0 - elif position == 'center': - dx = -(patch_size / 2) * scale_x - dy = -(patch_size / 2) * scale_y - elif position == 'top-right': - dx = -patch_size * scale_x - dy = 0 - elif position == 'bottom-left': - dx = 0 - dy = -patch_size * scale_y - elif position == 'bottom-right': - dx = -patch_size * scale_x - dy = -patch_size * scale_y - else: - raise ValueError(f"Invalid position '{position}'. Must be one of: {sorted(_valid_positions)}") + dx, dy = sampler._position_offset(position, patch_size, scale_x, scale_y) records = [{**r, 'x_topleft': r['x'] + dx, 'y_topleft': r['y'] + dy} for r in input_records] if align_transform is not None: x0, y0, sx, sy = sampler._parse_transform(align_transform) @@ -409,6 +392,8 @@ def grid_and_run_pipeline( buffer_distance: float = 0, validation_ratio: float = 0, random_seed: int = 0, + position: str = 'center', + tile_coverage: str = 'clip', **kwargs ) -> None: """Sample points from regular grid and then run a Beam pipeline to download image chips from Earth Engine. @@ -433,6 +418,16 @@ def grid_and_run_pipeline( Can be used to ensure complete coverage at edges of sampling_region. validation_ratio: Fraction of points to mark as validation (can be 0.0). random_seed: Seed for random sampling + position: Where each sampling point falls within its patch. One of 'center' + (default), 'top-left', 'top-right', 'bottom-left', 'bottom-right'. Used + for tile_coverage check. + tile_coverage: Which grid tiles to keep relative to ``sampling_region``. + 'clip' (default): keep a tile only if its sampling point (the patch reference + point set by ``position``) falls inside the region. + 'center_clip': keep a tile only if its patch center falls inside the region, + regardless of ``position``. + 'intersect': keep any tile whose full ``patch_size`` footprint touches the region, + guaranteeing gap-free coverage at the edges. **kwargs: Additional keyword arguments are documented in :meth:`pipeline.run_pipeline`. """ @@ -446,6 +441,9 @@ def grid_and_run_pipeline( scale=scale, buffer_distance=buffer_distance, align_transform=align_transform, + patch_size=patch_size, + position=position, + tile_coverage=tile_coverage, ) if validation_ratio > 0: @@ -463,4 +461,5 @@ def grid_and_run_pipeline( patch_size=patch_size, scale=scale, align_transform=align_transform, + position=position, **kwargs) diff --git a/src/geebeam/sampler.py b/src/geebeam/sampler.py index d58fe9a..516fef1 100644 --- a/src/geebeam/sampler.py +++ b/src/geebeam/sampler.py @@ -44,6 +44,45 @@ def _snap_to_grid(x, y, x0, y0, scale_x, scale_y): return float(xs), float(ys) return xs, ys +def _position_offset(position, patch_size, scale_x, scale_y): + """Return (dx, dy) from a sampling point to its patch top-left corner. + + ``position`` names where the sampling point sits within the patch; adding the + returned offset to the point gives the patch's top-left corner (the translateX/ + translateY used to extract pixels). + """ + _valid_positions = {'center', 'top-left', 'top-right', 'bottom-left', 'bottom-right'} + if position == 'top-left': + dx, dy = 0, 0 + elif position == 'center': + dx = -(patch_size / 2) * scale_x + dy = -(patch_size / 2) * scale_y + elif position == 'top-right': + dx = -patch_size * scale_x + dy = 0 + elif position == 'bottom-left': + dx = 0 + dy = -patch_size * scale_y + elif position == 'bottom-right': + dx = -patch_size * scale_x + dy = -patch_size * scale_y + else: + raise ValueError(f"Invalid position '{position}'. Must be one of: {sorted(_valid_positions)}") + return dx, dy + +def _pad_locs(locs, step, n): + """Extend a 1-D array of grid node coords by ``n`` extra steps on each side. + + The added nodes lie on the same lattice as ``locs`` (same ``step``), so the original + nodes keep their exact float values. Used to widen the candidate grid for the + 'center_clip'/'intersect' selectors without shifting the interior grid. + """ + if n <= 0: + return locs + left = locs[0] - np.arange(n, 0, -1)*step + right = locs[-1] + np.arange(1, n+1)*step + return np.concatenate([left, locs, right]) + def _get_roi( sampling_region: str | ee.Geometry | gpd.GeoDataFrame, target_crs: str @@ -123,12 +162,23 @@ def sample_region_random( buffer_distance: float = 0, align_transform: Affine | tuple[float] | list[float] | None = None, ) -> gpd.GeoDataFrame: - """Get random points within region of interest. + """Sample random points within a region of interest. - If align_transform (a rasterio Affine or list/tuple length 6 in the target crs) is provided, - the sampled points are snapped onto that transform's pixel grid. Note that snapping to a grid - that is coarse relative to the sampling density can collapse distinct random points - onto the same location; a warning is emitted if this happens. + Args: + roi: Region to sample from. May be a path to a vector file (str), an + ``ee.Geometry``, a ``shapely`` geometry, or a ``gpd.GeoDataFrame``; + a differing CRS is reprojected to ``crs``. + crs: Target CRS for the returned points (e.g. 'EPSG:4326'). + n_sample: Number of random points to sample/draw. + random_seed: Seed for reproducible sampling. Defaults to 0. + buffer_distance: Distance in meters to buffer ``roi`` before sampling, so + points can fall slightly outside the region. Defaults to 0 (no buffer). + align_transform: Optional rasterio ``Affine`` or length-6 tuple/list (in + ``crs`` units). If given, each point is snapped onto that transform's + pixel grid; a warning is emitted if snapping produces duplicate locations. + + Returns: + GeoDataFrame of point geometries in ``crs``. """ rng = np.random.default_rng(random_seed) roi = _get_roi(roi, crs) @@ -139,9 +189,9 @@ def sample_region_random( sampled_points = gpd.GeoDataFrame(geometry=roi.sample_points(n_sample, rng=rng).geometry.explode(), crs=crs) if align_transform is not None: - x0, y0, sx, sy = _parse_transform(align_transform) + x0, y0, scale_x, scale_y = _parse_transform(align_transform) xs, ys = _snap_to_grid(sampled_points.geometry.x.values, - sampled_points.geometry.y.values, x0, y0, sx, sy) + sampled_points.geometry.y.values, x0, y0, scale_x, scale_y) n_unique = np.unique(np.column_stack([xs, ys]), axis=0).shape[0] n_collisions = len(xs) - n_unique if n_collisions > 0: @@ -164,27 +214,65 @@ def sample_region_grid( scale: float | None = None, align_transform: Affine | tuple[float] | list[float] | None = None, buffer_distance: float = 0, + patch_size: int | None = None, + position: str = 'center', + tile_coverage: str = 'clip', ) -> gpd.GeoDataFrame: - """Get a regular grid of points covering region of interest. + """Build a regular grid of sampling points covering a region of interest. - ``scale`` (grid spacing in meters, as scale*stride) is required unless ``align_transform`` - is provided. If transform (a rasterio Affine or list/tuple in the target crs) is - provided, the grid uses that transform's pixel size (``scale`` is ignored) and its origin - is anchored to the transform, so every location falls on that transform's pixel grid. + Args: + roi: Region to sample from. May be a path to a vector file (str), an + ``ee.Geometry``, a ``shapely`` geometry, or a ``gpd.GeoDataFrame``; + a differing CRS is reprojected to ``crs``. + crs: Target CRS for the returned points (e.g. 'EPSG:4326'). + stride: Spacing between adjacent points, in pixels (the spacing in ``crs`` + units is ``stride`` times the pixel size). + scale: Pixel size / export resolution in meters. Required unless + ``align_transform`` is given. + align_transform: Optional rasterio ``Affine`` or length-6 tuple/list (in + ``crs`` units). If given, the grid uses this transform's pixel size + (``scale`` is ignored) and is anchored to its origin, so every point + lands on the transform's pixel grid. Defaults to None. + buffer_distance: Distance in meters to buffer ``roi`` before gridding, to + help ensure coverage at the edges. Defaults to 0 (no buffer). + patch_size: Patch size in pixels. Only required when ``tile_coverage`` is not + 'clip'. Defaults to None. + position: Where each point sits within its patch: 'center' (default), + 'top-left', 'top-right', 'bottom-left', or 'bottom-right'. Only + matters when ``tile_coverage`` is not 'clip'. + tile_coverage: Which tiles to keep relative to ``roi``. 'clip' (default) + keeps a tile if its sampling point falls inside ``roi``; 'center_clip' + keeps a tile if its patch center falls inside ``roi``, regardless of + ``position``; 'intersect' keeps a tile if any part of its patch footprint + touches ``roi``. The latter two require ``patch_size``. + + Returns: + GeoDataFrame of grid point geometries in ``crs``. + + Raises: + ValueError: If ``tile_coverage`` is invalid, if 'center_clip'/'intersect' is + requested without ``patch_size``, or if neither ``scale`` nor + ``align_transform`` is provided. """ + if tile_coverage not in ('clip', 'center_clip', 'intersect'): + raise ValueError(f"Invalid tile_coverage '{tile_coverage}'. " + "Must be 'clip', 'center_clip', or 'intersect'.") + if tile_coverage in ('center_clip', 'intersect') and patch_size is None: + raise ValueError(f"tile_coverage='{tile_coverage}' requires patch_size.") if align_transform is None and scale is None: raise ValueError('`scale` is required unless align_transform is provided.') roi = _get_roi(roi, crs) if align_transform is not None: - x0, y0, sx, sy = _parse_transform(align_transform) + x0, y0, scale_x, scale_y = _parse_transform(align_transform) if buffer_distance != 0: scale_proj_1m = _get_crs_scale(roi.crs.to_string(), 1) roi = roi.dissolve().buffer(scale_proj_1m*buffer_distance) xmin, ymin, xmax, ymax = roi.total_bounds - step_x, step_y = sx*stride, sy*stride - # Anchor the grid to the reference transform - x_start = x0 + np.floor((xmin - x0)/sx)*sx - y_start = y0 + np.floor((ymin - y0)/sy)*sy + step_x, step_y = scale_x*stride, scale_y*stride + + # Align the grid to the reference transform + x_start = x0 + np.floor((xmin - x0)/scale_x)*scale_x + y_start = y0 + np.floor((ymin - y0)/scale_y)*scale_y x_locs = np.arange(x_start, xmax+step_x, step_x) y_locs = np.arange(y_start, ymax+step_y, step_y) else: @@ -192,16 +280,39 @@ def sample_region_grid( if buffer_distance != 0: scale_proj_1m = scale_proj/scale roi = roi.dissolve().buffer(scale_proj_1m*buffer_distance) + scale_x = scale_y = scale_proj xmin, ymin, xmax, ymax = roi.total_bounds - x_locs = np.arange(xmin, xmax+scale_proj*stride, scale_proj*stride) - y_locs = np.arange(ymin, ymax+scale_proj*stride, scale_proj*stride) + step_x = step_y = scale_proj*stride + x_locs = np.arange(xmin, xmax+step_x, step_x) + y_locs = np.arange(ymin, ymax+step_y, step_y) + + if tile_coverage in ('center_clip', 'intersect'): + # Pad the candidate grid by a full patch on all sides + n_pad = int(np.ceil(patch_size/stride)) + x_locs = _pad_locs(x_locs, step_x, n_pad) + y_locs = _pad_locs(y_locs, step_y, n_pad) + meshgrid = np.array(np.meshgrid(x_locs, y_locs)).T.reshape(-1, 2) x_all, y_all = meshgrid[:,0], meshgrid[:,1] points_gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x_all, y_all), crs=crs) - # Clip to region - points_gdf = gpd.clip(points_gdf, roi) + if tile_coverage == 'clip': + # Keep points whose location (the patch reference point) falls inside the region. + points_gdf = gpd.clip(points_gdf, roi) + else: + roi_geom = roi.union_all() if hasattr(roi, 'union_all') else roi.unary_union + dx, dy = _position_offset(position, patch_size, scale_x, scale_y) + if tile_coverage == 'center_clip': + # Keep points whose patch center falls inside the region (position-independent). + cx, cy = dx + (patch_size/2)*scale_x, dy + (patch_size/2)*scale_y + selectors = gpd.points_from_xy(x_all + cx, y_all + cy) + else: # 'intersect' + # Keep points whose full patch footprint touches the region. + w, h = patch_size * scale_x, patch_size * scale_y + selectors = shapely.box(x_all + dx, y_all + dy, x_all + dx + w, y_all + dy + h) + mask = gpd.GeoSeries(selectors, crs=crs).intersects(roi_geom) + points_gdf = points_gdf[np.asarray(mask)] points_gdf.index = np.arange(points_gdf.shape[0]) return points_gdf diff --git a/tests/test_sampler.py b/tests/test_sampler.py index baff7b0..2c54e03 100644 --- a/tests/test_sampler.py +++ b/tests/test_sampler.py @@ -10,6 +10,7 @@ from geebeam.sampler import ( _get_roi, _parse_transform, + _position_offset, _process_sampling_points, _snap_to_grid, sample_region_grid, @@ -217,6 +218,104 @@ def test_sample_region_grid_requires_scale(): with pytest.raises(ValueError, match='scale'): sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1) +def test_sample_region_grid_intersect_more_points(): + """tile_coverage='intersect' keeps edge tiles the default clip rule drops.""" + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with patch('geebeam.sampler._get_crs_scale', return_value=0.1): + clip = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, + scale=1000.0, patch_size=2) + intersect = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, + scale=1000.0, patch_size=2, + tile_coverage='intersect') + # Edge tiles (center outside ROI but footprint overlapping) are retained. + assert len(intersect) > len(clip) + # Every retained tile's footprint really does touch the ROI. + roi_geom = roi_gdf.union_all() + half = 2 * 0.1 / 2 # patch_size * scale_proj / 2, centered position + for pt in intersect.geometry: + tile = box(pt.x - half, pt.y - half, pt.x + half, pt.y + half) + assert tile.intersects(roi_geom) + +def test_sample_region_grid_intersect_requires_patch_size(): + """intersect mode needs patch_size to build the tile footprint.""" + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with pytest.raises(ValueError, match='patch_size'): + sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + tile_coverage='intersect') + +def test_sample_region_grid_invalid_tile_coverage(): + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with pytest.raises(ValueError, match='tile_coverage'): + sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, tile_coverage='nonsense') + +def test_sample_region_grid_intersect_position(): + """The tile footprint is placed according to `position`.""" + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with patch('geebeam.sampler._get_crs_scale', return_value=0.1): + center = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, tile_coverage='intersect', position='center') + top_left = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, tile_coverage='intersect', position='top-left') + # For 'top-left' the patch extends up/right of the anchor, so retained anchors + # sit further toward the lower-left than for 'center'. + assert top_left.geometry.x.mean() < center.geometry.x.mean() + assert top_left.geometry.y.mean() < center.geometry.y.mean() + +def test_sample_region_grid_center_clip_matches_clip_for_center_position(): + """With position='center' the sampling point is the patch center, so + 'center_clip' and the default 'clip' select the same tiles.""" + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with patch('geebeam.sampler._get_crs_scale', return_value=0.1): + clip = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2) + center_clip = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, tile_coverage='center_clip', + position='center') + clip_pts = {(round(p.x, 6), round(p.y, 6)) for p in clip.geometry} + cc_pts = {(round(p.x, 6), round(p.y, 6)) for p in center_clip.geometry} + assert clip_pts == cc_pts + +def test_sample_region_grid_clip_ignores_position(): + """The default 'clip' rule is position-independent (it clips the sampling point).""" + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with patch('geebeam.sampler._get_crs_scale', return_value=0.1): + clip_center = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, position='center') + clip_top_left = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, position='top-left') + assert clip_center.geometry.x.mean() == clip_top_left.geometry.x.mean() + assert clip_center.geometry.y.mean() == clip_top_left.geometry.y.mean() + +def test_sample_region_grid_center_clip_position_aware(): + """Unlike 'clip', 'center_clip' shifts the selection with `position`.""" + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with patch('geebeam.sampler._get_crs_scale', return_value=0.1): + center = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, tile_coverage='center_clip', position='center') + top_left = sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + patch_size=2, tile_coverage='center_clip', position='top-left') + # 'top-left' places the patch up/right of the anchor, so anchors whose *center* lands + # in the ROI sit further toward the lower-left than for 'center'. + assert top_left.geometry.x.mean() < center.geometry.x.mean() + assert top_left.geometry.y.mean() < center.geometry.y.mean() + +def test_sample_region_grid_center_clip_requires_patch_size(): + roi_gdf = gpd.GeoDataFrame(geometry=[box(0.0, 0.0, 1.0, 1.0)], crs='EPSG:4326') + with pytest.raises(ValueError, match='patch_size'): + sample_region_grid(roi=roi_gdf, crs='EPSG:4326', stride=1, scale=1000.0, + tile_coverage='center_clip') + +def test_position_offset(): + # center: anchor is half a patch in from the top-left corner + assert _position_offset('center', 4, 0.5, 0.5) == (-1.0, -1.0) + # top-left: anchor is the corner, no offset + assert _position_offset('top-left', 4, 0.5, 0.5) == (0, 0) + # bottom-right: corner is a full patch back in both directions + assert _position_offset('bottom-right', 4, 0.5, 0.5) == (-2.0, -2.0) + with pytest.raises(ValueError, match='Invalid position'): + _position_offset('middle', 4, 0.5, 0.5) + def test_sample_region_random_transform_snaps_points(): mock_roi = MagicMock(spec=gpd.GeoDataFrame) mock_points = gpd.points_from_xy([12.3, 27.6], [51.1, 63.9], crs='EPSG:4326')