Add solar intensity based on declination angle (single parameter "latitude")

openskistats/analyze.py

@@ -13,7 +13,9 @@
     add_spatial_metric_columns,
     get_bearing_histograms,
     get_bearing_summary_stats,
+    get_solar_intensity_batch,
 )
+
 from openskistats.models import BearingStatsModel, SkiAreaModel
 from openskistats.openskimap_utils import (
     load_downhill_ski_areas_from_download_pl,
@@ -93,6 +95,11 @@ def analyze_all_ski_areas_polars(skip_runs: bool = False) -> None:
             combined_vertical=pl.col("distance_vertical_drop").sum(),
             combined_distance=pl.col("distance_3d").sum(),
             latitude=pl.col("latitude").mean(),
+            solar_intensity = pl.col("latitude")
+                .map_batches(
+                    lambda batch: get_solar_intensity_batch(batch),
+                    return_dtype=pl.Float64
+            ).mean(),
             longitude=pl.col("longitude").mean(),
             min_elevation=pl.col("elevation").min(),
             max_elevation=pl.col("elevation").max(),
@@ -221,6 +228,7 @@ def aggregate_ski_areas_pl(
             max_elevation=pl.max("max_elevation"),
             vertical_drop=pl.max("max_elevation") - pl.min("min_elevation"),
             latitude=pl.mean("latitude"),
+            solar_intensity=pl.mean("solar_intensity"),
             longitude=pl.mean("longitude"),
             _bearing_stats=pl.struct(
                 pl.col("bearing_mean").alias("bearing"),

openskistats/bearing.py

@@ -315,3 +315,100 @@ def get_bearing_summary_stats(
         bearing_magnitude_net=round(net_magnitude, 7),
         bearing_magnitude_cum=round(cum_magnitude, 7),
     )
+
+
+# def get_solar_intensity(
+#     latitudes: list[float] | npt.NDArray[np.float64]
+# ) -> pl.DataFrame:
+
+
+HOURS_IN_DAY = 24
+HOURS_START = 9  # Start time (9 AM)
+HOURS_END = 16   # End time (4 PM)
+
+def solar_declination(day_of_year: int) -> float:
+    """
+    Calculate the solar declination for a given day of the year.
+    This uses the approximation formula for solar declination.
+    """
+    # Calculate the solar declination using the approximation formula
+    # δ = 23.44 * sin(360° * (N + 10) / 365)
+    # Where N is the day of the year (1 to 365)
+    return 23.44 * np.sin(np.radians(360 * (day_of_year + 10) / 365))
+
+def solar_intensity_at_time(latitude: float, declination: float, hour: int) -> float:
+    """
+    Calculate the solar intensity at a given time of day (hour) for a given latitude and solar declination.
+    """
+    # Convert latitude and declination to radians for trigonometric calculations
+    latitude_rad = np.radians(latitude)
+    declination_rad = np.radians(declination)
+
+    # Hour angle (ω) for a given time of day (simple model)
+    # Hour angle increases by 15 degrees per hour, starting from 12:00 PM.
+    hour_angle = np.radians(15 * (hour - 12))
+
+    # Calculate the solar zenith angle using the formula
+    # cos(θ) = sin(δ) * sin(φ) + cos(δ) * cos(φ) * cos(ω)
+    # Where:
+    # δ = solar declination, φ = latitude, ω = hour angle
+    solar_zenith = np.arccos(
+        np.sin(declination_rad) * np.sin(latitude_rad) +
+        np.cos(declination_rad) * np.cos(latitude_rad) * np.cos(hour_angle)
+    )
+
+    intensity = np.cos(solar_zenith)
+    return max(intensity, 0)
+
+def integrate_solar_intensity(latitude: float, day_of_year: int) -> float:
+    """
+    Integrate solar intensity from 9 AM to 4 PM.
+    The function returns the average solar intensity over this period.
+    """
+    declination = solar_declination(day_of_year)
+    intensities = [solar_intensity_at_time(latitude, declination, hour) for hour in range(HOURS_START, HOURS_END)]
+
+    return np.mean(intensities)
+
+def get_hemisphere(latitude: float) -> Literal["north", "south"] | None:
+    if latitude is None:
+        return None
+    elif latitude >= 0:
+        return "north"
+    elif latitude < 0:
+        return "south"
+
+
+def get_solar_intensity(
+    latitude: float
+) -> float:
+    """
+    Calculate the average solar intensity from 9 AM to 4 PM for each latitude.
+
+    Args:
+        latitudes (list or ndarray): A list or numpy array of latitudes.
+        hemisphere: lets us know if we are in the northern or southern hemisphere
+
+    Returns:
+        polars.DataFrame: A DataFrame with latitudes and their corresponding average solar intensity.
+    """
+    # For each latitude, calculate the solar intensity (assuming a fixed day of the year, e.g., 172nd day for mid-winter in southern hemisphere)
+    # and the 355th day for mid-winter in northern hemisphere
+    import logging
+    logging.info(f"latitude is {latitude}")
+    hemisphere: Literal["north", "south"] | None = None
+    hemisphere = get_hemisphere(latitude)
+    if hemisphere == "north":
+        day_of_year = 355 # Mid-winter day (e.g., December 21)
+    elif hemisphere == "south":
+        day_of_year = 172 # Mid-winter day (e.g., June 21)
+    else:
+        logging.info(f"hemisphere has value: {hemisphere}")
+        return None
+
+    # TODO we could dynamically use the open and close dates of each resort
+    avg_intensity = integrate_solar_intensity(latitude, day_of_year)
+    return avg_intensity
+
+def get_solar_intensity_batch(latitudes: pl.Series) -> pl.Series:
+    return pl.Series([get_solar_intensity(latitude) for latitude in latitudes])

openskistats/display.py

@@ -100,6 +100,7 @@ def get_ski_area_frontend_table(story: bool = False) -> pl.DataFrame:
             "country_code",
             "region",
             "locality",
+            "solar_intensity",
             "latitude",
             "run_count",
             "lift_count",
@@ -204,6 +205,7 @@ def _percent_diverging_style(ci: reactable.CellInfo) -> dict[str, Any] | None:
     # add custom descriptions including overwriting SkiAreaModel descriptions
     "latitude": "Hemisphere where the ski area is located, either ℕ for north or 𝕊 for south, "
     "as well as the latitude (φ) in decimal degrees.",
+    "solar_intensity"
     "bin_proportion_4_N": "Proportion of vertical-weighted run segments oriented with a northern cardinal direction.",
     "bin_proportion_4_E": "Proportion of vertical-weighted run segments oriented with an eastern cardinal direction.",
     "bin_proportion_4_S": "Proportion of vertical-weighted run segments oriented with a southern cardinal direction.",
@@ -385,6 +387,11 @@ def _ski_area_metric_style(ci: reactable.CellInfo) -> dict[str, Any] | None:
                 filter_method=reactable.JS("filterLatitude"),
                 min_width=60,
             ),
+            reactable.Column(
+                id="solar_intensity",
+                name=f"Solar",
+                **_column_kwargs_location_str,
+            ),
             reactable.Column(
                 id="run_count",
                 name="Runs",
@@ -511,7 +518,7 @@ def _ski_area_metric_style(ci: reactable.CellInfo) -> dict[str, Any] | None:
                 cell=reactable.JS("cellRose"),
                 # max_width=45,
                 class_="border-left",
-            ),
+            )
         ],
         column_groups=[
             reactable.ColGroup(
@@ -521,6 +528,7 @@ def _ski_area_metric_style(ci: reactable.CellInfo) -> dict[str, Any] | None:
                 name="Location",
                 columns=[
                     "latitude",
+                    "solar_intensity",
                     "country_emoji",
                     "country",
                     "region",

openskistats/models.py

@@ -252,6 +252,11 @@ class SkiAreaModel(Model):  # type: ignore [misc]
     hemisphere: Literal["north", "south"] | None = Field(
         description="Hemisphere of the ski area.",
     )
+    solar_intensity: float | None = Field(
+        description="Solar intensity of the ski resort calculated by the declination at that latitude at the solstice",
+        ge=0.0,
+        le=1.0,
+    )
     min_elevation: Annotated[
         float | None,
         Field(