Add trees as shading for radiation script

cea/datamanagement/trees_helper.py

@@ -0,0 +1,58 @@
+import os
+import warnings
+
+import fiona
+import geopandas as gpd
+from osgeo import gdal, osr
+
+import cea.config
+import cea.inputlocator
+from cea.schemas import schemas
+
+
+def check_terrain_bounds(tree_geometries, terrain_raster):
+    terrian_projection = terrain_raster.GetProjection()
+    proj4_str = osr.SpatialReference(wkt=terrian_projection).ExportToProj4()
+
+    # minx, miny, maxx, maxy
+    geometry_bounds = tree_geometries.to_crs(proj4_str).total_bounds
+
+    (upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = terrain_raster.GetGeoTransform()
+    minx = upper_left_x
+    maxy = upper_left_y
+    maxx = minx + x_size * terrain_raster.RasterXSize
+    miny = maxy + y_size * terrain_raster.RasterYSize
+
+    if minx > geometry_bounds[0] or miny > geometry_bounds[1] or maxx < geometry_bounds[2] or maxy < \
+            geometry_bounds[3]:
+        warnings.warn("Terrain and trees geometries do not overlap. Please ")
+
+
+def verify_tree_properties(tree_df):
+    required_columns = schemas(plugins=[])["get_tree_geometry"]["schema"]["columns"].keys()
+    diff = required_columns - tree_df.columns
+
+    if len(diff) > 0:
+        raise ValueError(f"{diff} columns are missing for tree properties.")
+
+
+def main(config):
+    locator = cea.inputlocator.InputLocator(config.scenario)
+
+    trees_df = gpd.read_file(config.trees_helper.trees)
+    terrain_raster = gdal.Open(locator.get_terrain())
+
+    # Set trees to zone crs
+    with fiona.open(locator.get_zone_geometry(), 'r') as f:
+        zone_crs = f.crs
+    trees_df.to_crs(zone_crs, inplace=True)
+
+    verify_tree_properties(trees_df)
+    check_terrain_bounds(trees_df.geometry, terrain_raster)
+
+    os.makedirs(locator.get_tree_geometry_folder(), exist_ok=True)
+    trees_df.to_file(locator.get_tree_geometry())
+
+
+if __name__ == '__main__':
+    main(cea.config.Configuration())

cea/default.config

@@ -123,6 +123,13 @@ fix-overlapping-geometries = true
 fix-overlapping-geometries.type = BooleanParameter
 fix-overlapping-geometries.help = Tries to eliminate overlapping building volumes (and subsequent double counting) by fixing overlaps of ground-surface polygons and retaining only the tallest of the overlapping buildings.
 
+[trees-helper]
+trees =
+trees.type = FileParameter
+trees.extensions = shp
+trees.nullable = true
+trees.help = Path to shapefile containing the geometries of tree canopies in polygon format and attributes of Name (unique tree ID), height_tc (tree canopy height above ground in metres), and density_la (leaf area density).
+
 [weather-helper]
 weather =
 weather.type = WeatherPathParameter

cea/inputlocator.py

@@ -536,6 +536,12 @@ def get_building_properties_folder(self):
     def get_terrain_folder(self):
         return os.path.join(self.scenario, 'inputs', 'topography')
 
+    def get_tree_geometry_folder(self):
+        return os.path.join(self.scenario, 'inputs', 'tree-geometry')
+
+    def get_tree_geometry(self):
+        return os.path.join(self.scenario, 'inputs', 'tree-geometry', 'trees.shp')
+
     def get_zone_geometry(self):
         """scenario/inputs/building-geometry/zone.shp"""
         shapefile_path = os.path.join(self.get_building_geometry_folder(), 'zone.shp')

cea/interfaces/dashboard/api/inputs.py

@@ -41,7 +41,8 @@
     ('indoor-comfort', 'get_building_comfort'),
     ('air-conditioning-systems', 'get_building_air_conditioning'),
     ('supply-systems', 'get_building_supply'),
-    ('surroundings', 'get_surroundings_geometry')
+    ('surroundings', 'get_surroundings_geometry'),
+    ('trees', "get_tree_geometry")
 ]
 
 
@@ -61,7 +62,7 @@ def get_input_database_schemas():
 
 INPUTS = get_input_database_schemas()
 INPUT_KEYS = INPUTS.keys()
-GEOJSON_KEYS = ['zone', 'surroundings', 'streets', 'dc', 'dh']
+GEOJSON_KEYS = ['zone', 'surroundings', 'trees', 'streets', 'dc', 'dh']
 NETWORK_KEYS = ['dc', 'dh']
 
 
