Skip to content

Instantly share code, notes, and snippets.

@chrisleaman
Last active May 8, 2020 09:15
Show Gist options
  • Save chrisleaman/a6b189d3406a19454a873b6490a56439 to your computer and use it in GitHub Desktop.
Save chrisleaman/a6b189d3406a19454a873b6490a56439 to your computer and use it in GitHub Desktop.
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment