From 411198c9b02b50a8b8391ffe78d0dda4ee4f9181 Mon Sep 17 00:00:00 2001 From: Roeland Trommel Date: Mon, 25 Nov 2024 22:29:24 +0100 Subject: [PATCH 1/3] initial rewrite --- src/shadowfinder/shadowfinder.py | 161 +++++++++++++++++++++++-------- 1 file changed, 123 insertions(+), 38 deletions(-) diff --git a/src/shadowfinder/shadowfinder.py b/src/shadowfinder/shadowfinder.py index 63f67ce..f38b72a 100644 --- a/src/shadowfinder/shadowfinder.py +++ b/src/shadowfinder/shadowfinder.py @@ -1,9 +1,11 @@ from pytz import timezone, utc +import datetime import pandas as pd from suncalc import get_position import numpy as np import matplotlib.pyplot as plt import matplotlib.colors as colors +from matplotlib.figure import Figure from mpl_toolkits.basemap import Basemap from timezonefinder import TimezoneFinder import json @@ -39,15 +41,15 @@ def __init__( self.min_lon = -180 self.max_lon = 180 - def set_details( - self, - date_time, - object_height=None, - shadow_length=None, - time_format=None, - sun_altitude_angle=None, - ): - + def set_datetime( + self, date_time: datetime.datetime = None, time_format=None + ) -> None: + if time_format is not None: + assert time_format in [ + "utc", + "local", + ], "time_format must be 'utc' or 'local'" + self.time_format = time_format if date_time is not None and date_time.tzinfo is not None: warn( "date_time is expected to be timezone naive (i.e. tzinfo=None). Any timezone information will be ignored." @@ -55,12 +57,15 @@ def set_details( date_time = date_time.replace(tzinfo=None) self.date_time = date_time - if time_format is not None: - assert time_format in [ - "utc", - "local", - ], "time_format must be 'utc' or 'local'" - self.time_format = time_format + def set_details( + self, + date_time, + object_height=None, + shadow_length=None, + time_format=None, + sun_altitude_angle=None, + ) -> None: + self.set_datetime(date_time=date_time, time_format=time_format) # height and length must have the same None-ness # either height or angle must be set (but not both or neither) @@ -93,14 +98,14 @@ def set_details( # Lengths and angle are None and we use the same values as before pass - def quick_find(self, timezone_grid="timezone_grid.json"): + def quick_find(self, timezone_grid="timezone_grid.json") -> None: # try to load timezone grid from file, generate if not found try: self.load_timezone_grid(timezone_grid) except FileNotFoundError: self.generate_timezone_grid() - self.find_shadows() + _ = self.find_shadows() fig = self.plot_shadows() if self.sun_altitude_angle is not None: @@ -110,7 +115,7 @@ def quick_find(self, timezone_grid="timezone_grid.json"): fig.savefig(file_name) - def generate_timezone_grid(self): + def generate_timezone_grid(self) -> None: lats = np.arange(self.min_lat, self.max_lat, self.angular_resolution) lons = np.arange(self.min_lon, self.max_lon, self.angular_resolution) @@ -124,7 +129,7 @@ def generate_timezone_grid(self): ] ) - def save_timezone_grid(self, filename="timezone_grid.json"): + def save_timezone_grid(self, filename="timezone_grid.json") -> None: data = { "min_lat": self.min_lat, "max_lat": self.max_lat, @@ -136,7 +141,7 @@ def save_timezone_grid(self, filename="timezone_grid.json"): json.dump(data, open(filename, "w")) - def load_timezone_grid(self, filename="timezone_grid.json"): + def load_timezone_grid(self, filename="timezone_grid.json") -> None: data = json.load(open(filename, "r")) self.min_lat = data["min_lat"] @@ -151,14 +156,26 @@ def load_timezone_grid(self, filename="timezone_grid.json"): self.lons, self.lats = np.meshgrid(lons, lats) self.timezones = np.array(data["timezones"]) - def find_shadows(self): + def find_shadows( + self, + object_height: float = None, + shadow_length: float = None, + sun_altitude_angle: float = None, + timestamp: datetime.datetime = None, + ) -> np.ndarray: # Evaluate the sun's length at a grid of points on the Earth's surface + object_height = object_height if object_height else self.object_height + shadow_length = shadow_length if shadow_length else self.shadow_length + sun_altitude_angle = ( + sun_altitude_angle if sun_altitude_angle else self.sun_altitude_angle + ) + timestamp = timestamp if timestamp else self.date_time if self.lats is None or self.lons is None or self.timezones is None: self.generate_timezone_grid() if self.time_format == "utc": - valid_datetimes = utc.localize(self.date_time) + valid_datetimes = utc.localize(timestamp) valid_lats = self.lats.flatten() valid_lons = self.lons.flatten() elif self.time_format == "local": @@ -168,7 +185,7 @@ def find_shadows(self): None if tz is None else timezone(tz) - .localize(self.date_time) + .localize(timestamp) .astimezone(utc) .timestamp() ) @@ -194,9 +211,9 @@ def find_shadows(self): # If object height and shadow length are set the sun altitudes are used # to calculate the shadow lengths across the world and then compared to # the expected shadow length. - if self.object_height is not None and self.shadow_length is not None: + if object_height is not None and shadow_length is not None: # Calculate the shadow length - shadow_lengths = self.object_height / np.apply_along_axis( + shadow_lengths = object_height / np.apply_along_axis( np.tan, 0, valid_sun_altitudes ) @@ -204,17 +221,15 @@ def find_shadows(self): shadow_lengths[valid_sun_altitudes <= 0] = np.nan # Show the relative difference between the calculated shadow length and the observed shadow length - location_likelihoods = ( - shadow_lengths - self.shadow_length - ) / self.shadow_length + location_likelihoods = (shadow_lengths - shadow_length) / shadow_length # If the sun altitude angle is set then this value is directly compared # to the sun altitudes across the world. - elif self.sun_altitude_angle is not None: + elif sun_altitude_angle is not None: # Show relative difference between sun altitudes location_likelihoods = ( - np.array(valid_sun_altitudes) - radians(self.sun_altitude_angle) - ) / radians(self.sun_altitude_angle) + np.array(valid_sun_altitudes) - radians(sun_altitude_angle) + ) / radians(sun_altitude_angle) # Replace points where the sun is below the horizon location_likelihoods[valid_sun_altitudes <= 0] = np.nan @@ -237,12 +252,29 @@ def find_shadows(self): self.location_likelihoods = np.reshape( self.location_likelihoods, np.shape(self.lons), order="A" ) + return self.location_likelihoods def plot_shadows( self, figure_args={"figsize": (12, 6)}, basemap_args={"projection": "cyl", "resolution": "c"}, - ): + object_height: float = None, + shadow_length: float = None, + sun_altitude_angle: float = None, + timestamp: datetime.datetime = None, + location_likelihoods: np.ndarray = None, + plt_title: str = None, + ) -> Figure: + object_height = object_height if object_height else self.object_height + shadow_length = shadow_length if shadow_length else self.shadow_length + sun_altitude_angle = ( + sun_altitude_angle if sun_altitude_angle else self.sun_altitude_angle + ) + timestamp = timestamp if timestamp else self.date_time + location_likelihoods = ( + location_likelihoods if location_likelihoods else self.location_likelihoods + ) + plt_title = plt_title if plt_title else None fig = plt.figure(**figure_args) @@ -262,19 +294,72 @@ def plot_shadows( m.pcolormesh( x, y, - np.abs(self.location_likelihoods), + np.abs(location_likelihoods), cmap=cmap, norm=norm, alpha=0.7, ) # plt.colorbar(label='Relative Shadow Length Difference') - - if self.sun_altitude_angle is not None: - plt_title = f"Possible Locations at {self.date_time.strftime('%Y-%m-%d %H:%M:%S')} {self.time_format.title()}\n(sun altitude angle: {self.sun_altitude_angle})" - else: - plt_title = f"Possible Locations at {self.date_time.strftime('%Y-%m-%d %H:%M:%S')} {self.time_format.title()}\n(object height: {self.object_height}, shadow length: {self.shadow_length})" + if plt_title is None: + if sun_altitude_angle is not None: + plt_title = f"Possible Locations at {timestamp.strftime('%Y-%m-%d %H:%M:%S')} {self.time_format.title()}\n(sun altitude angle: {sun_altitude_angle})" + else: + plt_title = f"Possible Locations at {timestamp.strftime('%Y-%m-%d %H:%M:%S')} {self.time_format.title()}\n(object height: {object_height}, shadow length: {shadow_length})" plt.title(plt_title) self.fig = fig return fig + + def find_multiple_shadows( + self, + object_heights: list[float], + shadow_lengths: list[float], + timestamps: list[datetime.datetime], + ) -> Figure: + """ + The plot title string assumes the timestamps are ordered ascending; + """ + if len(shadow_lengths) != len(timestamps): + raise ValueError( + f"Argument lists are expected to have the same length. Provided {len(shadow_lengths)},{len(timestamps)}." + ) + elif len(object_heights) != len(timestamps): + if len(object_heights) == 1: + object_heights = [object_heights[0] for _ in timestamps] + else: + raise ValueError( + f"object_heights must be a list of the same length as the shadow lengths and timestamps, or contain only 1 value" + ) + + list_location_likelihoods = [ + self.find_shadows( + object_height=obj_height, + shadow_length=shadow_length, + timestamp=timestamp, + ) + for obj_height, shadow_length, timestamp in zip( + object_heights, shadow_lengths, timestamps + ) + ] + + # merge location_likelihoods: strict "AND" operation? + location_likelihoods = list_location_likelihoods[0] + max_val_for_rescaling = np.amax(location_likelihoods) + for loc_likelihoods in list_location_likelihoods: + max_val_for_rescaling = max(np.amax(loc_likelihoods), max_val_for_rescaling) + location_likelihoods = np.multiply(location_likelihoods, loc_likelihoods) + + # rescale for plotting + cur_max = np.amax(location_likelihoods) + if cur_max > 0: + location_likelihoods = location_likelihoods * ( + max_val_for_rescaling / cur_max + ) + + plt_title = f"Possible Locations for measurements between {timestamps[0].strftime('%Y-%m-%d %H:%M:%S')} and {timestamps[-1].strftime('%Y-%m-%d %H:%M:%S')} {self.time_format.title()}\n)" + fig = self.plot_shadows( + location_likelihoods=location_likelihoods, plt_title=plt_title + ) + + return fig From 0a1c5ae2628234e3a7c1bb6a2924c27c6ecf0adb Mon Sep 17 00:00:00 2001 From: Roeland Trommel Date: Mon, 25 Nov 2024 22:43:41 +0100 Subject: [PATCH 2/3] add test for multiple shadow find --- src/shadowfinder/shadowfinder.py | 8 +++++++- tests/test_shadowfinder.py | 16 +++++++++++++++- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/shadowfinder/shadowfinder.py b/src/shadowfinder/shadowfinder.py index f38b72a..9e8727c 100644 --- a/src/shadowfinder/shadowfinder.py +++ b/src/shadowfinder/shadowfinder.py @@ -22,6 +22,10 @@ def __init__( time_format="utc", sun_altitude_angle=None, ): + self.sun_altitude_angle = None + self.object_height = None + self.shadow_length = None + self.set_details( date_time, object_height, shadow_length, time_format, sun_altitude_angle ) @@ -272,7 +276,9 @@ def plot_shadows( ) timestamp = timestamp if timestamp else self.date_time location_likelihoods = ( - location_likelihoods if location_likelihoods else self.location_likelihoods + location_likelihoods + if isinstance(location_likelihoods, np.ndarray) + else self.location_likelihoods ) plt_title = plt_title if plt_title else None diff --git a/tests/test_shadowfinder.py b/tests/test_shadowfinder.py index 531c61e..f821683 100644 --- a/tests/test_shadowfinder.py +++ b/tests/test_shadowfinder.py @@ -1,4 +1,5 @@ -from datetime import datetime +from datetime import datetime, timedelta + from shadowfinder import ShadowFinder @@ -15,3 +16,16 @@ def test_creation_with_valid_arguments_should_pass(): ShadowFinder( object_height=object_height, shadow_length=shadow_length, date_time=date_time ) + +def test_find_multiple_shadows(): + """Baseline test to assert that we can create an instance of ShadowFinder with only object height, shadow length, + and a datetime object.""" + # GIVEN + object_heights = [6,6,6,6] + shadow_lengths = [3.2,3.0,2.9,2.8] + time_offsets = range(4) + timestamps = [datetime.now()+timedelta(hours=f) for f in time_offsets] + + shadow_finder = ShadowFinder() + # WHEN / THEN + figure = shadow_finder.find_multiple_shadows(object_heights=object_heights, shadow_lengths=shadow_lengths, timestamps=timestamps) From 7b29326ea786a2fb2324d4543e3c99ef0da1b9a1 Mon Sep 17 00:00:00 2001 From: Roeland Trommel Date: Wed, 27 Nov 2024 12:23:24 +0100 Subject: [PATCH 3/3] alpha version find multiple shadows --- src/shadowfinder/shadowfinder.py | 141 ++++++++++++++++++------------- tests/test_shadowfinder.py | 3 +- 2 files changed, 82 insertions(+), 62 deletions(-) diff --git a/src/shadowfinder/shadowfinder.py b/src/shadowfinder/shadowfinder.py index 9e8727c..ad7637b 100644 --- a/src/shadowfinder/shadowfinder.py +++ b/src/shadowfinder/shadowfinder.py @@ -160,24 +160,9 @@ def load_timezone_grid(self, filename="timezone_grid.json") -> None: self.lons, self.lats = np.meshgrid(lons, lats) self.timezones = np.array(data["timezones"]) - def find_shadows( - self, - object_height: float = None, - shadow_length: float = None, - sun_altitude_angle: float = None, - timestamp: datetime.datetime = None, - ) -> np.ndarray: - # Evaluate the sun's length at a grid of points on the Earth's surface - - object_height = object_height if object_height else self.object_height - shadow_length = shadow_length if shadow_length else self.shadow_length - sun_altitude_angle = ( - sun_altitude_angle if sun_altitude_angle else self.sun_altitude_angle - ) - timestamp = timestamp if timestamp else self.date_time - if self.lats is None or self.lons is None or self.timezones is None: - self.generate_timezone_grid() + def generate_valid_datetimes_lats_lons_and_mask(self, timestamp): + mask = None if self.time_format == "utc": valid_datetimes = utc.localize(timestamp) valid_lats = self.lats.flatten() @@ -208,6 +193,29 @@ def find_shadows( # Convert the datetimes to pandas series of timestamps valid_datetimes = pd.to_datetime(valid_datetimes, unit="s", utc=True) + return valid_datetimes, valid_lons, valid_lats, mask + + def find_shadows( + self, + object_height: float = None, + shadow_length: float = None, + sun_altitude_angle: float = None, + timestamp: datetime.datetime = None, + ) -> np.ndarray: + # Evaluate the sun's length at a grid of points on the Earth's surface + + object_height = object_height if object_height else self.object_height + shadow_length = shadow_length if shadow_length else self.shadow_length + sun_altitude_angle = ( + sun_altitude_angle if sun_altitude_angle else self.sun_altitude_angle + ) + timestamp = timestamp if timestamp else self.date_time + if self.lats is None or self.lons is None or self.timezones is None: + self.generate_timezone_grid() + + valid_datetimes, valid_lons, valid_lats, mask = ( + self.generate_valid_datetimes_lats_lons_and_mask(timestamp=timestamp) + ) pos_obj = get_position(valid_datetimes, valid_lons, valid_lats) valid_sun_altitudes = pos_obj["altitude"] # in radians @@ -258,29 +266,39 @@ def find_shadows( ) return self.location_likelihoods + def get_plot_data(self, data: dict): + object_height = data.get("object_height", self.object_height) + shadow_length = data.get("shadow_length", self.shadow_length) + sun_altitude_angle = data.get("sun_altitude_angle", self.sun_altitude_angle) + timestamp = data.get("timestamp", self.date_time) + location_likelihoods = data.get( + "location_likelihoods", self.location_likelihoods + ) + plt_title = data.get("plt_title", None) + return ( + object_height, + shadow_length, + sun_altitude_angle, + timestamp, + location_likelihoods, + plt_title, + ) + def plot_shadows( self, figure_args={"figsize": (12, 6)}, basemap_args={"projection": "cyl", "resolution": "c"}, - object_height: float = None, - shadow_length: float = None, - sun_altitude_angle: float = None, - timestamp: datetime.datetime = None, - location_likelihoods: np.ndarray = None, - plt_title: str = None, + plot_data: dict = None, ) -> Figure: - object_height = object_height if object_height else self.object_height - shadow_length = shadow_length if shadow_length else self.shadow_length - sun_altitude_angle = ( - sun_altitude_angle if sun_altitude_angle else self.sun_altitude_angle - ) - timestamp = timestamp if timestamp else self.date_time - location_likelihoods = ( - location_likelihoods - if isinstance(location_likelihoods, np.ndarray) - else self.location_likelihoods - ) - plt_title = plt_title if plt_title else None + plot_data = plot_data if plot_data else {} + ( + object_height, + shadow_length, + sun_altitude_angle, + timestamp, + location_likelihoods, + plt_title, + ) = self.get_plot_data(data=plot_data) fig = plt.figure(**figure_args) @@ -324,7 +342,9 @@ def find_multiple_shadows( timestamps: list[datetime.datetime], ) -> Figure: """ - The plot title string assumes the timestamps are ordered ascending; + Find possible locations using multiple shadows as input. + + If the same object is used every time, the object height can be provided as a list with a single value """ if len(shadow_lengths) != len(timestamps): raise ValueError( @@ -337,35 +357,36 @@ def find_multiple_shadows( raise ValueError( f"object_heights must be a list of the same length as the shadow lengths and timestamps, or contain only 1 value" ) + if self.lats is None or self.lons is None or self.timezones is None: + self.generate_timezone_grid() - list_location_likelihoods = [ - self.find_shadows( + arr_location_likelihoods = np.zeros( + (len(timestamps), self.lons.shape[0], self.lons.shape[1]) + ) + for k, (obj_height, shadow_length, timestamp) in enumerate( + zip(object_heights, shadow_lengths, timestamps) + ): + arr_location_likelihoods[k, :, :] = self.find_shadows( object_height=obj_height, shadow_length=shadow_length, timestamp=timestamp, ) - for obj_height, shadow_length, timestamp in zip( - object_heights, shadow_lengths, timestamps - ) - ] - - # merge location_likelihoods: strict "AND" operation? - location_likelihoods = list_location_likelihoods[0] - max_val_for_rescaling = np.amax(location_likelihoods) - for loc_likelihoods in list_location_likelihoods: - max_val_for_rescaling = max(np.amax(loc_likelihoods), max_val_for_rescaling) - location_likelihoods = np.multiply(location_likelihoods, loc_likelihoods) - - # rescale for plotting - cur_max = np.amax(location_likelihoods) - if cur_max > 0: - location_likelihoods = location_likelihoods * ( - max_val_for_rescaling / cur_max - ) - plt_title = f"Possible Locations for measurements between {timestamps[0].strftime('%Y-%m-%d %H:%M:%S')} and {timestamps[-1].strftime('%Y-%m-%d %H:%M:%S')} {self.time_format.title()}\n)" - fig = self.plot_shadows( - location_likelihoods=location_likelihoods, plt_title=plt_title - ) + # merge/logic for combining the shadows. #TODO + location_likelihoods = np.max(arr_location_likelihoods, axis=0) + + plt_title = f"Possible Locations for measurements at [" + for ts in timestamps: + plt_title += f"{ts.strftime('%Y-%m-%d %H:%M:%S')}, " + plt_title = plt_title[:-2] + f"] ({self.time_format.title()}\n)" + + plot_data = { + "location_likelihoods": location_likelihoods, + "plt_title": plt_title, + "timestamps": timestamps, + "object_heights": object_heights, + "shadow_lengths": shadow_lengths, + } + fig = self.plot_shadows(plot_data=plot_data) return fig diff --git a/tests/test_shadowfinder.py b/tests/test_shadowfinder.py index f821683..3460d28 100644 --- a/tests/test_shadowfinder.py +++ b/tests/test_shadowfinder.py @@ -18,8 +18,7 @@ def test_creation_with_valid_arguments_should_pass(): ) def test_find_multiple_shadows(): - """Baseline test to assert that we can create an instance of ShadowFinder with only object height, shadow length, - and a datetime object.""" + """Test find_multiple_shadows""" # GIVEN object_heights = [6,6,6,6] shadow_lengths = [3.2,3.0,2.9,2.8]