Bounds filter fix

[lbferreira]

I also observed the same issue described by @giswqs in #118 related to the output raster with some zero-value on the image edge, resulting in some black strips, If the image footprint image.geometry() is used as the geometry.

I reproduced @giswqs issue and also tested using different image source and geometry. In all cases, the output raster has some weird patterns. Below I present some images showing the problem. After investigating, I realized that this problem was caused because, based on the geometry provided, the current code gets the bounds in EPSG:4326 and then convert the bounds to the target CRS provided by the user. It's problematic because it can result in some deformations. Thus, I change the logic to get the bounds directly in the target CRS provided by the user. In summary, instead of doing geometry.bounds() and later converting the bounds, I do geometry.bounds(maxError=1, proj=target_crs). I'm not sure if 1 is a good value for maxError but it worked well in my tests. If someone, know a better value for maxError, please let me know. At the end I present some images to show the raster after this fix.

Current problems

• Landsat image (as presented by @giswqs)

ee_image = ee.Image("LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318")
ic = ee.ImageCollection(ee_image)
ds = xarray.open_dataset(ic, crs='EPSG:32610', scale=300, engine='ee', geometry=ee_image.geometry())
image = ds.isel(time=0).rename({'Y': 'y', 'X': 'x'}).transpose('y', 'x').rio.write_crs("EPSG:32610")
image['B1'].plot.imshow(figsize=(10, 10))

ee_image = ee.Image("GOOGLE/DYNAMICWORLD/V1/20220119T185709_20220119T185710_T10SFG")
ic = ee.ImageCollection(ee_image)
ds = xarray.open_dataset(ic, crs='EPSG:32610', scale=50, engine='ee', geometry=ee_image.geometry())
image = ds.isel(time=0).rename({'Y': 'y', 'X': 'x'}).transpose('y', 'x').rio.write_crs("EPSG:32610")
image['label'].plot.imshow(figsize=(10, 10))

ee_image = ee.Image("GOOGLE/DYNAMICWORLD/V1/20220119T185709_20220119T185710_T10SFG")
ic = ee.ImageCollection(ee_image)

roi_coords = [(-121.1100533594971, 37.64753404740281),
 (-121.10776475768417, 37.73761854601581),
 (-121.22120321499558, 37.73938555546548),
 (-121.22335479806203, 37.649295359770846),
 (-121.1100533594971, 37.64753404740281)]
roi_crs = 'EPSG:4326'
bbox_geometry = ee.Geometry.Polygon(roi_coords, proj=roi_crs)

ds = xarray.open_dataset(ic, crs='EPSG:32610', scale=10, engine='ee', geometry=bbox_geometry)
image = ds.isel(time=0).rename({'Y': 'y', 'X': 'x'}).transpose('y', 'x').rio.write_crs("EPSG:32610")

image['label'].plot.imshow(figsize=(10, 10))
original_geometry_utm = gpd.GeoDataFrame(
    {"geometry": [Polygon(roi_coords)]}, crs='EPSG:4326'
    ).to_crs('EPSG:32610')
original_geometry_utm.plot(ax=plt.gca(), edgecolor='red', facecolor='none', linewidth=4)

After the proposed FIX

The described fix was implemented for 2 of 3 possible cases handling geometry. According to the docs:
"geometry (optional): Specify an ee.Geometry to define the regional bounds when opening the data. When not set, the bounds are defined by the CRS's 'area_of_use boundaries. If those aren't present, the bounds re derived from the geometry of the first image of the collection." The proposed fix is applied only for the geometry provided by the user or obtained in the first image. if the bounds are defined by the CRS's 'area_of_use, the current behavior was kept (get bounds and later reproject them).

[jdbcode]

Test the example once #275 is included in a release.

[jdbcode]

We believe the v0.1.0 refactor has resolved the edge artifact and "black strip" issues described in this thread. As part of our due diligence, we recently validated this behavior against the updated API.

Why the problem existed:
Earlier versions relied on a "geometry-first" approach. The engine often calculated bounding boxes in the geometry's native CRS (usually EPSG:4326) and then re-projected those bounds to the user's target CRS. This reprojection process caused edge deformations and sub-pixel misalignments, which led to the zero-value gaps because the requested window didn't perfectly map to the underlying Earth Engine pixel grid.

How v0.1.0 addresses it:
The engine now uses an explicit "grid-first" approach via xee.helpers.fit_geometry.