@@ -126,6 +127,7 @@ def get(self):
         store['geojsons']['zone'], store['crs']['zone'] = df_to_json(locator.get_zone_geometry())
         store['geojsons']['surroundings'], store['crs']['surroundings'] = df_to_json(
             locator.get_surroundings_geometry())
+        store['geojsons']['trees'], store['crs']['trees'] = df_to_json(locator.get_tree_geometry())
         store['geojsons']['streets'], store['crs']['streets'] = df_to_json(locator.get_street_network())
         store['geojsons']['dc'], store['connected_buildings']['dc'], store['crs']['dc'] = get_network(config, 'dc')
         store['geojsons']['dh'], store['connected_buildings']['dh'],  store['crs']['dh'] = get_network(config, 'dh')
@@ -328,6 +330,9 @@ def df_to_json(file_location):
             lat, lon = get_lat_lon_projected_shapefile(table_df)
             crs = get_projected_coordinate_system(lat, lon)
 
+        if "Name" in table_df.columns:
+            table_df['Name'] = table_df['Name'].astype('str')
+
         # make sure that the geojson is coded in latitude / longitude
         out = table_df.to_crs(get_geographic_coordinate_system())
         out = json.loads(out.to_json())

cea/interfaces/dashboard/api/project.py

@@ -18,6 +18,7 @@
 import cea.inputlocator
 import cea.api
 import cea.config
+from cea.plots.colors import color_to_rgb
 from cea.utilities.standardize_coordinates import get_geographic_coordinate_system
 
 api = Namespace('Project', description='Current project for CEA')
@@ -296,20 +297,20 @@ def get(self, scenario):
                         zone_df = zone_df.to_crs(get_geographic_coordinate_system())
                         polygons = zone_df['geometry']
 
-                        m = StaticMap(256, 160)
+                        m = StaticMap(256, 160, url_template='http://a.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png')
                         if len(polygons) <= building_limit:
                             polygons = [list(polygons.geometry.exterior[row_id].coords) for row_id in
                                         range(polygons.shape[0])]
                             for polygon in polygons:
-                                out = Polygon(polygon, 'blue', 'black', False)
+                                out = Polygon(polygon, color_to_rgb('purple'), 'black', False)
                                 m.add_polygon(out)
                         else:
                             print(f'Number of buildings({len(polygons)}) exceed building limit({building_limit}): '
                                   f'Generating simplified image')
                             # Generate only the shape outline of the zone area
                             convex_hull = polygons.unary_union.convex_hull
                             polygon = convex_hull.exterior.coords
-                            out = Polygon(polygon, None, 'blue', False)
+                            out = Polygon(polygon, None, color_to_rgb('purple'), False)
                             m.add_polygon(out)
 
                         image = m.render()

cea/resources/radiation/geometry_generator.py

@@ -576,6 +576,28 @@ def check_terrain_bounds(zone_df, surroundings_df, terrain_raster):
         raise ValueError('Terrain provided does not cover all building geometries')
 
 
+def tree_geometry_generator(tree_df, terrain_raster):
+    terrian_projection = terrain_raster.GetProjection()
+    proj4_str = osr.SpatialReference(wkt=terrian_projection).ExportToProj4()
+    tree_df = tree_df.to_crs(proj4_str)
+
+    elevation_map = ElevationMap.read_raster(terrain_raster)
+
+    from multiprocessing.pool import Pool
+    from multiprocessing import cpu_count
+
+    with Pool(cpu_count() - 1) as pool:
+        surfaces = [
+            fetch.faces_frm_solid(result) for result in pool.starmap(
+                process_geometries, (
+                    (geom, elevation_map, (0, 1), z) for geom, z in zip(tree_df['geometry'], tree_df['height_tc'])
+                )
+            )
+        ]
+
+    return surfaces
+
+
 def geometry_main(config, zone_df, surroundings_df, terrain_raster, architecture_wwr_df, geometry_pickle_dir):
     print("Standardizing coordinate systems")
     zone_df, surroundings_df, terrain_raster = standardize_coordinate_systems(zone_df, surroundings_df, terrain_raster)

cea/resources/radiation/main.py

@@ -159,7 +159,7 @@ def main(config):
                                                                                                           architecture_wwr_df,
                                                                                                           geometry_staging_location)
 
-    daysim_staging_location = os.path.join(locator.get_temporary_folder(), 'cea_radiation')
+    daysim_staging_location = os.path.join(locator.get_solar_radiation_folder(), 'cea_radiation')
     cea_daysim = CEADaySim(daysim_staging_location, daysim_bin_path, daysim_lib_path)
 
     # create radiance input files
