fill_value incompatible with output dtype (stackstac.stack)

[ZZMitch]

Hello,

I recently updated to stackstac 0.5.1 and found a bug/incorrect documentation with stackstac.stack() and the fill_value and dtype parameters. This is related to #250 and tracks back to numpy.can_cast().

Take this example, where I want to grab just the fmask band from HLS L30 over an area at the smallest possible dtype ('unit8'):

Example:

import stackstac as ss

# Temporal bounds
year = 2018
start = str(year) + '-01-01'
end =  str(year) + '-12-31'

# Spatial bounds
bboxLL = (np.float64(-149.53623682903373), np.float64(68.55062861124203), np.float64(-148.30140593046852), np.float64(69.2185427988887))
epsg = 32606
bbox = (np.float64(399535.48481308576), np.float64(7606266.209112197), np.float64(446924.35747663863), np.float64(7679311.477803782))

# STAC settings
gdalEnv = ss.DEFAULT_GDAL_ENV.updated(dict(
    GDAL_DISABLE_READDIR_ON_OPEN = 'TRUE',
    GDAL_HTTP_COOKIEFILE = os.path.expanduser('~/cookies.txt'),
    GDAL_HTTP_COOKIEJAR = os.path.expanduser('~/cookies.txt'),
    GDAL_HTTP_MAX_RETRY = 10,
    GDAL_HTTP_RETRY_DELAY = 15, 
    GDAL_HTTP_UNSAFESSL = 'YES')) # Will pass to Dask cluster workers

catalog = pc.Client.open('https://cmr.earthdata.nasa.gov/stac/LPCLOUD')

# Create stack
items = catalog.search(bbox = bboxLL, datetime = f'{start}/{end}', collections = ['HLSL30.v2.0'], limit = 100).item_collection() 

stack = ss.stack(items, assets = ['Fmask'], epsg = epsg, resolution = 30, bounds = bbox, dtype = 'uint8', fill_value = 255, #np.uint8(255)
    rescale = False, gdal_env = gdalEnv)
stack

The above will fail with:
ValueError: The fill_value 255 is incompatible with the output dtype uint8. Either use dtype='int64', or pick a different fill_value.
Providing another value between 0-255 results in the same error.

Providing a fill_value of np.uint8(255) works and creates a stack with the correct dtype (verified with stack.load():

Checking numpy.can_cast(), which the above error points to, we see:

import numpy as np
np.can_cast(type(255), 'uint8') # False
np.can_cast(type(np.uint8(255)), 'uint8') # True

Previously, you could just provide an int like 255. So I am not sure if this is a bug or the documentation for stackstac.stack() needs to be updated to indicate that fill_value should specifically be a numpy dtype.

[fgassert]

Hi, I've also run into this issue.

I think this should correctly test if the fill_value will fit into the dtype in both numpy 1.x and 2.x.

np.can_cast(np.min_scaler_type(fill_value), dtype)

E.g.

import numpy as np
np.can_cast(type(0), 'float32') # False
np.can_cast(np.min_scaler_type(0), 'float32') # True

Thanks!

[fkroeber]

To provide a bit more of related background: The new behaviour of np.can_cast() as introduced in numpy>2 does not apply any value-based logic for 0-D arrays anymore (further background info here). I guess it remains a design choice if...

a) ...we want to stick with the guiding principle of "Values must never influence result type" underlying the change in numpy>2, which would implicate that users actually have to specify the fill value as a numpy scalar with a defined dtype (as carried out by @ZZMitch)
b) ...we prioritise the ease-of-use of simple Python scalars and find some workaround (as proposed by @fgassert)

Regarding b): np.min_scalar_type is a good starting point to determine the type for Python scalars but unfortunately still doesn't provide us with the same evaluation logic as in numpy 1. While the example above assessing compabatility of 0 and float32 works the same in both python versions, the following example doesn't yield the same results:

import numpy as np

# numpy version 1.x 
np.can_cast(0, "int8") # evaluates to True

# numpy version 2.x 
np.can_cast(type(0), "int8") # evaluates to False since type(0) = int
np.can_cast(np.min_scalar_type(0), "int8") # evaluates to False since min_scalar_type(0) = uint8