Skip to content

mcda.constrain() drops raster attrs when exclusion masks are applied #3147

@brendancol

Description

@brendancol

Describe the bug

constrain() in xrspatial/mcda/constrain.py returns a DataArray with empty attrs whenever exclude contains at least one mask. The loop rewrites the result with xr.where(mask, fill, result), and xr.where takes its attrs from the first value argument, which is the scalar fill. So res, crs, and nodatavals are gone from the constrained surface. With exclude=[] the attrs survive, since the function just returns a copy.

import numpy as np
import xarray as xr
from xrspatial.mcda import constrain

suit = xr.DataArray(
    np.full((2, 3), 0.5), dims=["y", "x"],
    attrs={"res": (1.0, 1.0), "crs": "EPSG:32611", "nodatavals": (-9999.0,)},
)
mask = xr.DataArray(np.zeros((2, 3), dtype=bool), dims=["y", "x"])
print(constrain(suit, exclude=[]).attrs)      # attrs kept
print(constrain(suit, exclude=[mask]).attrs)  # {} -- attrs gone

Expected behavior

The output keeps the input's attrs regardless of how many masks are passed. constrain is usually the last step before exporting the suitability surface, so losing crs/res here means the exported GeoTIFF has no georeferencing.

Additional context

Affects numpy, dask+numpy, and dask+cupy the same way (verified on all three). Bare cupy currently fails earlier inside xr.where because of the known cupy 13.6 / xarray 2025.12 incompatibility, which is a separate issue.

Found by the metadata propagation sweep against the mcda module.

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