@@ -169,6 +169,14 @@ def main(config):
     cea_daysim.create_radiance_geometry(geometry_terrain, building_surface_properties, zone_building_names,
                                         surroundings_building_names, geometry_staging_location)
 
+    trees_path = locator.get_tree_geometry()
+    if os.path.exists(trees_path):
+        trees_df = gpd.GeoDataFrame.from_file(trees_path)
+        tree_surfaces = geometry_generator.tree_geometry_generator(trees_df, terrain_raster)
+        print("Creating radiance shading file")
+        tree_lad = trees_df["density_la"]
+        cea_daysim.create_radiance_shading(tree_surfaces, tree_lad)
+
     print("Converting files for DAYSIM")
     weather_file = locator.get_weather_file()
     print('Transforming weather files to daysim format')

cea/resources/radiation/radiance.py

@@ -50,6 +50,7 @@ def __init__(self, staging_path, daysim_dir, daysim_lib):
         self.daysim_material_path = os.path.join(self.common_inputs, 'daysim_material.rad')
         self.daysim_geometry_path = os.path.join(self.common_inputs, 'daysim_geometry.rad')
         self.wea_weather_path = os.path.join(self.common_inputs, 'weather_60min.wea')
+        self.daysim_shading_path = os.path.join(self.common_inputs, 'daysim_shading.rad')
 
         # Header Properties
         self.site_info = None
@@ -67,7 +68,7 @@ def initialize_daysim_project(self, project_name):
         """
         return DaySimProject(project_name, self.projects_dir, self.daysim_dir, self.daysim_lib,
                              self.daysim_material_path, self.daysim_geometry_path, self.wea_weather_path,
-                             self.site_info)
+                             self.site_info, self.daysim_shading_path)
 
     def create_radiance_material(self, building_surface_properties):
         add_rad_mat(self.rad_material_path, building_surface_properties)
@@ -77,6 +78,27 @@ def create_radiance_geometry(self, geometry_terrain, building_surface_properties
         create_rad_geometry(self.rad_geometry_path, geometry_terrain, building_surface_properties, zone_building_names,
                             surroundings_building_names, geometry_pickle_dir)
 
+    def create_radiance_shading(self, tree_surfaces, leaf_area_densities):
+        def tree_to_radiance(tree_id, tree_surface_list):
+            for num, occ_face in enumerate(tree_surface_list):
+                surface_name = f"tree_surface_{tree_id}_{num}"
+                yield RadSurface(surface_name, occ_face, f"tree_material_{tree_id}")
+
+        with open(self.daysim_shading_path, "w") as rad_file:
+            # # write material for trees
+            for i, tree in enumerate(tree_surfaces):
+                # light would pass through at least 2 surfaces so we divide the effect by half
+                transmissivity = math.sqrt(1 - leaf_area_densities[i])
+                string = f"void glass tree_material_{i}\n" \
+                         "0\n" \
+                         "0\n" \
+                         f"3 {transmissivity} {transmissivity} {transmissivity}"
+                rad_file.writelines(string + '\n')
+
+            for i, tree in enumerate(tree_surfaces):
+                for tree_surface_rad in tree_to_radiance(i, tree):
+                    rad_file.write(tree_surface_rad.rad())
+
     @staticmethod
     def run_cmd(cmd, daysim_dir, daysim_lib):
         print(f'Running command `{cmd}`')
@@ -164,7 +186,7 @@ def execute_radfiles2daysim(self):
 class DaySimProject(object):
     def __init__(self, project_name, project_directory, daysim_bin_directory, daysim_lib_directory,
                  daysim_material_path, daysim_geometry_path, wea_weather_path,
-                 site_info):
+                 site_info, daysim_shading_path):
 
         # Project info
         self.project_name = project_name
@@ -180,6 +202,7 @@ def __init__(self, project_name, project_directory, daysim_bin_directory, daysim
         self.daysim_geometry_path = daysim_geometry_path
         self.wea_weather_path = wea_weather_path
         self.sensor_path = os.path.join(self.project_path, "sensors.pts")
+        self.daysim_shading_path = daysim_shading_path
 
         self.hea_path = os.path.join(self.project_path, f"{project_name}.hea")
         # Header Properties
@@ -218,6 +241,10 @@ def _create_project_header_file(self):
 
             hea_file.write(header)
 
+    @property
+    def shading_exists(self):
+        return os.path.exists(self.daysim_shading_path)
+
     def cleanup_project(self):
         shutil.rmtree(self.project_path)
 
@@ -296,19 +323,87 @@ def write_radiance_parameters(self, rad_ab, rad_ad, rad_as, rad_ar, rad_aa, rad_
 
             hea_file.write(radiance_parameters)
 
-    def write_static_shading(self):
+    def generate_shading_profile(self):
+        """
+        The external shading profile should be a comma seperate file with a three line header and the format:
+        month, day, hour, shading fraction [0,1] (0=fully opened;1=fully closed) in each line.
+
+        Creates shading profile in the temp directory of the project.
+        Uses date values found in the weather file for DAYSIM
+        """
+        shading_profile = os.path.join(self.tmp_directory, "shading_profile.csv")
+        with open(self.wea_weather_path) as wea, open(shading_profile, "w") as f:
+            csv_writer = csv.writer(f, delimiter=',')
+            wea_reader = csv.reader(wea, delimiter=' ', )
+
+            # Write headers
+            csv_writer.writerows([
+                ["# Daysim annual blind schdule", "", "", ""],
+                ["# time_step 60", "comment:", "", ""],
+                ["# month", "day", "time", "shading fraction"]
+            ])
+
+            # Skip first 6 rows (headers) in weather file
+            for i in range(6):
+                next(wea_reader)
+
+            for row in wea_reader:
+                date_columns = row[:3]
+                csv_writer.writerow(date_columns + [1])
+
+        return shading_profile
+
+    def write_shading_parameters(self):
         """
