Feature: support efficient per-cell capacity factor time series extraction

[MaykThewessen]

Summary

When generating capacity factor profiles for many individual grid cells (e.g., per-bus profiles for a power system model), the current recommended workflow involves creating a separate single-pixel Cutout for each location. This results in O(N) file I/O operations and O(N) calls to convert_wind/convert_pv, which is very slow for hundreds of locations.

We found that calling convert_wind and convert_pv directly on a multi-cell dataset returns per-cell (time, y, x) DataArrays, enabling fully vectorized computation — a ~50x speedup for our use case (345 bus locations, 152 unique grid pixels → 85 seconds for 3 weather years instead of ~60 minutes).

Current workaround

from atlite.convert import convert_wind, convert_pv, get_windturbineconfig, get_solarpanelconfig, get_orientation

# 1. Open cutout and select sub-region covering all locations
ds = xr.open_dataset(cutout_path, chunks="auto")
sub = ds.sel(x=slice(x_min, x_max), y=slice(y_min, y_max)).compute()

# 2. Vectorized conversion — ALL grid cells at once
turbine = get_windturbineconfig("Vestas_V112_3MW")
wind_cf = convert_wind(sub, turbine, interpolation_method="logarithmic")
# → DataArray with dims (time, y, x), shape e.g. (8760, 13, 16)

panel = get_solarpanelconfig("CSi")
orientation = get_orientation("latitude_optimal")
solar_cf = convert_pv(sub, panel, orientation, tracking=None)
# → DataArray with dims (y, time, x)

# 3. Extract per-location profiles via coordinate lookup
meeden_wind = wind_cf.sel(x=6.9, y=53.1, method="nearest").values  # shape (8760,)

This works but relies on internal functions (convert_wind, convert_pv) that are not part of the documented public API.

Proposal

I see two possible improvements:

Option A: Document the capacity_factor_timeseries=True workflow for multi-cell cutouts

Looking at convert_and_aggregate, when capacity_factor_timeseries=True is used without layout/shapes/matrix, it already returns the per-cell DataArray directly (the no_args path at line 138-142). This could be documented as the recommended approach for per-cell profile extraction:

cutout = atlite.Cutout(path="sub_region.nc")
# Returns per-cell (time, y, x) DataArray when no layout/shapes/matrix given
wind_cf = cutout.wind(turbine="Vestas_V112_3MW", capacity_factor_timeseries=True)
solar_cf = cutout.pv(panel="CSi", orientation="latitude_optimal", capacity_factor_timeseries=True)

# Extract individual locations
meeden = wind_cf.sel(x=6.9, y=53.1, method="nearest")

Option B: Add a convenience method for multi-point extraction

A new Cutout.profiles_at_points() method that takes a DataFrame of (lat, lon) coordinates and returns per-location time series:

locations = pd.DataFrame({
    "name": ["Meeden", "Borssele"],
    "x": [6.9, 3.0],
    "y": [53.1, 51.7],
})
profiles = cutout.profiles_at_points(locations, turbine="Vestas_V112_3MW")
# → DataFrame with columns ["Meeden", "Borssele"] and DatetimeIndex

Use case

This is common in power system modeling where each bus/substation needs location-specific renewable capacity factors. In PyPSA-based models, generators at different buses should have distinct wind/solar profiles reflecting local weather conditions, not a single country-level average.

Performance comparison

For our model (345 bus locations across 8 European countries, SARAH3+ERA5 cutouts):

Approach	Time per year	Method
Per-pixel Cutout (current)	~20-30 min	152 × Cutout + cutout.wind()
Vectorized convert_wind/convert_pv	~28 seconds	8 country sub-regions, parallel

The vectorized approach is ~50x faster because it eliminates per-pixel file I/O and leverages numpy's vectorized operations across the entire grid.

[euronion]

I don't understand the problem statement. You can pass a 3D matrix with a 2D capacity layout and an index to convert_and_aggregate already. That's how we create e.g. profiles for each bus region in PyPSA-Eur.

[MaykThewessen]

Thanks for the pointer! You're right — the matrix path with convert_and_aggregate works well for regional aggregation, and I can see how PyPSA-Eur uses indicator matrices from Voronoi regions to produce per-bus profiles.

My issue is about a slightly different use case: point extraction — where you have (lat, lon) coordinates and want the capacity factor at the nearest grid cell, without needing to construct region polygons or indicator matrices first.

I looked into whether switching to Voronoi-based regional averaging would be better for my use case (345 bus locations across 8 countries), but at ERA5's 0.3° resolution (~1000 km² per cell), the median Voronoi region for my buses is only ~250-400 km² — smaller than a single weather pixel. So regional averaging wouldn't meaningfully change the profiles; the weather data resolution is the bottleneck, not the extraction method.

For point extraction, the "no_args" path already works cleanly:

cf = cutout.wind(turbine="Vestas_V112_3MW", capacity_factor_timeseries=True)
# Returns (time, y, x) DataArray — full grid, no aggregation
meeden = cf.sel(x=6.9, y=53.1, method="nearest")

This path exists (lines 137-143 in convert.py) but isn't documented. Would a documentation PR be welcome covering this workflow alongside the existing matrix-based approach?

[FabianHofmann]

This path exists (lines 137-143 in convert.py) but isn't documented. Would a documentation PR be welcome covering this workflow alongside the existing matrix-based approach?

yes, always welcome! and I see that this is a gap in the documentation (only docstrings atm).

The sparse operation that only applies the conversion at selected points, I actually think that is a nice idea! My intuition says we should further overload the conversion functions (and later in future break it down). cutout.profiles_at_points as you wrote it, does not specify the kind of conversion. so something like cutout.pv(points=...) could do. what do you think?

[MaykThewessen]

Great idea cutout.pv(points=...) is much cleaner than a separate method. It keeps the API surface minimal and the intent obvious.

The accessor pattern (cutout.pv.profiles_at_points) is interesting for future extensibility, but for the immediate use case I think a points parameter on the existing conversion methods is the simplest path:

• it's backward-compatible, discoverable, and avoids introducing a new abstraction layer.

Under the hood it could work like:

1. Call the existing conversion with capacity_factor_timeseries=True (no matrix/layout)
2. Use .sel(x=..., y=..., method="nearest") on the result for each point
3. Return a DataFrame with point names as columns

Happy to open a separate issue to discuss the API design for this and split it from the docs PR. I'll also update the PR to be docs-only as suggested in #481.