`solshade.solar`¶

compute_solar_ephem ¶

compute_solar_ephem(lat, lon, startutc=None, stoputc=None, timestep=3600, cache_dir='data/skyfield')

Compute the Sun’s apparent altitude/azimuth and the unit ENU direction vector over a time range at a given site.

The ENU vector is derived directly from the local topocentric (alt, az) returned by Skyfield—avoiding intermediate ECEF rotations and their numerical pitfalls.

Parameters:

Name	Type	Description	Default
`lat`	`float`	Geographic latitude in degrees (+N).	required
`lon`	`float`	Geographic longitude in degrees (+E).	required
`startutc`	`datetime`	UTC start time. If None, uses start of current UTC year. If naive, assumed UTC.	`None`
`stoputc`	`datetime`	UTC stop time. If None, defaults to one year after `startutc`. Must be strictly after `startutc`. If naive, assumed UTC.	`None`
`timestep`	`int`	Sampling step in seconds. Must be > 0.	`3600`
`cache_dir`	`str or Path`	Skyfield cache directory (timescale + DE440s kernel). Default: `./data/skyfield`.	`'data/skyfield'`

Returns:

Name	Type	Description
`times_utc`	`ndarray of datetime64[ns], shape (N,)`	UTC timestamps for each sample.
`alt_deg`	`ndarray of float, shape (N,)`	Apparent altitude (degrees). No atmospheric refraction applied.
`az_deg`	`ndarray of float, shape (N,)`	Apparent azimuth (degrees), clockwise from true north, wrapped to [0, 360).
`dist_au`	`ndarray of float, shape (N,)`	Apparent distance between observer and Sun in astronomical units
`enu_unit`	`ndarray of float, shape (N, 3)`	Unit vectors pointing from the site toward the Sun in local ENU coordinates (columns: E, N, U). Each row has norm ~= 1.

Raises:

Type	Description
`ValueError`	If `timestep` <= 0 or if `stoputc` <= `startutc`.
`FileNotFoundError`	If Skyfield cache files (timescale or ephemeris) are missing.

Notes

Uses the compact DE440s kernel.
Apparent alt/az include light-time/aberration, no refraction.
ENU unit vector is computed via:

E = cos(alt) * sin(az) N = cos(alt) * cos(az) U = sin(alt)

where azimuth is CW from north (Skyfield’s default for topocentric alt/az).

Logging

INFO: parameters (lat/lon, duration, timestep).
DEBUG: step counts, times constructed, ENU unit vector shape.

Source code in src/solshade/solar.py