-        This function writes the static shading into the header file.
+        This function writes the shading properties into the header file.
+        If no shading is found, static shading mode is used.
 
         `shading 1 <descriptive_string> <file_name.dc> <file_name.ill>`
 
         The integer 1 represents static/no shading
+
+        If shading if found, generated required files and load shading geometries
+
+        `shading -n
+        <base_file_name.dc> <base_file_name_no_blinds.ill>
+        [followed by n shading group definitions]
+        <shading_group_1_name>
+        m
+        control_keyword <shading_group_opened.rad> [followed by m lines]
+        <shading_group_1_state1.rad> <shading_group_1_state1.dc> <shading_group_1_state1.ill>`
+
+        with n = 1 or 2 = number of shading groups, m = number of states in shading group
         """
         dc_file = f"{self.project_name}.dc"
         ill_file = f"{self.project_name}.ill"
+
+        if not self.shading_exists:
+            # Use static system
+            shading_parameters = f"shading 1 static_system {dc_file} {ill_file}\n"
+        else:
+            # # Create empty shading file for base case
+            # empty_shading_file = "no_shading.rad"
+            # with open(os.path.join(self.project_path, empty_shading_file), 'w') as f:
+            #     pass
+            #
+            # # Generate shading schedule
+            # shading_profile = self.generate_shading_profile()
+            #
+            # shading_parameters = (f"shading -1\n"
+            #                       f"{dc_file} {ill_file}\n"
+            #                       f"tree_shading_group\n"
+            #                       f"1\n"
+            #                       f"AnnualShadingSchedule {shading_profile} {empty_shading_file}\n"
+            #                       f"{self.daysim_shading_path} shading_{dc_file} shading_{ill_file}")
+
+            shading_parameters = (f"shading -1\n"
+                                  f"{dc_file} {ill_file}\n"
+                                  f"tree_shading_group\n"
+                                  f"0\n"
+                                  f"ManualControl {self.daysim_shading_path}\n")
+
         with open(self.hea_path, "a") as hea_file:
-            static_shading = f"shading 1 static_system {dc_file} {ill_file}\n"
-            hea_file.write(static_shading)
+            hea_file.write(shading_parameters)
 
     def execute_gen_dc(self):
         """
@@ -321,7 +416,7 @@ def execute_gen_dc(self):
         -paste  pastes direct and diffuse daylight coefficient output files into a single complete file
         """
         # write the shading header
-        self.write_static_shading()
+        self.write_shading_parameters()
 
         command1 = f'gen_dc "{self.hea_path}" -dir'
         command2 = f'gen_dc "{self.hea_path}" -dif'
@@ -347,6 +442,8 @@ def eval_ill(self):
         """
 
         ill_path = os.path.join(self.project_path, f"{self.project_name}.ill")
+        # if self.shading_exists:
+        #     ill_path = os.path.join(self.project_path, f"shading_{self.project_name}.ill")
         with open(ill_path) as f:
             reader = csv.reader(f, delimiter=' ')
             data = np.array([np.array(row[4:], dtype=np.float32) for row in reader]).T
@@ -470,8 +567,6 @@ def add_rad_mat(daysim_mat_file, ageometry_table):
                 write_file.writelines('\n' + string + '\n')
                 written_mat_name_list.add(mat_name)
 
-        write_file.close()
-
 
 def create_rad_geometry(file_path, geometry_terrain, building_surface_properties, zone_building_names,
                         surroundings_building_names, geometry_pickle_dir):