Last active
May 8, 2020 09:15
-
-
Save chrisleaman/a6b189d3406a19454a873b6490a56439 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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