Target-CRS Math: It calculates the bounding box and affine transform locally, directly in the target CRS, before requesting any data from Earth Engine.

Snapped Grids: It creates an integer-snapped affine grid (e.g., locking 10m or 300m pixels to clean global integer coordinates) and uses explicit tuple scales (e.g., (10, -10)) to strictly maintain a North-up raster orientation.

This ensures the requested grid perfectly aligns with the data, eliminating edge gaps. (Note: Because Xarray coordinates are now strictly center-referenced, you will see a mathematically correct "half-pixel bleed" when the raster is overlaid with raw floating-point vector boundaries).

Below are three migrated scripts demonstrating the fix in v0.1.0 for custom geometries and native image extents:

Landsat Image Extent (Downsampled to 300m)

import ee
import xarray as xr
import shapely.geometry
from xee import helpers
import matplotlib.pyplot as plt

# 1. Define Image and Collection
ee_image = ee.Image("LANDSAT/LC08/C02/T1_TOA/LC08_044034_20140318")
ic = ee.ImageCollection(ee_image)

# 2. Extract EE Geometry to a Shapely Object
footprint_dict = ee_image.geometry().getInfo()
roi_shapely = shapely.geometry.shape(footprint_dict)

# 3. Generate Grid Parameters
grid_params = helpers.fit_geometry(
    geometry=roi_shapely,
    geometry_crs='EPSG:4326',  # EE geometries are WGS84 by default
    grid_crs='EPSG:32610',     # Target CRS
    grid_scale=(300, -300)     # Target downsampled scale in meters
)

# 4. Open Dataset and Plot
ds = xr.open_dataset(ic, engine='ee', **grid_params)
image = ds.isel(time=0).rio.write_crs("EPSG:32610")

image['B1'].plot.imshow(figsize=(10, 10))
plt.show()

Dynamic World Image Extent (Downsampled to 50m)

import ee
import xarray as xr
import shapely.geometry
from xee import helpers
import matplotlib.pyplot as plt

# 1. Define Image and Collection
ee_image = ee.Image("GOOGLE/DYNAMICWORLD/V1/20220119T185709_20220119T185710_T10SFG")
ic = ee.ImageCollection(ee_image)

# 2. Extract EE Geometry to a Shapely Object
footprint_dict = ee_image.geometry().getInfo()
roi_shapely = shapely.geometry.shape(footprint_dict)

# 3. Generate Grid Parameters
grid_params = helpers.fit_geometry(
    geometry=roi_shapely,
    geometry_crs='EPSG:4326',  # EE geometries are WGS84 by default
    grid_crs='EPSG:32610',     # Target CRS
    grid_scale=(50, -50)       # Target downsampled scale in meters
)

# 4. Open Dataset and Plot
ds = xr.open_dataset(ic, engine='ee', **grid_params)
image = ds.isel(time=0).rio.write_crs("EPSG:32610")

image['label'].plot.imshow(figsize=(10, 10))
plt.show()

Custom Coordinates Box (10m Scale)

import ee
import xarray as xr
from shapely.geometry import Polygon
from xee import helpers
import matplotlib.pyplot as plt

# 1. Define Image and Collection
ee_image = ee.Image("GOOGLE/DYNAMICWORLD/V1/20220119T185709_20220119T185710_T10SFG")
ic = ee.ImageCollection(ee_image)

# 2. Define Custom ROI as a Shapely Object
roi_coords = [
    (-121.1100533594971, 37.64753404740281),
    (-121.10776475768417, 37.73761854601581),
    (-121.22120321499558, 37.73938555546548),
    (-121.22335479806203, 37.649295359770846),
    (-121.1100533594971, 37.64753404740281)
]
roi_shapely = Polygon(roi_coords)

# 3. Generate Grid Parameters
grid_params = helpers.fit_geometry(
    geometry=roi_shapely,
    geometry_crs='EPSG:4326',  # Native CRS of the coordinates
    grid_crs='EPSG:32610',     # Target CRS
    grid_scale=(10, -10)       # Target scale in meters (negative Y for North-up)
)

# 4. Open Dataset and Plot
ds = xr.open_dataset(ic, engine='ee', **grid_params)
image = ds.isel(time=0).rio.write_crs("EPSG:32610")

image['label'].plot.imshow(figsize=(10, 10))
plt.show()