Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add dtype validation when rescale=True #209

Merged
merged 13 commits into from
Sep 19, 2023
1 change: 1 addition & 0 deletions .pdm-python
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
/Users/gabe/dev/stackstac/.venv/bin/python
37 changes: 32 additions & 5 deletions stackstac/prepare.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,9 @@
from .stac_types import ItemSequence
from . import accumulate_metadata, geom_utils

ASSET_TABLE_DT = np.dtype([("url", object), ("bounds", "float64", 4), ("scale_offset", "float64", 2)])
ASSET_TABLE_DT = np.dtype(
[("url", object), ("bounds", "float64", 4), ("scale_offset", "float64", 2)]
)


class Mimetype(NamedTuple):
Expand Down Expand Up @@ -64,6 +66,8 @@ def prepare_items(
bounds: Optional[Bbox] = None,
bounds_latlon: Optional[Bbox] = None,
snap_bounds: bool = True,
rescale: bool = True,
dtype: np.dtype = np.dtype("float64"),
) -> Tuple[np.ndarray, RasterSpec, List[str], ItemSequence]:

if bounds is not None and bounds_latlon is not None:
Expand Down Expand Up @@ -144,7 +148,7 @@ def prepare_items(
asset_bbox = asset.get("proj:bbox", item_bbox)
asset_shape = asset.get("proj:shape", item_shape)
asset_transform = asset.get("proj:transform", item_transform)
raster_bands = asset.get('raster:bands')
raster_bands = asset.get("raster:bands")

if raster_bands is not None:
if len(raster_bands) != 1:
Expand All @@ -154,12 +158,31 @@ def prepare_items(
"If you don't care about this asset, you can skip it by giving a list "
"of asset IDs you *do* want in `assets=`, and leaving this one out."
)
asset_scale = raster_bands[0].get('scale', 1)
asset_offset = raster_bands[0].get('offset', 0)
asset_scale = raster_bands[0].get("scale", 1)
asset_offset = raster_bands[0].get("offset", 0)
else:
asset_scale = 1
asset_offset = 0

if rescale:
if not np.can_cast(asset_scale, dtype):
raise ValueError(
f"`rescale=True`, but safe casting cannot be completed between "
f"asset scale value {asset_scale} and output dtype {dtype}.\n"
f"To continue using `{dtype=}`, pass `rescale=False` to retrieve "
"data in its raw, unscaled values. Or, if you want rescaled "
"values, pass a different `dtype=` (typically `float`)."
)

if not np.can_cast(asset_offset, dtype):
raise ValueError(
f"`rescale=True`, but safe casting cannot be completed between "
f"asset offset value {asset_offset} and output dtype {dtype}.\n"
f"To continue using `{dtype=}`, pass `rescale=False` to retrieve "
"data in its raw, unscaled values. Or, if you want rescaled "
"values, pass a different `dtype=` (typically `float`)."
)

asset_affine = None

# Auto-compute CRS
Expand Down Expand Up @@ -339,7 +362,11 @@ def prepare_items(
continue

# Phew, we figured out all the spatial stuff! Now actually store the information we care about.
asset_table[item_i, asset_i] = (asset["href"], asset_bbox_proj, (asset_scale, asset_offset))
asset_table[item_i, asset_i] = (
asset["href"],
asset_bbox_proj,
(asset_scale, asset_offset),
)
# ^ NOTE: If `asset_bbox_proj` is None, NumPy automatically converts it to NaNs

# At this point, everything has been set (or there was as error)
Expand Down
2 changes: 2 additions & 0 deletions stackstac/stack.py
Original file line number Diff line number Diff line change
Expand Up @@ -293,6 +293,8 @@ def stack(
bounds=bounds,
bounds_latlon=bounds_latlon,
snap_bounds=snap_bounds,
rescale=rescale,
dtype=dtype,
)
arr = items_to_dask(
asset_table,
Expand Down