PyPSA / atlite

Atlite: A Lightweight Python Package for Calculating Renewable Power Potentials and Time Series
https://atlite.readthedocs.io
265 stars 89 forks source link

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

Open irm-codebase opened 1 month ago

irm-codebase commented 1 month ago

Version Checks (indicate both or one)

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 commented 1 month ago

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

fneum commented 1 month ago

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 commented 1 month ago

@fneum did the test, int16 fixed it!