plot.png

reflector.py

import datetime
from dataclasses import dataclass

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pytz
import utm
from matplotlib import colors
from matplotlib.patches import Patch
from pysolar.solar import get_altitude_fast, get_azimuth_fast


@dataclass
class Point:
    """
    Point class which helps us convert between lat,lon and x,y easily.
    """

    lat: float
    lon: float
    z: float

    @property
    def x(self):
        x, _, _, _ = utm.from_latlon(self.lat, self.lon)
        return x

    @property
    def y(self):
        _, y, _, _ = utm.from_latlon(self.lat, self.lon)
        return y


class ReflectionCalculator:
    def __init__(self, mirror_pt1, mirror_pt2, focus_pt):

        self.mirror_pt1 = mirror_pt1
        self.mirror_pt2 = mirror_pt2
        self.focus_pt = focus_pt

        self.valid_sun_azimuths, self.valid_sun_elevations = self.find_sun_positions()

    @property
    def mirror_azimuth(self):
        """
        Returns the azimuth of the mirror in degrees. Needed to calculate the range
        of sun positions which result in a reflection.
        """

        return 90 - np.degrees(
            np.arctan(
                (self.mirror_pt2.y - self.mirror_pt1.y)
                / (self.mirror_pt2.x - self.mirror_pt1.x)
            )
        )

    def find_sun_positions(self):
        """
        Find the minimum and maximum azimuth and elevations that the sun must be
        positioned at to reflect off the building mirror and into the focus point
        """

        # Check both mirror points
        mirror_pts = [self.mirror_pt1, self.mirror_pt2]

        # Initialize lists to store the calculated azimuths and elevations
        sun_azimuths = []
        sun_elevations = []

        for mirror_pt in mirror_pts:

            elevation, azimuth = self.angles_between_points(mirror_pt, self.focus_pt)

            # Calculate the position of where the sun needs to be to result in a
            # reflection
            sun_azimuths.append(90 + azimuth + 2 * self.mirror_azimuth)
            sun_elevations.append(-elevation)

        valid_sun_azimuths = (min(sun_azimuths), max(sun_azimuths))
        valid_sun_elevations = (min(sun_elevations), max(sun_elevations))
        return valid_sun_azimuths, valid_sun_elevations

    def calculate_reflective_times(
        self, start, end, freq, tz, hour_start=14, hour_end=18
    ):
        """
        Creates a pandas dataframe containing whether the sun will be reflecting into
        the apartment at a each time or not.

        :param start: datetime of the start time to analyse
        :param end: datetime of the end time to analyse
        :param freq: str of the frequency of items between start and end
        :param tz: pytz.timezone of the timezone of the location
        :param hour_start: int of the minimum hour of day to check
        :param hour_end: int of the maxmimum hour of day to check
        """

        # Use position at the center of the mirror
        mirror_lat = np.mean([self.mirror_pt1.lat, self.mirror_pt2.lat])
        mirror_lon = np.mean([self.mirror_pt1.lon, self.mirror_pt2.lon])

        # Build timeseries
        index = pd.date_range(start, end, freq=freq, tz=tz)

        # Make dataframe, but only interested in certain times during the day
        df = pd.DataFrame(index=index)

        df = df.loc[
            (df.index.to_series().dt.hour >= hour_start)
            & (df.index.to_series().dt.hour < hour_end)
        ]

        # Calculate position of sun at each time step
        df["sun_altitudes"] = df.index.to_series().apply(
            lambda x: get_altitude_fast(mirror_lat, mirror_lon, x)
        )
        df["sun_azimuths"] = df.index.to_series().apply(
            lambda x: get_azimuth_fast(mirror_lat, mirror_lon, x)
        )

        # Determine whether sun is in correct position to reflect off mirror
        df["in_elevation"] = False
        df["in_azimuth"] = False
        df["will_reflect"] = False

        df.loc[
            df.sun_altitudes.between(*self.valid_sun_elevations), "in_elevation"
        ] = True
        df.loc[df.sun_azimuths.between(*self.valid_sun_azimuths), "in_azimuth"] = True
        df.loc[(df.in_elevation) & (df.in_azimuth), "will_reflect"] = True

        return df

    @staticmethod
    def angles_between_points(pt1, pt2):
        """
        Returns the elevation and azimuth in degrees at pt1, looking at pt2
        """

        dx = pt2.x - pt1.x
        dy = pt2.y - pt1.y
        dz = pt2.z - pt1.z

        azimuth = np.degrees(np.arctan2(dy, dx))

        if azimuth < 0:
            azimuth += 360

        elevation = np.degrees(np.arctan(dz / np.sqrt(dx ** 2 + dy ** 2)))

        return elevation, azimuth


