Cython Implementation of the NREL SPA

[leaver2000]

This is not so much a Feature Request, rather a conversation opener. I have written a cython implementation of the NREL SPA that was heavily influenced by the pvlib.spa.solar_position_numpy and I just wanted to share and gauge interest. I've linked the repo below and provided some timeit's

https://github.com/leaver2000/fast_spa/blob/master/tests/benchmark.ipynb

import numpy as np
import pvlib.spa
import fast_spa

def spa(
    obj,
    lat,
    lon,
    elevation=0.0,
    pressure=1013.25,
    temp=12.0,
    refraction=0.5667,
):
    delta_t = fast_spa.pe4dt(obj)
    dt = np.asanyarray(obj, dtype="datetime64[ns]")
    unix_time = dt.astype(np.float64) // 1e9

    return np.stack(
        [
            np.stack(
                pvlib.spa.solar_position_numpy(
                    ut,
                    lat=lat,
                    lon=lon,
                    elev=elevation,
                    pressure=pressure,
                    temp=temp,
                    delta_t=delta_t[i : i + 1],
                    atmos_refract=refraction,
                    numthreads=None,
                    sst=False,
                    esd=False,
                )[:-1]
            )
            for i, ut in enumerate(unix_time)
        ],
        axis=1,
    )

lats = np.linspace(30, 31, 100)
lons = np.linspace(-80, -79, 100)
lats, lons = np.meshgrid(lats, lons)
datetime_obj = (
    np.arange("2023-01-01", "2023-01-02", dtype="datetime64[h]")
    .astype("datetime64[ns]")
    .astype(str)
    .tolist()
)

%timeit fast_spa.fast_spa(datetime_obj, lats, lons)
# 29.1 ms ± 299 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit spa(datetime_obj, lats, lons)
# 65 ms ± 498 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

[kandersolar]

I am surprised that the runtime difference is only ~2x! Have you compared runtime with pvlib's numba-based SPA function (I guess single-threaded, to be apples-to-apples?).

In general I often daydream about speeding up a few particular parts of pvlib (not just the SPA) using cython or similar, but it has to be worth the cost of added difficulty in maintenance and distribution. For distribution, we'd either have to supply pre-built wheels for each platform, or require cython and a compiler during installation, or make this new code optional, right? And for maintenance, we would need to become proficient in cython, which isn't too difficult, but also not trivial. In those respects, numba seems like less of an investment than cython does, for comparable benefit.

For solar position calculations specifically, I wonder if the developer time would be better spent on a more efficient algorithm (https://github.com/pvlib/pvlib-python/discussions/1308) rather than a more efficient implementation of the SPA.

[leaver2000]

I am surprised that the runtime difference is only ~2x! Have you compared runtime with pvlib's numba-based SPA function (I guess single-threaded, to be apples-to-apples?).

The performance boost is dependent on the size of the grid and time period. Which I suspect is closely tied to the terms tables. A larger gap here over a longer period.

lats = np.linspace(30, 31, 5)
lons = np.linspace(-80, -79, 5)
lats, lons = np.meshgrid(lats, lons)
datetime_obj = (
    np.arange("2023-01-01", "2023-01-07", dtype="datetime64[h]")
    .astype("datetime64[ns]")
    .astype(str)
    .tolist()
)
# 4.44 ms ± 24.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
# 166 ms ± 765 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

For solar position calculations specifically, I wonder if the developer time would be better spent on a more efficient algorithm (#1308) rather than a more efficient implementation of the SPA.

I have worked with the USNO SLAC (solar lunar almanac core) in the past, mostly unit test and c bindings. I don't think it's public domain, however.

https://books.google.com/books/about/The_Solar_lunar_Almanac_Core_SLAC.html?id=ZAlGywAACAAJ

[kurt-rhee]

Probably not useful at all, but I wrote a Rust implementation of the NREL 2008 algorithm. I get similar benchmarking numbers to @leaver2000's C implementation (about 2x speedup compared to method='nrel_numpy') over 1 year of hourly time-steps.

Happy to open source it and write python bindings for pvlib if anybody is interested.

[leaver2000]

Probably not useful at all, but I wrote a Rust implementation of the NREL 2008 algorithm. I get similar benchmarking numbers to @leaver2000's C implementation (about 2x speedup compared to method='nrel_numpy') over 1 year of hourly time-steps.

Happy to open source it and write python bindings for pvlib if anybody is interested.

I'm interested to take a look. In the past I've built some very basic rust functions with python/numpy bindings. I'm sure I could learn something from your work.