hub/city_model_structure/building_demand/surface.py

@@ -46,6 +46,7 @@ class Surface:
     self._vegetation = None
     self._percentage_shared = None
     self._solar_collectors_area_reduction_factor = None
+    self._global_irradiance_tilted = {}
 
   @property
   def name(self):
@@ -384,3 +385,19 @@ class Surface:
     :param value: float
     """
     self._solar_collectors_area_reduction_factor = value
+
+  @property
+  def global_irradiance_tilted(self) -> dict:
+    """
+    Get global irradiance on a tilted surface in J/m2
+    :return: dict
+    """
+    return self._global_irradiance_tilted
+
+  @global_irradiance_tilted.setter
+  def global_irradiance_tilted(self, value):
+    """
+    Set global irradiance on a tilted surface in J/m2
+    :param value: dict
+    """
+    self._global_irradiance_tilted = value

hub/city_model_structure/city_object.py

@@ -41,9 +41,10 @@ class CityObject:
     self._ground_temperature = {}
     self._global_horizontal = {}
     self._diffuse = {}
-    self._beam = {}
+    self._direct_normal = {}
     self._sensors = []
     self._neighbours = None
+    self._beam = {}
 
   @property
   def level_of_detail(self) -> LevelOfDetail:
@@ -238,20 +239,20 @@ class CityObject:
     self._diffuse = value
 
   @property
-  def beam(self) -> dict:
+  def direct_normal(self) -> dict:
     """
     Get beam radiation surrounding the city object in J/m2
     :return: dict{dict{[float]}}
     """
-    return self._beam
+    return self._direct_normal
 
-  @beam.setter
-  def beam(self, value):
+  @direct_normal.setter
+  def direct_normal(self, value):
     """
     Set beam radiation surrounding the city object in J/m2
     :param value: dict{dict{[float]}}
     """
-    self._beam = value
+    self._direct_normal = value
 
   @property
   def lower_corner(self):
@@ -302,3 +303,19 @@ class CityObject:
     Set the list of neighbour_objects and their properties associated to the current city_object
     """
     self._neighbours = value
+
+  @property
+  def beam(self) -> dict:
+    """
+    Get beam radiation surrounding the city object in J/m2
+    :return: dict{dict{[float]}}
+    """
+    return self._beam
+
+  @beam.setter
+  def beam(self, value):
+    """
+    Set beam radiation surrounding the city object in J/m2
+    :param value: dict{dict{[float]}}
+    """
+    self._beam = value

hub/city_model_structure/energy_systems/pv_generation_system.py

@@ -26,6 +26,10 @@ class PvGenerationSystem(GenerationSystem):
     self._width = None
     self._height = None
     self._electricity_power = None
+    self._tilt_angle = None
+    self._surface_azimuth = None
+    self._solar_altitude_angle = None
+    self._solar_azimuth_angle = None
 
   @property
   def nominal_electricity_output(self):
