ShadePath.py

# Calculate the position where the line pointing to the sun reaches a given altitude above the Earth's surface for a given target location & datetime (in UTC)
# Target location supplied in latitude and longitude coordinates in decimal degrees format.

# Reference: Trujillo, J. H. S. (2014). Calculation of the shadow-penumbra relation and its application on efficient architectural design. Solar energy, 110, 139–150.

# ToDo: 
# - calculate velocities
# - Filter to max velocity, and insert return trajectory
# - Add velocity vectors
# - Add insolation rate to data and map labels
# - Make width of path proportional to insolation rate?

import sys
import csv
import numpy as np
import datetime
from pysolar.solar import * # https://pysolar.readthedocs.io/en/latest/
import great_circle_calculator.great_circle_calculator as gcc
import pandas as pd

import plotly.express as px

# set up matplotlib
#import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

from enum import Enum

class ModeOptions(Enum):
	MOBILE = 'mobile',
	STATIC = 'static'

# class TestClass:    
#     def __init__(self, value: ValueOptions) -> None:
#         self.value = value



class ShadePath:

    # declare types, to support checking with mypy etc https://www.kite.com/blog/python/type-hinting/
    latitude_degrees: float
    longitude_degrees: float
    altitude: float
    results: dict

    mode: ModeOptions

    def __init__(self, latitude_degrees: float, longitude_degrees: float, altitude: float , mode: ModeOptions):
        self.type = type        
        self.latitude_degrees = latitude_degrees
        self.longitude_degrees = longitude_degrees
        self.altitude = altitude
        self.mode = mode

        self.results = dict()
        self.results['datetime'] = []
        self.results['local-time'] = []
        self.results['aboveHorizon'] = []
        self.results['solarAltitude'] = []
        self.results['azimuth'] = []
        self.results['insolation'] = []
        self.results['chord'] = []
        self.results['distance'] = []
        self.results['shading-longitude'] = [] 
        self.results['shading-latitude'] = []
        self.results['hour-of-day'] = []
        self.results['circuit'] = []

    # Calculate the position where the line reaches a given altitude
    def calculate_position(self, current_datetime):

        self.results['datetime'].append(current_datetime)
        self.results['local-time'].append(current_datetime.strftime('%H:%M:%S'))

        # Get position of the sun https://pysolar.readthedocs.io/en/latest/
        solar_altitude_degrees = get_altitude(self.latitude_degrees, self.longitude_degrees, current_datetime)  
        azimuth = get_azimuth(self.latitude_degrees, self.longitude_degrees, current_datetime)
        insolation = radiation.get_radiation_direct(current_datetime, solar_altitude_degrees) # https://pysolar.readthedocs.io/en/latest/#estimate-of-clear-sky-radiation

        # If sun above the horizon, calculate shading position
        if (solar_altitude_degrees > 0): # Sun above the horizon - https://rhodesmill.org/pyephem/quick.html#body-compute-target

            distance, chord_length = distanceFromAltitude(self.altitude, solar_altitude_degrees)

            if (self.mode == 'mobile'):
                shading_longitude, shading_latitude = locatePointB(self.longitude_degrees,self.latitude_degrees,azimuth,distance)
            elif (self.mode == 'static'): # bearing is 180 degrees opposite
                shading_longitude, shading_latitude = locatePointB(self.longitude_degrees,self.latitude_degrees,(azimuth - 180),distance)
            else:
                exit("ShadePath ERROR: Supplied mode '" + self.mode + "'not recognised")


            # append results
            self.results['aboveHorizon'].append(True)
            self.results['solarAltitude'].append(solar_altitude_degrees)
            self.results['azimuth'].append(azimuth)
            self.results['insolation'].append(insolation)
            self.results['chord'].append(chord_length)
            self.results['distance'].append(distance)
            self.results['shading-longitude'].append(shading_longitude) 
            self.results['shading-latitude'].append(shading_latitude)
            #self.results['checkDistance'].append( gcc.distance_between_points(target_point, position, unit='kilometers', haversine=True) )

        else:
            # append NULL results
            self.results['aboveHorizon'].append(False)
            self.results['solarAltitude'].append(0)
            self.results['azimuth'].append(0)
            self.results['insolation'].append(0)
            self.results['chord'].append(0)
            self.results['distance'].append(0)
            self.results['shading-longitude'].append(0)
            self.results['shading-latitude'].append(0)

            return True


    def dayCircuit(self,start_datetime,interval,circuit):
        # loop over one day + one interval (which completes the loop)

        current_datetime = start_datetime
        end_datetime = start_datetime + datetime.timedelta(days = 1, minutes = interval)
        total_minutes = 0

        while current_datetime < end_datetime:
            self.calculate_position(current_datetime)        
            self.results['circuit'].append(circuit)
            self.results['hour-of-day'].append(total_minutes/60)
            total_minutes += interval
            current_datetime += datetime.timedelta(minutes = interval)

        self.result_df = pd.DataFrame.from_dict(self.results)

        #return self.results
    
    def runCircuits(self,dates,timezone_string,interval):
        date_format = '%Y-%m-%d %H:%M:%S%z'
        circuit = 0
        for date in dates:
            date_string = date + ' 00:00:00' + timezone_string
            start_datetime = datetime.datetime.strptime(date_string, date_format)

            self.dayCircuit(start_datetime,interval,circuit)
            circuit += 1
            # print(result, file=sys.stderr)


    def interactiveMap(self,df):
        # Interactive map with plotly
        # see 
        # - https://plotly.com/python/lines-on-maps/
        # - https://plotly.com/python-api-reference/generated/plotly.express.line_geo.html
        # - https://plotly.com/python/map-configuration/

        #projection (str) – One of 'equirectangular', 'mercator', 'orthographic', 'natural earth', 'kavrayskiy7', 'miller', 'robinson', 'eckert4', 'azimuthal equal area', 'azimuthal equidistant', 'conic equal area', 'conic conformal', 'conic equidistant', 'gnomonic', 'stereographic', 'mollweide', 'hammer', 'transverse mercator', 'albers usa', 'winkel tripel', 'aitoff', or 'sinusoidal'
        #https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection
        center = {}
        center['lat'] = self.latitude_degrees
        center['lon'] = self.longitude_degrees
        fig = px.line_geo(data_frame = df , projection='azimuthal equidistant' ,lat= 'shading-latitude', lon='shading-longitude', hover_name='local-time', fitbounds = 'locations', center = center)
        fig.show()


    def staticMap(self,line_colours,dot_colour):
        # Static map with matplotlib
        # See: 
        # - https://matplotlib.org/basemap/api/basemap_api.html
        # - https://matplotlib.org/basemap/users/examples.html
        # - https://basemaptutorial.readthedocs.io/en/latest/subplots.html

        # width	- width of desired map domain in projection coordinates (meters).
        # height - height of desired map domain in projection coordinates (meters).
        # lon_0	- center of desired map domain (in degrees).
        # lat_0	- center of desired map domain (in degrees).

        # llcrnrlon	- longitude of lower left hand corner of the desired map domain (degrees).
        # llcrnrlat	- latitude of lower left hand corner of the desired map domain (degrees).
        # urcrnrlon	- longitude of upper right hand corner of the desired map domain (degrees).
        # urcrnrlat	- latitude of upper right hand corner of the desired map domain (degrees).

        df = self.result_df[self.result_df['aboveHorizon']] # filter to when sun above the horizon
   
        # ToDo: Need to limit to reasonable latitude range, otherwise breaks with low solar angle
        #df = df[df['distance'] < 200] # filter to max distance from target - https://www.geeksforgeeks.org/ways-to-filter-pandas-dataframe-by-column-values/

        longs = df['shading-longitude']
        lats = df['shading-latitude']

        # define the size of the map box
        margin_degrees = 0.5
        scale_offset = margin_degrees / 2
        box = {}
        box['BL_lon'] = min(longs) - margin_degrees
        box['BL_lat'] = min(lats)- margin_degrees - 0.25
        box['TR_lon'] = max(longs) + margin_degrees
        box['TR_lat'] = max(lats) + margin_degrees
        # print("bottom left:", box['BL_lon'], box['BL_lat'], file=sys.stderr)
        # print("top right", box['TR_lon'], box['TR_lat'], file=sys.stderr)

        mean_latitude = np.average(lats)

        # Lambert Conformal Conic Projection
        # map = Basemap(llcrnrlon=box['BL_lon'],llcrnrlat=box['BL_lat'],urcrnrlon=box['TR_lon'],urcrnrlat=box['TR_lat'],
        #             projection='lcc',lon_0=self.longitude_degrees, lat_0=self.latitude_degrees,
        #             resolution ='f',area_thresh=1)

        # npaeqd -	North-Polar Azimuthal Equidistant Projection - https://www.axismaps.com/guide/map-projections
        # nice for overview over the north pole
        # map = Basemap(width=200000,height=400000,
        #              projection='npaeqd',lon_0=self.longitude_degrees, boundinglat=box['BL_lat'],
        #              resolution ='f')


        # Azimuthal Equidistant Projection
        map = Basemap(#width=250000,height=350000,
                    llcrnrlon=box['BL_lon'],llcrnrlat=box['BL_lat'],urcrnrlon=box['TR_lon'],urcrnrlat=box['TR_lat'],
                    projection='aeqd',lon_0=self.longitude_degrees, lat_0=mean_latitude, # self.latitude_degrees+0.75
                    resolution ='f') 

        # Transverse Mercator Projection
        # map = Basemap(width=200000,height=400000,
        #             projection='tmerc',lon_0=self.longitude_degrees, lat_0=self.latitude_degrees,
        #             resolution ='f') 

        # draw coastlines and fill-in land 
        map.drawcoastlines()
        map.drawmapboundary(fill_color='#99ffff')
        map.fillcontinents(color='#cc9966',lake_color='#99ffff')


        # --- Plot the flight / shadow path ---
        # Loop over circuits, filtering data on 'circuit' attribute, and reading colour from an array

        num_circuits = max(df['circuit'])
        print(num_circuits,file=sys.stderr)

        for c in range(num_circuits+1):

            circuit_df = df[df['circuit'] == c] # filter dataframe on circuit number

            # create array of flight path data
            df_lat_long = circuit_df[['shading-longitude','shading-latitude']]
            flight_path = list(df_lat_long.itertuples(index=False, name=None)) # https://stackoverflow.com/a/34551914

            # Convert flight path coordinates to map projection
            x, y = map(*zip(*flight_path)) # create x and y lists in map projection from array of tuples

            # plot the flight path / shade pattern
            map.plot(x, y, linewidth = 1.5, color = line_colours[c])

            # Add velocity vectors -- see https://matplotlib.org/basemap/users/examples.html
            # uproj,vproj,xx,yy = map.transform_vector(u,v,x,y,31,31,returnxy=True,masked=True) # transform vectors to projection grid.
            # Q = map.quiver(xx,yy,uproj,vproj,scale=700) # now plot.
            # qk = plt.quiverkey(Q, 0.1, 0.1, 20, '20 m/s', labelpos='W') # make quiver key.


        # --- Finalise the map ---

        # add marker for target
        marker_x, marker_y = map([self.longitude_degrees],[self.latitude_degrees])
        map.scatter(marker_x,marker_y,30,marker='o',color=dot_colour) # https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.scatter.html

        # add scale,see https://basemaptutorial.readthedocs.io/en/latest/utilities.html
        map.drawmapscale(self.longitude_degrees,box['BL_lat'] + scale_offset, self.longitude_degrees , self.latitude_degrees,100, barstyle='fancy')
 
        return map

    def outputCSV(self,result):
        # Set up CSV output
        header = ['local-time','solar altitude', 'azimuth', 'chord-length' ,'arc-distance','shading-latitude', 'shading-longitude',]
        writer = csv.writer(sys.stdout)
        writer.writerow(header)

        #for i in range( int(24 * 60 / interval) ):
        for i in range(len(result['aboveHorizon'])):
            if (result['aboveHorizon'][i]): 
                data = [result['local-time'][i],result['solarAltitude'][i], result['azimuth'][i], result['chord'][i], result['distance'][i], result['shading-latitude'][i], result['shading-longitude'][i]]
                writer.writerow(data)


def distanceFromAltitude(altitude, solar_altitude_degrees):

    # Calculate the distance from the target to the point where the line reaches the given altitude
    chord_length = altitude / np.tan(solar_altitude_degrees * np.pi / 180)  #  tan(sun.alt) = altitude / chord_length

    # account for the curve of the earth (minimal effect except over large distances)
    earth_radius = 6371 # km
    segment_angle_radians = 2 * np.arcsin( chord_length / (2 * earth_radius) ) #  https://owlcation.com/stem/How-to-Calculate-the-Arc-Length-of-a-Circle-Segment-and-Sector-Area
    distance = segment_angle_radians * earth_radius  # length of arc, see steps at https://www.omnicalculator.com/math/arc-length

    return distance, chord_length



def locatePointB(longitude_a,latitude_a,bearing,distance):
    # Calculate the coordinates of the shading position using "GCC" https://pypi.org/project/great-circle-calculator/#point_given_start_and_bearing
    point_a = (longitude_a, latitude_a) # For "gcc" a point is a tuple in the form (lon, lat)
    return gcc.point_given_start_and_bearing(point_a, bearing, distance, unit='kilometers') # return longitude_b, latitude_b