`compute_shape_availability` triggers numpy failure: Python integer 255 out of bounds for int8

[irm-codebase]

Version Checks (indicate both or one)

• [x] I have confirmed this bug exists on the lastest release of Atlite.

• [ ] I have confirmed this bug exists on the current master branch of Atlite.

Issue Description

While trying to reproduce the land availability example, the following error is returned.

It seems like this is caused by a depreciated functionality in numpy: https://github.com/numpy/numpy/issues/26596

Either numpy or something else needs to be pinned in the environment, maybe?

File ~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/atlite/gis.py:586, in ExclusionContainer.compute_shape_availability(self, geometry, dst_transform, dst_crs, dst_shape)
    [582](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/atlite/gis.py:582)     return shape_availability_reprojected(
    [583](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/atlite/gis.py:583)         geometry, self, dst_transform, dst_crs, dst_shape
    [584](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/atlite/gis.py:584)     )
    [585](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/atlite/gis.py:585) else:
--> [586](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/atlite/gis.py:586)     return shape_availability(geometry, self)
...
    [494](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/numpy/ma/core.py:494)             err_msg = "Cannot convert fill_value %s to dtype %s"
--> [495](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/numpy/ma/core.py:495)             raise TypeError(err_msg % (fill_value, ndtype)) from e
    [496](https://file+.vscode-resource.vscode-cdn.net/home/ivanruizmanuel/Documents/git/ec_modules/tmp/~/miniforge3/envs/ec-atlite-pv-series/lib/python3.12/site-packages/numpy/ma/core.py:496) return np.array(fill_value)

TypeError: Cannot convert fill_value 255 to dtype int8

Reproducible Example

import atlite
import xarray as xr
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.plot import show
from atlite.gis import shape_availability, ExclusionContainer
import cartopy.io.shapereader as shpreader
import numpy as np

shp = shpreader.Reader(
    shpreader.natural_earth(
        resolution="110m", category="cultural", name="admin_0_countries"
    )
)


countries = shp.records()
country = next(countries)
country.attributes

records = list(
    filter(lambda r: r.attributes["NAME"] in ["Serbia", "Croatia", "Bosnia and Herz.", "Montenegro"], shp.records())
)
shapes = (
    gpd.GeoDataFrame([{**r.attributes, "geometry": r.geometry} for r in records])
    .set_index("NAME")
    .set_crs(4236)
)
shapes.plot()

b_min = shapes.total_bounds[:2] -1
b_max = shapes.total_bounds[2:] + 1
cutout = atlite.Cutout(
    "balkans", module="era5", bounds=np.concatenate([b_min, b_max]), time=slice("2013-01-01", "2013-01-02")
)

plt.rc("figure", figsize=[10, 7])
fig, ax = plt.subplots()
shapes.plot(ax=ax)
cutout.grid.plot(ax=ax, edgecolor="grey", color="None")

CORINE = "tmp/U2018_CLC2018_V2020_20u1.tif"
excluder = ExclusionContainer()
excluder.add_raster(CORINE, codes=range(20))


croatia = shapes.loc[["Croatia"]].geometry.to_crs(excluder.crs)

masked, transform = excluder.compute_shape_availability(croatia)

Expected Behavior

Reproducing the example should be possible.

Notes: the .tif was downloaded from the CORINE website.

Installed Versions

0.2.13

[irm-codebase]

After digging around, I confirmed that the issue is an unpinned numpy version. Pinning numpy = 1.23 will solve the issue.

[fneum]

https://github.com/PyPSA/atlite/blob/16c02fd4aa1816f73d96ed5ae52b1be69f671b18/atlite/gis.py#L378

Could you try setting this to int16? That might solve it.

(anyway I observed a potential performance regression with numpy 2 and we should pin the version for now)

[irm-codebase]

@fneum did the test, int16 fixed it!

[fneum]

Picked this one up again to tackle numpy>=2 compatibility.

There was actually no problem per se with numpy here. The new version just handles datatype overflows more rigorously.

The file in question reads as int8 (because it has values within -128 and 127):

import rasterio
rasterio.open("U2018_CLC2018_V2020_20u1.tif").read(1)

array([[-128, -128, -128, ..., -128, -128, -128],
       [-128, -128, -128, ..., -128, -128, -128],
       [-128, -128, -128, ..., -128, -128, -128],
       ...,
       [-128, -128, -128, ..., -128, -128, -128],
       [-128, -128, -128, ..., -128, -128, -128],
       [-128, -128, -128, ..., -128, -128, -128]],
      shape=(5662, 3404), dtype=int8)

The ExclusionContainer assumes for raster data a default for "nodata" values of 255 (which does not fit in int8).

The correct (and recommended way) to add this raster to the exclusion container is with the correct "nodata" value of the dataset, which is -128:

from atlite.gis import ExclusionContainer
excluder = ExclusionContainer(crs=3035)
excluder.add_raster("U2018_CLC2018_V2020_20u1.tif", codes=range(20), nodata=-128)