def list_day_start_end_times(df, output_csv):
    """
    Summarizes a pandas dataframe of times where sun will be reflecting. Output
    shows each day, the start and end times of reflection as well as the duration of
    reflection.
    """

    # Add day to dataframe
    df["day"] = df.index.to_series().dt.floor("d").dt.date
    df["datetime"] = df.index.to_series()

    # Drop times where it's not reflecting
    df = df[df.will_reflect]

    def process_day(x):
        return pd.Series(
            {
                "reflection_start_time": x.datetime.min().time(),
                "reflection_end_time": x.datetime.max().time(),
                "reflection_mins": (x.datetime.max() - x.datetime.min()).seconds / 60,
            }
        )

    df_days = df.groupby("day").apply(process_day)
    df_days.to_csv(output_csv)


def plot_time_heatmap(df, output_file, use_tex=False):
    """
    Plot a time heat map of the times where the sun is reflecting
    """

    # Let's pivot the table so it's a bit more useful
    # Add day to dataframe
    df["day"] = df.index.to_series().dt.floor("d")
    df["time"] = df.index.to_series().dt.time

    # Create another column with the reason if it will reflect or not
    df["reflect_reason"] = 3
    df.loc[(df.in_azimuth) & (~df.in_elevation), "reflect_reason"] = 2
    df.loc[(~df.in_azimuth) & (df.in_elevation), "reflect_reason"] = 1
    df.loc[(df.in_azimuth) & (df.in_elevation), "reflect_reason"] = 0

    # Pivot the table, se we have time as the index and each day as a column with if
    # it will reflect as the values
    df = pd.pivot_table(df, columns="day", index="time", values="reflect_reason")

    # Data to plot.
    plot_data = df.values

    # Use custom color map
    cmap = colors.ListedColormap(["#9e0142", "#fee08b", "#abdda4", "#5e4fa2"])
    bounds = [-0.5, 0.5, 1.5, 2.5, 3.5]
    norm = colors.BoundaryNorm(bounds, cmap.N)

    # Create y_lims
    y_vals = [datetime.datetime.combine(datetime.date(2000, 1, 1), x) for x in df.index]
    y_vals = mdates.date2num(y_vals)

    # Create x-lims
    x_vals = mdates.date2num(df.columns)

    # Add settings to customize look of plot
    plt.rc("grid", linestyle="--", color="black", alpha=0.3)
    plt.rc("font", family="sans-serif")

    # Note that LaTeX needs to be install on the machine for this to run properly.
    # It's not necessary to use LaTeX, but you get nicer fonts
    if use_tex:
        plt.rc("text", usetex=True)
        plt.rc(
            "text.latex",
            preamble=[
                r"\usepackage{siunitx}",
                r"\sisetup{detect-all}",
                r"\usepackage[default]{sourcesanspro}",
                r"\usepackage{amsmath}",
                r"\usepackage{sansmath}",
                r"\sansmath",
            ],
        )

    # Create figure
    fig, ax = plt.subplots(figsize=(5, 4))

    _ = ax.imshow(
        plot_data,
        cmap=cmap,
        norm=norm,
        extent=[x_vals[0], x_vals[-1], y_vals[-1], y_vals[0]],
        origin="upper",
        aspect="auto",
    )

    ax.xaxis_date()
    ax.yaxis_date()

    month_year_format = mdates.DateFormatter("%b %Y")
    hour_min_format = mdates.DateFormatter("%H:%M %p")

    ax.xaxis.set_major_formatter(month_year_format)
    ax.yaxis.set_major_formatter(hour_min_format)

    # Rotate ticks
    plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")

    # Make custom legend
    legend_elements = [
        Patch(
            facecolor="#9e0142",
            edgecolor="#9e0142",
            label="Correct elevation\n& azimuth",
        ),
        Patch(facecolor="#fee08b", edgecolor="#fee08b", label="Correct elevation"),
        Patch(facecolor="#abdda4", edgecolor="#abdda4", label="Correct azimuth"),
        Patch(
            facecolor="#5e4fa2", edgecolor="#5e4fa2", label="Wrong elevation\n& azimuth"
        ),
    ]
    ax.legend(handles=legend_elements, prop={"size": 6})

    plt.grid(True)
    ax.set_title("When will sun reflect from building to home?")
    plt.tight_layout()

    fig.savefig(output_file, dpi=300)


if __name__ == "__main__":

    # Coordinates redacted here, but feel free to substitue your own.

    # Mirror points define the plane of the building which is reflecting the sun
    mirror_pt1 = Point(lat=-33.8, lon=151.2, z=200)
    mirror_pt2 = Point(lat=-33.8, lon=151.2, z=100)

    # Focus point is the apartment
    focus_pt = Point(lat=-33.8, lon=151.2, z=0)

    # Initalize our reflectoion calculator class
    reflector = ReflectionCalculator(mirror_pt1, mirror_pt2, focus_pt)

    # Create a dataframe for the period we're interested in which calculates whether
    # the sun will reflect into the apartment
    df = reflector.calculate_reflective_times(
        start=datetime.datetime(2020, 3, 1),
        end=datetime.datetime(2020, 10, 31),
        freq="1min",
        tz=pytz.timezone("Australia/Brisbane"),
    )

    print(df.head())

    # Generate a .csv with list of times for each day and plot a heatmap
    list_day_start_end_times(df, output_csv="day_start_end_times.csv")
    plot_time_heatmap(df, output_file="plot.png")