def compute_solar_ephem(
    lat: float,
    lon: float,
    startutc: Optional[datetime] = None,
    stoputc: Optional[datetime] = None,
    timestep: int = 3600,
    cache_dir: Optional[Union[str, Path]] = "data/skyfield",
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Compute the Sun’s apparent altitude/azimuth and the **unit ENU direction vector**
    over a time range at a given site.

    The ENU vector is derived *directly* from the local topocentric (alt, az) returned
    by Skyfield—avoiding intermediate ECEF rotations and their numerical pitfalls.

    Parameters
    ----------
    lat : float
        Geographic latitude in degrees (+N).
    lon : float
        Geographic longitude in degrees (+E).
    startutc : datetime, optional
        UTC start time. If None, uses start of current UTC year. If naive, assumed UTC.
    stoputc : datetime, optional
        UTC stop time. If None, defaults to one year after ``startutc``.
        Must be strictly after ``startutc``. If naive, assumed UTC.
    timestep : int, default=3600
        Sampling step in seconds. Must be > 0.
    cache_dir : str or Path, optional
        Skyfield cache directory (timescale + DE440s kernel). Default: ``./data/skyfield``.

    Returns
    -------
    times_utc : ndarray of datetime64[ns], shape (N,)
        UTC timestamps for each sample.
    alt_deg : ndarray of float, shape (N,)
        Apparent altitude (degrees). No atmospheric refraction applied.
    az_deg : ndarray of float, shape (N,)
        Apparent azimuth (degrees), clockwise from true north, wrapped to [0, 360).
    dist_au : ndarray of float, shape (N,)
        Apparent distance between observer and Sun in astronomical units
    enu_unit : ndarray of float, shape (N, 3)
        Unit vectors pointing from the site toward the Sun in local ENU coordinates
        (columns: E, N, U). Each row has norm ~= 1.

    Raises
    ------
    ValueError
        If ``timestep`` <= 0 or if ``stoputc`` <= ``startutc``.
    FileNotFoundError
        If Skyfield cache files (timescale or ephemeris) are missing.

    Notes
    -----
    - Uses the compact DE440s kernel.
    - Apparent alt/az include light-time/aberration, **no refraction**.
    - ENU unit vector is computed via:

          E = cos(alt) * sin(az)
          N = cos(alt) * cos(az)
          U = sin(alt)

      where azimuth is CW from north (Skyfield’s default for topocentric alt/az).

    Logging
    -------
    - INFO: parameters (lat/lon, duration, timestep).
    - DEBUG: step counts, times constructed, ENU unit vector shape.
    """

    def ensure_utc(dt: datetime) -> datetime:
        if dt.tzinfo is None:
            return dt.replace(tzinfo=timezone.utc)
        return dt.astimezone(timezone.utc)

    start_of_year_utc = datetime(datetime.now(timezone.utc).year, 1, 1, tzinfo=timezone.utc)
    start_dt = ensure_utc(startutc) if startutc else start_of_year_utc
    stop_dt = ensure_utc(stoputc) if stoputc else (start_of_year_utc + timedelta(days=365))

    if timestep <= 0:
        log.error("timestep must be positive")
        raise ValueError("timestep must be a positive number of seconds")
    if stop_dt <= start_dt:
        log.error("stoputc must be after startutc")
        raise ValueError("stoputc must be after startutc")

    total_seconds = int((stop_dt - start_dt).total_seconds())
    steps = total_seconds // timestep
    log.info(f"Computing solar ephemeris lat={lat:.3f}, lon={lon:.3f}, span={stop_dt - start_dt}, steps={steps}")

    offsets = np.arange(0, steps * timestep + 1, timestep, dtype="int64")
    times_py = [start_dt + timedelta(seconds=int(s)) for s in offsets]
    times_utc = np.array([t.replace(tzinfo=None) for t in times_py], dtype="datetime64[ns]")

    sun, earth, ts = load_sun_ephemeris(cache_dir)
    t = ts.from_datetimes(times_py)

    observer = wgs84.latlon(latitude_degrees=lat, longitude_degrees=lon)
    apparent = (earth + observer).at(t).observe(sun).apparent()  # type: ignore
    dist_au = apparent.distance().au
    alt, az, _ = apparent.altaz()

    alt_deg = np.asarray(alt.degrees, dtype=float)
    az_deg = np.asarray(np.mod(az.degrees, 360.0), dtype=float)

    alt_rad = np.deg2rad(alt_deg)
    az_rad = np.deg2rad(az_deg)

    c_alt = np.cos(alt_rad)
    enu_unit = np.stack(
        (
            c_alt * np.sin(az_rad),  # E
            c_alt * np.cos(az_rad),  # N
            np.sin(alt_rad),  # U
        ),
        axis=1,
    ).astype(float)

    log.debug(f"Computed ENU unit vectors of shape {enu_unit.shape}")
    return times_utc, alt_deg, az_deg, dist_au, enu_unit

load_sun_ephemeris ¶

load_sun_ephemeris(cache_dir='data/skyfield')

Load the Sun and Earth ephemeris segments from the DE440s file and a Timescale.

Parameters:

Name	Type	Description	Default
`cache_dir`	`str or Path`	Directory for ephemeris and timescale cache. Defaults to ./data/skyfield.	`'data/skyfield'`

Returns:

Name	Type	Description
`sun`	`object`	Skyfield segment for the Sun (use with `.observe()`).
`earth`	`object`	Skyfield segment for the Earth (use with `earth + topos`).
`ts`	`Timescale`	Timescale object for constructing Skyfield times.

Raises:

Type	Description
`FileNotFoundError`	If required ephemeris or timescale files are missing from the cache.

Notes

Uses the compact DE440s kernel (sufficient for Sun/Earth work).
You must run once with internet access to populate the Skyfield cache, or manually place 'de440s.bsp' and timescale data in the cache_dir.

Logging

INFO: confirms successful loading of ephemeris/timescale.
ERROR: missing files raise exceptions.

Source code in src/solshade/solar.py

def load_sun_ephemeris(
    cache_dir: Optional[Union[str, Path]] = "data/skyfield",
) -> Tuple[SunSegment, EarthSegment, Timescale]:
    """
    Load the Sun and Earth ephemeris segments from the DE440s file and a Timescale.

    Parameters
    ----------
    cache_dir : str or Path, optional
        Directory for ephemeris and timescale cache. Defaults to ./data/skyfield.

    Returns
    -------
    sun : object
        Skyfield segment for the Sun (use with `.observe()`).
    earth : object
        Skyfield segment for the Earth (use with `earth + topos`).
    ts : skyfield.timelib.Timescale
        Timescale object for constructing Skyfield times.

    Raises
    ------
    FileNotFoundError
        If required ephemeris or timescale files are missing from the cache.

    Notes
    -----
    - Uses the compact DE440s kernel (sufficient for Sun/Earth work).
    - You must run once with internet access to populate the Skyfield cache,
      or manually place 'de440s.bsp' and timescale data in the cache_dir.

    Logging
    -------
    - INFO: confirms successful loading of ephemeris/timescale.
    - ERROR: missing files raise exceptions.
    """
    eph_name = "de440s.bsp"
    cache_root = Path(cache_dir) if cache_dir is not None else Path("data/skyfield")
    cache_root.mkdir(parents=True, exist_ok=True)
    loader = Loader(str(cache_root))

    log.debug(f"Loading Skyfield ephemeris from cache {cache_root}")

    # Load timescale
    try:
        ts = loader.timescale()
    except Exception as exc:
        log.error(f"Timescale not found in Skyfield cache {cache_root}")
        raise FileNotFoundError(
            f"Timescale data not found in Skyfield cache ({cache_root}). "
            "Run once with internet or manually copy the required files."
        ) from exc

    # Load ephemeris
    try:
        ephem = loader(eph_name)
    except Exception as exc:
        log.error(f"Ephemeris {eph_name} not found in {cache_root}")
        raise FileNotFoundError(
            f"Ephemeris '{eph_name}' not found in Skyfield cache ({cache_root}). "
            "Run once with internet or manually copy the BSP file."
        ) from exc

    log.info(f"Loaded ephemeris {eph_name} and timescale from {cache_root}")
    sun = ephem["sun"]
    earth = ephem["earth"]
    return sun, earth, ts

nearest_horizon_indices ¶

nearest_horizon_indices(sun_az_deg, horizon_az_deg)

Map each solar azimuth to the nearest horizon azimuth index.

Parameters:

Name	Type	Description	Default
`sun_az_deg`	`(ndarray, shape(M))`	Solar azimuths in degrees (any real values). Wrapped to [0, 360).	required
`horizon_az_deg`	`(ndarray, shape(N))`	Uniformly spaced horizon azimuths in degrees over [0, 360). Typically from `np.linspace(0, 360, n_directions, endpoint=False)`.	required

Returns:

Name	Type	Description
`idx`	`(ndarray, shape(M))`	Indices into `horizon_az_deg` of the nearest azimuth for each sun azimuth.

Logging

DEBUG: reports horizon spacing and input count.

Source code in src/solshade/solar.py

def nearest_horizon_indices(
    sun_az_deg: np.ndarray,
    horizon_az_deg: np.ndarray,
) -> np.ndarray:
    """
    Map each solar azimuth to the nearest horizon azimuth index.

    Parameters
    ----------
    sun_az_deg : np.ndarray, shape (M,)
        Solar azimuths in degrees (any real values). Wrapped to [0, 360).
    horizon_az_deg : np.ndarray, shape (N,)
        Uniformly spaced horizon azimuths in degrees over [0, 360).
        Typically from ``np.linspace(0, 360, n_directions, endpoint=False)``.

    Returns
    -------
    idx : np.ndarray, shape (M,)
        Indices into ``horizon_az_deg`` of the nearest azimuth for each sun azimuth.

    Logging
    -------
    - DEBUG: reports horizon spacing and input count.
    """
    n = horizon_az_deg.size
    step = 360.0 / n
    log.debug(f"Mapping {sun_az_deg.size} solar azimuths to {n} horizon directions (step={step:.2f}°)")
    sun_wrapped = np.mod(sun_az_deg, 360.0)
    idx = np.rint(sun_wrapped / step).astype(int) % n
    return idx

solar_envelope_by_folding ¶

solar_envelope_by_folding(times_utc, alt_deg, az_deg, smooth_n=360)

Build daily solar envelopes (min/max altitude vs azimuth) by folding a uniformly sampled time series into UTC slots and taking per-slot extrema.

Parameters:

Name	Type	Description	Default
`times_utc`	`np.ndarray of datetime64[*]`	UTC timestamps at uniform cadence covering >= 1 full day. Any datetime64 unit is accepted and internally normalized to seconds.	required
`alt_deg`	`(ndarray, shape(N))`	Solar altitude samples in degrees.	required
`az_deg`	`(ndarray, shape(N))`	Solar azimuth samples in degrees. Values are wrapped to [0, 360).	required
`smooth_n`	`int`	Number of equally spaced azimuth points (degrees) for the smoothed envelope. Must be >= 4.	`360`

Returns:

Name	Type	Description
`az_plot`	`(ndarray, shape(smooth_n + 1))`	Azimuth grid in degrees, wrapped to include 360 for closed plotting.
`alt_min_plot`	`(ndarray, shape(smooth_n + 1))`	Smoothed minimum-altitude envelope at az_plot.
`alt_max_plot`	`(ndarray, shape(smooth_n + 1))`	Smoothed maximum-altitude envelope at az_plot.

Raises:

Type	Description
`ValueError`	Non-1D or mismatched input lengths Non-uniform sampling cadence, or cadence that doesn't divide 86400 s Not enough samples for one full day <4 unique azimuth samples for cubic spline

Logging

DEBUG: input sizes, cadence, slots per day.
INFO: completion and output grid length.

Source code in src/solshade/solar.py

def solar_envelope_by_folding(
    times_utc: np.ndarray,
    alt_deg: np.ndarray,
    az_deg: np.ndarray,
    smooth_n: int = 360,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Build daily solar envelopes (min/max altitude vs azimuth) by folding a
    uniformly sampled time series into UTC slots and taking per-slot extrema.

    Parameters
    ----------
    times_utc : np.ndarray of datetime64[*]
        UTC timestamps at uniform cadence covering >= 1 full day. Any datetime64
        unit is accepted and internally normalized to seconds.
    alt_deg : np.ndarray, shape (N,)
        Solar altitude samples in degrees.
    az_deg : np.ndarray, shape (N,)
        Solar azimuth samples in degrees. Values are wrapped to [0, 360).
    smooth_n : int, default=360
        Number of equally spaced azimuth points (degrees) for the smoothed
        envelope. Must be >= 4.

    Returns
    -------
    az_plot : np.ndarray, shape (smooth_n+1,)
        Azimuth grid in degrees, wrapped to include 360 for closed plotting.
    alt_min_plot : np.ndarray, shape (smooth_n+1,)
        Smoothed minimum-altitude envelope at az_plot.
    alt_max_plot : np.ndarray, shape (smooth_n+1,)
        Smoothed maximum-altitude envelope at az_plot.

    Raises
    ------
    ValueError
        - Non-1D or mismatched input lengths
        - Non-uniform sampling cadence, or cadence that doesn't divide 86400 s
        - Not enough samples for one full day
        - <4 unique azimuth samples for cubic spline

    Logging
    -------
    - DEBUG: input sizes, cadence, slots per day.
    - INFO: completion and output grid length.
    """
    log.debug("Building solar envelope by folding")
    if alt_deg.shape != az_deg.shape:
        raise ValueError("alt_deg and az_deg must have the same shape.")
    if alt_deg.ndim != 1:
        raise ValueError("alt_deg and az_deg must be 1-D arrays.")
    if times_utc.shape[0] != alt_deg.shape[0]:
        raise ValueError("times_utc, alt_deg, az_deg lengths must match.")

    n = alt_deg.size
    if n < 2:
        raise ValueError("Need at least two samples to infer cadence.")

    t_sec = times_utc.astype("datetime64[s]")
    diffs = np.diff(t_sec.astype("int64"))
    step_s = int(diffs[0])
    if step_s <= 0:
        raise ValueError("Non-positive timestep inferred from times_utc.")
    if not np.all(diffs == step_s):
        raise ValueError("times_utc must be uniformly sampled.")
    if 86400 % step_s != 0:
        raise ValueError(f"Inferred cadence {step_s}s does not divide 86400s.")

    slots_per_day = 86400 // step_s
    log.debug(f"Inferred cadence={step_s}s, slots_per_day={slots_per_day}")

    n_complete = (n // slots_per_day) * slots_per_day
    if n_complete == 0:
        raise ValueError("Not enough samples for a single complete day.")

    alt = np.asarray(alt_deg[:n_complete], dtype=float)
    az = np.mod(np.asarray(az_deg[:n_complete], dtype=float), 360.0)

    n_days = n_complete // slots_per_day
    log.debug(f"Using {n_days} full days of data ({n_complete} samples)")

    alt_2d = alt.reshape(n_days, slots_per_day)
    az_2d = az.reshape(n_days, slots_per_day)

    cols = np.arange(slots_per_day)
    idx_min = np.argmin(alt_2d, axis=0)
    idx_max = np.argmax(alt_2d, axis=0)

    slot_min_alt = alt_2d[idx_min, cols]
    slot_min_az = az_2d[idx_min, cols]
    slot_max_alt = alt_2d[idx_max, cols]
    slot_max_az = az_2d[idx_max, cols]

    if smooth_n < 4:
        raise ValueError("smooth_n must be >= 4 for cubic spline fitting.")

    az_smooth = np.linspace(0.0, 360.0, smooth_n, endpoint=False)

    def _periodic_cubic(az_in: np.ndarray, alt_in: np.ndarray) -> np.ndarray:
        order = np.argsort(az_in)
        x = az_in[order]
        y = alt_in[order]

        xu, idx = np.unique(x, return_index=True)
        yu = y[idx]
        if xu.size < 4:
            raise ValueError("Not enough unique azimuth samples for cubic spline.")

        xw = np.concatenate([xu, [xu[0] + 360.0]])
        yw = np.concatenate([yu, [yu[0]]])

        cs = CubicSpline(xw, yw, bc_type="periodic")
        return cs(az_smooth)

    min_alt_smooth = _periodic_cubic(slot_min_az, slot_min_alt)
    max_alt_smooth = _periodic_cubic(slot_max_az, slot_max_alt)

    az_plot = np.append(az_smooth, 360.0)
    alt_min_plot = np.append(min_alt_smooth, min_alt_smooth[0])
    alt_max_plot = np.append(max_alt_smooth, max_alt_smooth[0])

    log.info(f"Solar envelope built: smooth_n={smooth_n}, output grid length={az_plot.size}")
    return az_plot, alt_min_plot, alt_max_plot

solshade.solar¶

compute_solar_ephem ¶

load_sun_ephemeris ¶

nearest_horizon_indices ¶

solar_envelope_by_folding ¶

`solshade.solar`¶