@@ -202,3 +206,35 @@ class PvGenerationSystem(GenerationSystem):
     :param value: float
     """
     self._electricity_power = value
+
+  @property
+  def tilt_angle(self):
+    """
+    Get tilt angle of PV system in degrees
+    :return: float
+    """
+    return self._tilt_angle
+
+  @tilt_angle.setter
+  def tilt_angle(self, value):
+    """
+    Set PV system tilt angle in degrees
+    :param value: float
+    """
+    self._tilt_angle = value
+
+  @property
+  def surface_azimuth(self):
+    """
+    Get surface azimuth angle of PV system in degrees. 0 is North
+    :return: float
+    """
+    return self._surface_azimuth
+
+  @surface_azimuth.setter
+  def surface_azimuth(self, value):
+    """
+    Set PV system tilt angle in degrees
+    :param value: float
+    """
+    self._surface_azimuth = value

hub/exports/formats/simplified_radiosity_algorithm.py

@@ -67,7 +67,7 @@ class SimplifiedRadiosityAlgorithm:
             i = (total_days + day - 1) * 24 + hour - 1
           representative_building = self._city.buildings[0]
           _global = representative_building.diffuse[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES
-          _beam = representative_building.beam[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES
+          _beam = representative_building.direct_normal[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES
           content += f'{day}  {month}  {hour}  {_global}  {_beam}\n'
     with open(file, 'w', encoding='utf-8') as file:
       file.write(content)

hub/imports/energy_systems/montreal_future_energy_systems_parameters.py

@@ -82,8 +82,7 @@ class MontrealFutureEnergySystemParameters:
 
     return _generic_energy_systems
 
-  @staticmethod
-  def _create_generation_systems(archetype_system):
+  def _create_generation_systems(self, archetype_system):
     _generation_systems = []
     archetype_generation_systems = archetype_system.generation_systems
     if archetype_generation_systems is not None:
@@ -107,6 +106,7 @@ class MontrealFutureEnergySystemParameters:
           _generation_system.cell_temperature_coefficient = archetype_generation_system.cell_temperature_coefficient
           _generation_system.width = archetype_generation_system.width
           _generation_system.height = archetype_generation_system.height
+          _generation_system.tilt_angle = self._city.latitude
           _generic_storage_system = None
           if archetype_generation_system.energy_storage_systems is not None:
             _generic_storage_system = ElectricalStorageSystem()

hub/imports/results/simplified_radiosity_algorithm.py

@@ -34,7 +34,7 @@ class SimplifiedRadiosityAlgorithm:
     for key in self._results:
       _irradiance = {}
       header_name = key.split(':')
-      result = [x for x in self._results[key]]
+      result = [x * cte.WATTS_HOUR_TO_JULES for x in self._results[key]]
       city_object_name = header_name[1]
       building = self._city.city_object(city_object_name)
       surface_id = header_name[2]

hub/imports/weather/epw_weather_parameters.py

@@ -114,10 +114,14 @@ class EpwWeatherParameters:
                                               for x in self._weather_values['global_horizontal_radiation_wh_m2']]
       building.diffuse[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES
                                     for x in self._weather_values['diffuse_horizontal_radiation_wh_m2']]
-      building.beam[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES
-                                 for x in self._weather_values['direct_normal_radiation_wh_m2']]
+      building.direct_normal[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES
+                                          for x in self._weather_values['direct_normal_radiation_wh_m2']]
+      building.beam[cte.HOUR] = [building.global_horizontal[cte.HOUR][i] -
+                                 building.diffuse[cte.HOUR][i]
+                                 for i in range(len(building.global_horizontal[cte.HOUR]))]
       building.cold_water_temperature[cte.HOUR] = wh().cold_water_temperature(building.external_temperature[cte.HOUR])
 
+
     # create the monthly and yearly values out of the hourly
     for building in self._city.buildings:
       building.external_temperature[cte.MONTH] = \

main.py

@@ -1,7 +1,6 @@
 from scripts.geojson_creator import process_geojson
 from pathlib import Path
 import subprocess
-from scripts.ep_run_enrich import energy_plus_workflow
 from hub.imports.geometry_factory import GeometryFactory
 from hub.helpers.dictionaries import Dictionaries
 from hub.imports.construction_factory import ConstructionFactory
@@ -18,8 +17,11 @@ from scripts.costs.cost import Cost
 from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV
 import hub.helpers.constants as cte
 from hub.exports.exports_factory import ExportsFactory
+import csv
+from scripts.solar_angles import CitySolarAngles
+from scripts.radiation_tilted import RadiationTilted
 # Specify the GeoJSON file path
-# geojson_file = process_geojson(x=-73.5953602192335, y=45.492414530022515, diff=0.001)
+geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.001)
 file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson')
 # Specify the output path for the PDF file
 output_path = (Path(__file__).parent / 'out_files').resolve()
@@ -30,9 +32,8 @@ city = GeometryFactory('geojson',
                        year_of_construction_field='year_of_construction',
                        function_field='function',
                        function_to_hub=Dictionaries().montreal_function_to_hub_function).city
-# Enrich city data
+# # Enrich city data
 ConstructionFactory('nrcan', city).enrich()
-
 UsageFactory('nrcan', city).enrich()
 WeatherFactory('epw', city).enrich()
 energy_plus_workflow(city)

scripts/pv_sizing_and_simulation.py

@@ -0,0 +1,21 @@
+from scripts.radiation_tilted import RadiationTilted
+from scripts.solar_angles import CitySolarAngles
+import hub.helpers.constants as cte
+
+
+class PVSizingSimulation:
+  def __init__(self, city, tilt_angle):
+    self.city = city
+    self.tilt_angle = tilt_angle
+    self.solar_angles = CitySolarAngles(file_name=self.city.name,
+                                        location_latitude=self.city.latitude,
+                                        location_longitude=self.city.longitude,
+                                        tilt_angle=self.tilt_angle,
+                                        surface_azimuth_angle=180,
+                                        standard_meridian=-75).calculate
+    self.enrich_radiation_data()
+
+  def enrich_radiation_data(self):
+    for building in self.city.buildings:
+      roof_horizontal = [x / cte.WATTS_HOUR_TO_JULES for x in building.roofs[0].global_irradiance[cte.HOUR]]
+      RadiationTilted(building, self.solar_angles, tilt_angle=self.tilt_angle, ghi=roof_horizontal).enrich()

scripts/radiation_tilted.py

@@ -0,0 +1,109 @@
+import pandas as pd
+import math
+import hub.helpers.constants as cte
+from hub.helpers.monthly_values import MonthlyValues
+
+
+class RadiationTilted:
+  def __init__(self, building, solar_angles, tilt_angle, ghi, solar_constant=1366.1, maximum_clearness_index=1,
+               min_cos_zenith=0.065, maximum_zenith_angle=87):
+    self.building = building
+    self.ghi = ghi
+    self.tilt_angle = tilt_angle
+    self.zeniths = solar_angles['zenith'].tolist()[:-1]
+    self.incidents = solar_angles['incident angle'].tolist()[:-1]
+    self.date_time = solar_angles['DateTime'].tolist()[:-1]
+    data = {'DateTime': self.date_time, 'zenith': self.zeniths, 'incident angle': self.incidents, 'ghi': self.ghi}
+    self.df = pd.DataFrame(data)
+    self.df['DateTime'] = pd.to_datetime(self.df['DateTime'])
+    self.df.set_index('DateTime', inplace=True)
+    self.solar_constant = solar_constant
+    self.maximum_clearness_index = maximum_clearness_index
+    self.min_cos_zenith = min_cos_zenith
+    self.maximum_zenith_angle = maximum_zenith_angle
+    self.i_on = []
+    self.i_oh = []
+    self.k_t = []
+    self.fraction_diffuse = []
+    self.diffuse_horizontal = []
+    self.beam_horizontal = []
+    self.dni = []
+    self.beam_tilted = []
+    self.diffuse_tilted = []
+    self.total_radiation_tilted = []
+    self.calculate()
+
+  def dni_extra(self):
+    for i in range(len(self.df)):
+      self.i_on.append(self.solar_constant * (1 + 0.033 * math.cos(math.radians(360 * self.df.index.dayofyear[i] / 365))))
+
+    self.df['extraterrestrial normal radiation (Wh/m2)'] = self.i_on
+
+  def clearness_index(self):
+    for i in range(len(self.df)):
+      self.i_oh.append(self.i_on[i] * max(math.cos(math.radians(self.zeniths[i])), self.min_cos_zenith))
+      self.k_t.append(self.ghi[i] / self.i_oh[i])
+      self.k_t[i] = max(0, self.k_t[i])
+      self.k_t[i] = min(self.maximum_clearness_index, self.k_t[i])
+    self.df['extraterrestrial radiation on horizontal (Wh/m2)'] = self.i_oh
+    self.df['clearness index'] = self.k_t
+
+  def diffuse_fraction(self):
+    for i in range(len(self.df)):
+      if self.k_t[i] <= 0.22:
+        self.fraction_diffuse.append(1 - 0.09 * self.k_t[i])
+      elif self.k_t[i] <= 0.8:
+        self.fraction_diffuse.append(0.9511 - 0.1604 * self.k_t[i] + 4.388 * self.k_t[i] ** 2 -
+                                     16.638 * self.k_t[i] ** 3 + 12.336 * self.k_t[i] ** 4)
+      else:
+        self.fraction_diffuse.append(0.165)
+      if self.zeniths[i] > self.maximum_zenith_angle:
+        self.fraction_diffuse[i] = 1
+
+    self.df['diffuse fraction'] = self.fraction_diffuse
+
+  def radiation_components_horizontal(self):
+    for i in range(len(self.df)):
+      self.diffuse_horizontal.append(self.ghi[i] * self.fraction_diffuse[i])
+      self.beam_horizontal.append(self.ghi[i] - self.diffuse_horizontal[i])
+      self.dni.append((self.ghi[i] - self.diffuse_horizontal[i]) / math.cos(math.radians(self.zeniths[i])))
+      if self.zeniths[i] > self.maximum_zenith_angle or self.dni[i] < 0:
+        self.dni[i] = 0
+
+    self.df['diffuse horizontal (Wh/m2)'] = self.diffuse_horizontal
+    self.df['dni (Wh/m2)'] = self.dni
+    self.df['beam horizontal (Wh/m2)'] = self.beam_horizontal
+
+  def radiation_components_tilted(self):
+    for i in range(len(self.df)):
+      self.beam_tilted.append(self.dni[i] * math.cos(math.radians(self.incidents[i])))
+      self.beam_tilted[i] = max(self.beam_tilted[i], 0)
+      self.diffuse_tilted.append(self.diffuse_horizontal[i] * ((1 + math.cos(math.radians(self.tilt_angle))) / 2))
+      self.total_radiation_tilted.append(self.beam_tilted[i] + self.diffuse_tilted[i])
+
+    self.df['beam tilted (Wh/m2)'] = self.beam_tilted
+    self.df['diffuse tilted (Wh/m2)'] = self.diffuse_tilted
+    self.df['total radiation tilted (Wh/m2)'] = self.total_radiation_tilted
+
+  def calculate(self) -> pd.DataFrame:
+    self.dni_extra()
+    self.clearness_index()
+    self.diffuse_fraction()
+    self.radiation_components_horizontal()
+    self.radiation_components_tilted()
+    return self.df
+
+  def enrich(self):
+    tilted_radiation = self.total_radiation_tilted
+    print(len(tilted_radiation))
+    self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in
+                                                                 tilted_radiation]
+    self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in
+                                                                 tilted_radiation]
+    self.building.roofs[0].global_irradiance_tilted[cte.MONTH] = (
+      MonthlyValues.get_total_month(self.building.roofs[0].global_irradiance_tilted[cte.HOUR]))
+    self.building.roofs[0].global_irradiance_tilted[cte.YEAR] = \
+      [sum(self.building.roofs[0].global_irradiance_tilted[cte.MONTH])]
+
+
+

scripts/solar_angles.py

@@ -0,0 +1,146 @@
+import math
+import pandas as pd
+from datetime import datetime
+from pathlib import Path
+
+
+class CitySolarAngles:
+  def __init__(self, file_name, location_latitude, location_longitude, tilt_angle, surface_azimuth_angle,
+               standard_meridian=-75):
+    self.file_name = file_name
+    self.location_latitude = location_latitude
+    self.location_longitude = location_longitude
+    self.location_latitude_rad = math.radians(location_latitude)
+    self.surface_azimuth_angle = surface_azimuth_angle
+    self.surface_azimuth_rad = math.radians(surface_azimuth_angle)
+    self.tilt_angle = tilt_angle
+    self.tilt_angle_rad = math.radians(tilt_angle)
+    self.standard_meridian = standard_meridian
+    self.longitude_correction = (location_longitude - standard_meridian) * 4
+    self.timezone = 'Etc/GMT+5'
+
+    self.eot = []
+    self.ast = []
+    self.hour_angles = []
+    self.declinations = []
+    self.solar_altitudes = []
+    self.solar_azimuths = []
+    self.zeniths = []
+    self.incidents = []
+    self.beam_tilted = []
+    self.factor = []
+    self.times = pd.date_range(start='2023-01-01', end='2024-01-01', freq='H', tz=self.timezone)
+    self.df = pd.DataFrame(index=self.times)
+    self.day_of_year = self.df.index.dayofyear
+
+  def solar_time(self, datetime_val, day_of_year):
+    b = (day_of_year - 81) * 2 * math.pi / 364
+    eot = 9.87 * math.sin(2 * b) - 7.53 * math.cos(b) - 1.5 * math.sin(b)
+    self.eot.append(eot)
+
+    # Calculate Local Solar Time (LST)
+    lst_hour = datetime_val.hour
+    lst_minute = datetime_val.minute
+    lst_second = datetime_val.second
+    lst = lst_hour + lst_minute / 60 + lst_second / 3600
+
+    # Calculate Apparent Solar Time (AST) in decimal hours
+    ast_decimal = lst + eot / 60 + self.longitude_correction / 60
+    ast_hours = int(ast_decimal)
+    ast_minutes = round((ast_decimal - ast_hours) * 60)
+
+    # Ensure ast_minutes is within valid range
+    if ast_minutes == 60:
+      ast_hours += 1
+      ast_minutes = 0
+    elif ast_minutes < 0:
+      ast_minutes = 0
+    ast_time = datetime(year=datetime_val.year, month=datetime_val.month, day=datetime_val.day,
+                        hour=ast_hours, minute=ast_minutes)
+    self.ast.append(ast_time)
+    return ast_time
+
+  def declination_angle(self, day_of_year):
+    declination = 23.45 * math.sin(math.radians(360 / 365 * (284 + day_of_year)))
+    declination_radian = math.radians(declination)
+    self.declinations.append(declination)
+    return declination_radian
+
+  def hour_angle(self, ast_time):
+    hour_angle = ((ast_time.hour * 60 + ast_time.minute) - 720) / 4
+    hour_angle_radian = math.radians(hour_angle)
+    self.hour_angles.append(hour_angle)
+    return hour_angle_radian
+
+  def solar_altitude(self, declination_radian, hour_angle_radian):
+    solar_altitude_radians = math.asin(math.cos(self.location_latitude_rad) * math.cos(declination_radian) *
+                                       math.cos(hour_angle_radian) + math.sin(self.location_latitude_rad) *
+                                       math.sin(declination_radian))
+    solar_altitude = math.degrees(solar_altitude_radians)
+    self.solar_altitudes.append(solar_altitude)
+    return solar_altitude_radians
+
+  def zenith(self, solar_altitude_radians):
+    solar_altitude = math.degrees(solar_altitude_radians)
+    zenith_degree = 90 - solar_altitude
+    zenith_radian = math.radians(zenith_degree)
+    self.zeniths.append(zenith_degree)
+    return zenith_radian
+
+  def solar_azimuth_analytical(self, hourangle, declination, zenith):
+    numer = (math.cos(zenith) * math.sin(self.location_latitude_rad) - math.sin(declination))
+    denom = (math.sin(zenith) * math.cos(self.location_latitude_rad))
+    if math.isclose(denom, 0.0, abs_tol=1e-8):
+      cos_azi = 1.0
+    else:
+      cos_azi = numer / denom
+
+    cos_azi = max(-1.0, min(1.0, cos_azi))
+
+    sign_ha = math.copysign(1, hourangle)
+    solar_azimuth_radians = sign_ha * math.acos(cos_azi) + math.pi
+    solar_azimuth_degrees = math.degrees(solar_azimuth_radians)
+    self.solar_azimuths.append(solar_azimuth_degrees)
+    return solar_azimuth_radians
+
+  def incident_angle(self, solar_altitude_radians, solar_azimuth_radians):
+    incident_radian = math.acos(math.cos(solar_altitude_radians) *
+                                math.cos(abs(solar_azimuth_radians - self.surface_azimuth_rad)) *
+                                math.sin(self.tilt_angle_rad) + math.sin(solar_altitude_radians) *
+                                math.cos(self.tilt_angle_rad))
+    incident_angle_degrees = math.degrees(incident_radian)
+    self.incidents.append(incident_angle_degrees)
+    return incident_radian
+
+  @property
+  def calculate(self) -> pd.DataFrame:
+    for i in range(len(self.times)):
+      datetime_val = self.times[i]
+      day_of_year = self.day_of_year[i]
+      declination_radians = self.declination_angle(day_of_year)
+      ast_time = self.solar_time(datetime_val, day_of_year)
+      hour_angle_radians = self.hour_angle(ast_time)
+      solar_altitude_radians = self.solar_altitude(declination_radians, hour_angle_radians)
+      zenith_radians = self.zenith(solar_altitude_radians)
+      solar_azimuth_radians = self.solar_azimuth_analytical(hour_angle_radians, declination_radians, zenith_radians)
+      incident_angle_radian = self.incident_angle(solar_altitude_radians, solar_azimuth_radians)
+
+    self.df['DateTime'] = self.times
+    self.df['AST'] = self.ast
+    self.df['hour angle'] = self.hour_angles
+    self.df['eot'] = self.eot
+    self.df['declination angle'] = self.declinations
+    self.df['solar altitude'] = self.solar_altitudes
+    self.df['zenith'] = self.zeniths
+    self.df['solar azimuth'] = self.solar_azimuths
+    self.df['incident angle'] = self.incidents
+
+    return self.df
+
+
+
+
+
+
+
+

tests/test_construction_factory.py

@@ -81,7 +81,7 @@ class TestConstructionFactory(TestCase):
       self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated')
       self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated')
       self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated')
-      self.assertEqual(len(building.beam), 0, 'building beam is calculated')
+      self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated')
       self.assertIsNotNone(building.lower_corner, 'building lower corner is none')
       self.assertEqual(len(building.sensors), 0, 'building sensors are assigned')
       self.assertIsNotNone(building.internal_zones, 'no internal zones created')

tests/test_geometry_factory.py

@@ -52,7 +52,7 @@ class TestGeometryFactory(TestCase):
       self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated')
       self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated')
       self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated')
-      self.assertEqual(len(building.beam), 0, 'building beam is calculated')
+      self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated')
       self.assertIsNotNone(building.lower_corner, 'building lower corner is none')
       self.assertEqual(len(building.sensors), 0, 'building sensors are assigned')
       self.assertIsNotNone(building.internal_zones, 'no internal zones created')

tests/test_usage_factory.py

@@ -44,7 +44,7 @@ class TestUsageFactory(TestCase):
       self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated')
       self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated')
       self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated')
-      self.assertEqual(len(building.beam), 0, 'building beam is calculated')
+      self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated')
       self.assertIsNotNone(building.lower_corner, 'building lower corner is none')
       self.assertEqual(len(building.sensors), 0, 'building sensors are assigned')
       self.assertIsNotNone(building.internal_zones, 'no internal zones created')