diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 77d30c59da..1f28cfff7c 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -28,4 +28,6 @@ jobs: - name: Run tests shell: bash -l {0} - run: cea test --workflow slow + run: | + cea test --type unittest + cea test --type integration diff --git a/.github/workflows/pull.request.yml b/.github/workflows/pull.request.yml index e715fe394b..172cb8c669 100644 --- a/.github/workflows/pull.request.yml +++ b/.github/workflows/pull.request.yml @@ -66,4 +66,6 @@ jobs: - name: Run tests shell: bash -l {0} - run: cea test --workflow quick + run: | + cea test --type unittest + cea test --type integration diff --git a/.github/workflows/setup_build.yml b/.github/workflows/setup_build.yml index 49e0dd374a..a7bb49674b 100644 --- a/.github/workflows/setup_build.yml +++ b/.github/workflows/setup_build.yml @@ -114,9 +114,13 @@ jobs: path: gui repository: architecture-building-systems/CityEnergyAnalyst-GUI - - uses: mamba-org/setup-micromamba@v1 - with: - micromamba-binary-path: gui/dependencies/micromamba + - name: Fetch micromamba + run: | + cd gui + mkdir -p ./dependencies/arm64 + curl -Ls https://micro.mamba.pm/api/micromamba/osx-arm64/latest | tar -xvj -C ./dependencies/arm64 --strip-components=1 bin/micromamba + mkdir -p ./dependencies/x64 + curl -Ls https://micro.mamba.pm/api/micromamba/osx-64/latest | tar -xvj -C ./dependencies/x64 --strip-components=1 bin/micromamba - name: Save Apple API Key secret to file env: diff --git a/CHANGELOG.md b/CHANGELOG.md index f15a2945d4..4da3f66a7b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,5 @@ +- 2024-05-27 - 3.36.0 - #3586 Add nullable info for column schema +- 2024-05-21 - 3.36.0 - #3584 Release 3.36.0 - 2024-05-16 - 3.35.6 - #3583 Use MKL for Windows to solve memory errors - 2024-05-15 - 3.35.6 - #3574 Add trees as shading for radiation script - 2024-05-14 - 3.35.6 - #3576 Consider adjacent walls by default diff --git a/CREDITS.md b/CREDITS.md index c081abf8b3..6d1b95e70e 100644 --- a/CREDITS.md +++ b/CREDITS.md @@ -1,7 +1,7 @@ How to Cite =========== -The CEA Team. (2024). City Energy Analyst (v3.36.0). Zenodo. https://doi.org/10.5281/zenodo.10950383 +The CEA Team. (2024). City Energy Analyst (v3.36.1). Zenodo. https://doi.org/10.5281/zenodo.11382063 The CEA team ============ @@ -11,6 +11,25 @@ and split into the categories of developers, product owner, project sponsor and +- Version 3.36.1 - May 2024 + + Developers: + * Reynold Mok + * Martín Mosteiro Romero + * Mathias Niffeler + * Zhongming Shi + + Product lead: + * Zhongming Shi + + Product sponsor: + * [Arno Schlueter](https://systems.arch.ethz.ch/arno-schlueter) + * [Toni Piëch Foundation](https://www.tonipiechfoundation.org/) + + Others: + * https://cityenergyanalyst.com/people + + - Version 3.36.0 - May 2024 Developers: diff --git a/README.rst b/README.rst index eab8021677..061fe4b30d 100644 --- a/README.rst +++ b/README.rst @@ -4,10 +4,8 @@ :alt: GitHub license .. |repo_size| image:: https://img.shields.io/github/repo-size/architecture-building-systems/CityEnergyAnalyst :alt: Repo Size -.. |lines_of_code| image:: https://img.shields.io/tokei/lines/github/architecture-building-systems/CityEnergyAnalyst - :alt: Lines of code -.. |zenodo| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.10950383.svg - :target: https://doi.org/10.5281/zenodo.10950383 +.. |zenodo| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.11382063.svg + :target: https://doi.org/10.5281/zenodo.11382063 .. image:: cea_logo.png :scale: 25 % @@ -33,4 +31,4 @@ We invite all CEA users to get acquainted with the CEA Dashboard and CEA Console Cite us: -------- -The CEA Team. (2024). City Energy Analyst (v3.36.0). Zenodo. https://doi.org/10.5281/zenodo.10950383 +The CEA Team. (2024). City Energy Analyst (v3.36.1). Zenodo. https://doi.org/10.5281/zenodo.11382063 diff --git a/bin/create_trace_graphviz.py b/bin/create_trace_graphviz.py index 8f7aa13092..8b9fabd1d5 100644 --- a/bin/create_trace_graphviz.py +++ b/bin/create_trace_graphviz.py @@ -12,7 +12,7 @@ import yaml import cea.config -from cea.tests.trace_inputlocator import create_graphviz_output +from cea.utilities.trace_inputlocator.trace_inputlocator import create_graphviz_output def main(config): with open(config.trace_inputlocator.yaml_output_file, 'r') as f: diff --git a/cea/__init__.py b/cea/__init__.py index 04c3d20515..146852f41a 100644 --- a/cea/__init__.py +++ b/cea/__init__.py @@ -1,4 +1,4 @@ -__version__ = "3.36.0" +__version__ = "3.36.1" class ConfigError(Exception): diff --git a/cea/datamanagement/surroundings_helper.py b/cea/datamanagement/surroundings_helper.py index 05100515ce..0d2ee2f18e 100644 --- a/cea/datamanagement/surroundings_helper.py +++ b/cea/datamanagement/surroundings_helper.py @@ -10,6 +10,7 @@ import pandas as pd from geopandas import GeoDataFrame as gdf from geopandas.tools import sjoin as spatial_join +from shapely import MultiPolygon import cea.config import cea.inputlocator @@ -35,7 +36,11 @@ def calc_surrounding_area(zone_gdf, buffer_m): :param float buffer_m: Buffer to add to zone building geometries :return: Surrounding area GeoDataFrame """ - surrounding_area = gdf(geometry=[zone_gdf.geometry.buffer(buffer_m).unary_union], crs=zone_gdf.crs) + merged_zone = zone_gdf.geometry.unary_union + if isinstance(merged_zone, MultiPolygon): + merged_zone = merged_zone.convex_hull + + surrounding_area = gdf(geometry=[merged_zone.buffer(buffer_m)], crs=zone_gdf.crs) return surrounding_area @@ -135,7 +140,6 @@ def erase_no_surrounding_areas(all_surroundings, zone, area_with_buffer): :return: GeoDataFrame with surrounding buildings """ buffer_polygon = area_with_buffer.to_crs(zone.crs).geometry.values[0] - zone_area = gdf(geometry=[zone.geometry.unary_union], crs=zone.crs) within_buffer = all_surroundings.geometry.intersects(buffer_polygon) surroundings = all_surroundings[within_buffer] @@ -144,7 +148,7 @@ def erase_no_surrounding_areas(all_surroundings, zone, area_with_buffer): if not any([s_building_footprint.intersects(z_building_footprint) for z_building_footprint in zone.geometry])] footprints_gdf = gdf(geometry=footprints_without_overlaps, crs=surroundings.crs) - relevant_surroundings = spatial_join(surroundings, footprints_gdf, op='within') + relevant_surroundings = spatial_join(surroundings, footprints_gdf, predicate='within') return relevant_surroundings.copy() @@ -175,8 +179,7 @@ def geometry_extractor_osm(locator, config): # get footprints of all the surroundings print("Getting building footprints") area_with_buffer_polygon = area_with_buffer.to_crs(get_geographic_coordinate_system()).geometry.values[0] - all_surroundings = osmnx.geometries.geometries_from_polygon(polygon=area_with_buffer_polygon, - tags={"building": True}) + all_surroundings = osmnx.features_from_polygon(polygon=area_with_buffer_polygon, tags={"building": True}) all_surroundings = all_surroundings.to_crs(get_projected_coordinate_system(float(lat), float(lon))) # erase overlapping area diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index 9dfc93c197..ee72855162 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -1,145 +1,265 @@ -""" -This script extracts terrain elevation from NASA - SRTM -https://www2.jpl.nasa.gov/srtm/ -""" - - - - +import datetime +import io +import math import os +from typing import Iterable, Tuple, Dict, Union, List +import geopandas as gpd import numpy as np import pandas as pd +import rasterio import requests -from geopandas import GeoDataFrame as Gdf -from osgeo import gdal, ogr, osr -from shapely.geometry import Polygon +from pyproj import CRS +from rasterio import MemoryFile +from rasterio.mask import mask +from rasterio.merge import merge +from rasterio.warp import calculate_default_transform, Resampling, reproject +from shapely import box import cea.config import cea.inputlocator -from cea.datamanagement.surroundings_helper import get_zone_and_surr_in_projected_crs -from cea.utilities.standardize_coordinates import get_projected_coordinate_system, get_geographic_coordinate_system - -__author__ = "Jimeno Fonseca" -__copyright__ = "Copyright 2018, Architecture and Building Systems - ETH Zurich" -__credits__ = ["Jimeno Fonseca"] -__license__ = "MIT" -__version__ = "0.1" -__maintainer__ = "Daren Thomas" -__email__ = "cea@arch.ethz.ch" -__status__ = "Production" - - -def request_elevation(lon, lat): - # script for returning elevation from lat, long, based on open elevation data - # which in turn is based on SRTM - query = ('https://api.open-elevation.com/api/v1/lookup?locations=' + str(lat) + ',' + str(lon)) - r = requests.get(query).json() # json object, various ways you can extract value - # one approach is to use pandas json functionality: - elevation = pd.io.json.json_normalize(r, 'results')['elevation'].values[0] - return elevation - - -def calc_bounding_box_projected_coordinates(locator): - # connect both files and avoid repetition - data_zone, data_dis = get_zone_and_surr_in_projected_crs(locator) - data_dis = data_dis.loc[~data_dis["Name"].isin(data_zone["Name"])] - data = pd.concat([ - data_zone.to_crs(get_geographic_coordinate_system()), - data_dis.to_crs(get_geographic_coordinate_system()) - ], ignore_index=True, sort=True) - lon = data.geometry[0].centroid.coords.xy[0][0] - lat = data.geometry[0].centroid.coords.xy[1][0] - crs = get_projected_coordinate_system(float(lat), float(lon)) - data = data.to_crs(get_projected_coordinate_system(float(lat), float(lon))) - result = data.total_bounds - result = [np.float32(x) for x in result] # in float32 so the raster works - return result, crs, lon, lat - - -def terrain_elevation_extractor(locator, config): - """this is where the action happens if it is more than a few lines in ``main``. - NOTE: ADD YOUR SCRIPT'S DOCUMENTATION HERE (how) - NOTE: RENAME THIS FUNCTION (SHOULD PROBABLY BE THE SAME NAME AS THE MODULE) + +URL_FORMAT = "https://s3.amazonaws.com/elevation-tiles-prod/geotiff/{zoom}/{x}/{y}.tif" +TILE_CRS = CRS.from_epsg(3857) +DEFAULT_CRS = CRS.from_epsg(4326) + +DATA_SOURCE_URL = "https://github.com/tilezen/joerd/blob/master/docs/data-sources.md" +ATTRIBUTION_URL = "https://github.com/tilezen/joerd/blob/master/docs/attribution.md" + +ATTRIBUTION = """ +* ArcticDEM terrain data DEM(s) were created from DigitalGlobe, Inc., imagery and + funded under National Science Foundation awards 1043681, 1559691, and 1542736; +* Australia terrain data © Commonwealth of Australia (Geoscience Australia) 2017; +* Austria terrain data © offene Daten Österreichs – Digitales Geländemodell (DGM) + Österreich; +* Canada terrain data contains information licensed under the Open Government + Licence – Canada; +* Europe terrain data produced using Copernicus data and information funded by the + European Union - EU-DEM layers; +* Global ETOPO1 terrain data U.S. National Oceanic and Atmospheric Administration +* Mexico terrain data source: INEGI, Continental relief, 2016; +* New Zealand terrain data Copyright 2011 Crown copyright (c) Land Information New + Zealand and the New Zealand Government (All rights reserved); +* Norway terrain data © Kartverket; +* United Kingdom terrain data © Environment Agency copyright and/or database right + 2015. All rights reserved; +* United States 3DEP (formerly NED) and global GMTED2010 and SRTM terrain data + courtesy of the U.S. Geological Survey. +""" + + +def get_tile_number(lat: float, lon: float, zoom: int) -> Tuple[int, int]: + """ + Converts latitude and longitude to slippy map tile coordinates. + Reference: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames """ + lat_rad = math.radians(lat) - # local variables: - elevation = config.terrain_helper.elevation - grid_size = config.terrain_helper.grid_size - extra_border = np.float32(30) # adding extra 30 m to avoid errors of no data - raster_path = locator.get_terrain() - locator.ensure_parent_folder_exists(raster_path) - - # get the bounding box coordinates - assert os.path.exists( - locator.get_surroundings_geometry()), 'Get surroundings geometry file first or the coordinates of the area where' \ - ' to extract the terrain from in the next format: lon_min, lat_min, lon_max, lat_max' - print("generating terrain from Surroundings area") - bounding_box_surroundings_file, crs, lon, lat = calc_bounding_box_projected_coordinates(locator) - x_min = bounding_box_surroundings_file[0] - extra_border - y_min = bounding_box_surroundings_file[1] - extra_border - x_max = bounding_box_surroundings_file[2] + extra_border - y_max = bounding_box_surroundings_file[3] + extra_border - - # make sure output is a whole number when min-max is divided by grid size - x_extra = grid_size - ((x_max - x_min) % grid_size)/2 - y_extra = grid_size - ((y_max - y_min) % grid_size)/2 - x_min -= x_extra - y_min -= y_extra - x_max += x_extra - y_max += y_extra - - ##TODO: get the elevation from satellite data. Open-elevation was working, but the project is dying. - # if elevation is None: - # print('extracting elevation from satellite data, this needs connection to the internet') - # elevation = request_elevation(lon, lat) - # print("Proceeding to calculate terrain file with fixed elevation in m of ", elevation) - # else: - # print("Proceeding to calculate terrain file with fixed elevation in m of ",elevation) - - print("Proceeding to calculate terrain file with fixed elevation in m of ", elevation) - - # now calculate the raster with the fixed elevation - calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, - x_max, x_min, y_max, y_min) - - -def calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max, - y_min): - # local variables: - temp_shapefile = locator.get_temporary_file("terrain.shp") - cols = int((x_max - x_min) / grid_size) - rows = int((y_max - y_min) / grid_size) - shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]]) - geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes]) - geodataframe.to_file(temp_shapefile) - # 1) opening the shapefile - source_ds = ogr.Open(temp_shapefile) - source_layer = source_ds.GetLayer() - target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32) ##COMMENT 2 - target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size)) ##COMMENT 3 - # 5) Adding a spatial reference ##COMMENT 4 - target_dsSRS = osr.SpatialReference() - target_dsSRS.ImportFromProj4(crs) - target_ds.SetProjection(target_dsSRS.ExportToWkt()) - band = target_ds.GetRasterBand(1) - band.SetNoDataValue(-9999) ##COMMENT 5 - gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation]) ##COMMENT 6 - target_ds = None # closing the file + n = 2.0 ** zoom + x_tile = int((lon + 180.0) / 360.0 * n) + y_tile = int((1.0 - math.log(math.tan(lat_rad) + 1.0 / math.cos(lat_rad)) / math.pi) / 2.0 * n) + return x_tile, y_tile -def main(config): +def get_all_tile_numbers(min_x: float, min_y: float, max_x: float, max_y: float, + zoom: int) -> Iterable[Tuple[int, int]]: + """ + Gets tile numbers based on given bounds. + """ + # Get tile numbers of bounds + tile_numbers = {get_tile_number(min_y, min_x, zoom), + get_tile_number(min_y, max_x, zoom), + get_tile_number(max_y, max_x, zoom), + get_tile_number(max_y, min_x, zoom)} + + # Get max and min of tile numbers + min_x_tile, min_y_tile, max_x_tile, max_y_tile = (None, None, None, None) + for i, tile_number in enumerate(tile_numbers): + if i == 0: + min_x_tile = tile_number[0] + max_x_tile = tile_number[0] + + min_y_tile = tile_number[1] + max_y_tile = tile_number[1] + + min_x_tile = min(min_x_tile, tile_number[0]) + max_x_tile = max(max_x_tile, tile_number[0]) + + min_y_tile = min(min_y_tile, tile_number[1]) + max_y_tile = max(max_y_tile, tile_number[1]) + + # Generate all required tile numbers + all_tile_numbers = [] + for x in range(min_x_tile, max_x_tile + 1): + for y in range(min_y_tile, max_y_tile + 1): + all_tile_numbers.append((x, y)) + + return all_tile_numbers + + +def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.Affine, Dict, List[Dict]]: + """ + Merge raster files from tile urls into a single raster. + """ + rasters = [] + tile_info = [] + try: + for tile_url in tile_urls: + with requests.get(tile_url, stream=True) as r: + r.raise_for_status() + data = io.BytesIO(r.content) + + # Load raster data using rasterio + raster = rasterio.open(data) + rasters.append(raster) + + # Store tile info + tile_info.append({ + "url": tile_url, + "headers": r.headers + }) + + # Merge rasters + dest, transform = merge(rasters) + meta = rasters[0].meta.copy() + meta.update({ + "driver": "GTiff", + "height": dest.shape[1], + "width": dest.shape[2], + "transform": transform, + }) + finally: + for raster in rasters: + raster.close() + + return dest, transform, meta, tile_info + + +def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, + dst_crs: Union[CRS, dict], grid_size: int) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + """ + Reproject raster array to specified CRS. """ - Create the terrain.tif file + # Get bounds from transform + minx, miny = src_transform * (0, src_array.shape[2]) + maxx, maxy = src_transform * (src_array.shape[1], 0) + + transform, width, height = calculate_default_transform( + meta["crs"], dst_crs, src_array.shape[2], src_array.shape[1], + minx, miny, maxx, maxy, + resolution=(grid_size, grid_size) + ) - :param config: - :type config: cea.config.Configuration - :return: + new_meta = meta.copy() + new_meta.update({ + 'crs': dst_crs, + 'transform': transform, + 'width': width, + 'height': height + }) + + dst_array = np.empty((src_array.shape[0], height, width), dtype=src_array.dtype) + for i in range(src_array.shape[0]): + reproject( + source=src_array[i], + destination=dst_array[i], + src_transform=src_transform, + src_crs=meta["crs"], + dst_transform=transform, + dst_crs=dst_crs, + resampling=Resampling.nearest + ) + + return dst_array, transform, new_meta + + +def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12, grid_size: int = 30, + src_crs: Union[CRS, dict] = DEFAULT_CRS) -> Tuple[np.ndarray, rasterio.Affine, Dict, List[Dict]]: """ - assert os.path.exists(config.scenario), 'Scenario not found: %s' % config.scenario + Fetch raster data array based on bounds in the given source CRS (Default: lat, lon). + The resulting raster data would also be in the given source CRS. + Also returns the info of the tile that were fetched. + """ + + bounding_box = gpd.GeoDataFrame(geometry=[box(min_x, min_y, max_x, max_y)], crs=src_crs) + print(f"Generating raster based on bounds: ({min_x}, {min_y}, {max_x}, {max_y}), {src_crs}") + + # Fetch tile numbers based on lat lon bounds + reprojected_bounds = bounding_box.to_crs(DEFAULT_CRS).total_bounds + tile_numbers = get_all_tile_numbers(*reprojected_bounds, zoom) + + # Get merged raster array + tile_urls = [URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers] + dest, transform, meta, tile_info = merge_raster_tiles(tile_urls) + + # Reproject raster array to bounds crs + dest, transform, meta = reproject_raster_array(dest, transform, meta, src_crs, grid_size) + + # Crop raster based on bounds + with MemoryFile() as memfile: + with memfile.open(**meta) as dataset: + dataset.write(dest) + + out_dest, out_transform = mask(dataset, bounding_box.geometry, crop=True) + out_meta = dataset.meta.copy() + out_meta.update({ + "height": out_dest.shape[1], + "width": out_dest.shape[2], + "transform": out_transform + }) + + return out_dest, out_transform, out_meta, tile_info + + +def main(config): + grid_size = config.terrain_helper.grid_size locator = cea.inputlocator.InputLocator(config.scenario) - terrain_elevation_extractor(locator, config) + # Get total bounds + zone_df = gpd.read_file(locator.get_zone_geometry()) + surroundings_df = gpd.read_file(locator.get_surroundings_geometry()).to_crs(zone_df.crs) + total_df = gpd.GeoDataFrame(pd.concat([zone_df.geometry, surroundings_df.geometry])) + total_bounds = total_df.total_bounds + + # Add buffer to bounds in meters (using projected crs), to ensure overlap + buffer = 30 + projected_crs = total_df.estimate_utm_crs() + reprojected_df = gpd.GeoDataFrame(geometry=[box(*total_bounds)], crs=zone_df.crs).to_crs(projected_crs) + buffer_df = reprojected_df.buffer(buffer) + min_x, min_y, max_x, max_y = buffer_df.total_bounds + + # Fetch tiff data + dest, transform, meta, tile_info = fetch_tiff(min_x, min_y, max_x, max_y, + grid_size=grid_size, src_crs=buffer_df.crs) + + os.makedirs(os.path.dirname(locator.get_terrain()), exist_ok=True) + # Write reference file + reference_file = os.path.join(os.path.dirname(locator.get_terrain()), "reference.txt") + tile_info_string = '\n'.join([ + f"url: {info['url']}\n" + f"data sources: {info['headers']['x-amz-meta-x-imagery-sources']}\n" + f"last modified: {info['headers']['Last-Modified']}\n" for info in tile_info + ]) + content = (f"Citation:\n" + f"Terrain Tiles was accessed on {datetime.datetime.now().date()} " + f"from https://registry.opendata.aws/terrain-tiles.\n" + f"\n" + f"Information of tiles used:\n" + f"{tile_info_string}\n" + f"Acknowledgement of Data Sources:" + f"{ATTRIBUTION}\n" + f"For more information about the data:\n" + f"Data sources: {DATA_SOURCE_URL}\n" + f"Attribution: {ATTRIBUTION_URL}\n") + + with open(reference_file, "w") as f: + f.write(content) + print(content) + print(f"Reference file is written to: {reference_file}") + + # Write to disk + with rasterio.open(locator.get_terrain(), "w", **meta) as f: + f.write(dest) if __name__ == '__main__': diff --git a/cea/datamanagement/zone_helper.py b/cea/datamanagement/zone_helper.py index 7a0fb96331..78e6707e83 100644 --- a/cea/datamanagement/zone_helper.py +++ b/cea/datamanagement/zone_helper.py @@ -225,7 +225,6 @@ def assign_attributes_additional(shapefile): return shapefile - def fix_overlapping_geoms(buildings, zone): """ This function eliminates overlapping geometries. To decide which portions of two overlapping geometries to @@ -274,7 +273,7 @@ def fix_overlapping_geoms(buildings, zone): # FIX OVERLAYS IN THE BUILDING GEOMETRIES # overlay the grid with the zone polygon, retaining the overlapping grid cells, and ... - grid = zone.overlay(grid, how="intersection") + grid = zone.to_crs(grid.crs).overlay(grid, how="intersection") # iterate through the grid cells, overlaying them with the buildings, retaining the buildings for cell_index in range(grid.geometry.size): @@ -500,12 +499,12 @@ def polygon_to_zone(buildings_floors, buildings_floors_below_ground, buildings_h lon = poly.geometry[0].centroid.coords.xy[0][0] lat = poly.geometry[0].centroid.coords.xy[1][0] # get all footprints in the district tagged as 'building' or 'building:part' in OSM - shapefile = osmnx.features.features_from_polygon(polygon=poly['geometry'].values[0], tags={"building": True}) + shapefile = osmnx.features_from_polygon(polygon=poly['geometry'].values[0], tags={"building": True}) if include_building_parts: try: # get all footprints in the district tagged as 'building' or 'building:part' in OSM - building_parts = osmnx.features.features_from_polygon(polygon=poly['geometry'].values[0], - tags={"building": ["part"]}) + building_parts = osmnx.features_from_polygon(polygon=poly['geometry'].values[0], + tags={"building": ["part"]}) shapefile = pd.concat([shapefile, building_parts], ignore_index=True) # using building:part tags requires fixing overlapping polygons if not fix_overlapping: @@ -572,7 +571,7 @@ def flatten_geometries(gdf): DISCARDED_GEOMETRY_TYPES = ['Point', 'LineString'] # Explode MultiPolygons and GeometryCollections - gdf = gdf.explode() + gdf = gdf.explode(index_parts=True) # Drop geometry types that cannot be processed by CEA gdf = gdf.loc[~ gdf.geometry.geom_type.isin(DISCARDED_GEOMETRY_TYPES)] # Process individual geometries in MultiPolygon and GeometryCollection data types diff --git a/cea/default.config b/cea/default.config index a592a5c738..40a0052e42 100644 --- a/cea/default.config +++ b/cea/default.config @@ -66,11 +66,6 @@ grid-size.type = IntegerParameter grid-size.nullable = false grid-size.help = Grid size for the terrain file (do not use less than 10 meters). -elevation = 1 -elevation.type = IntegerParameter -elevation.nullable = false -elevation.help = Fixed elevation of the terrain (in meters). - [surroundings-helper] buffer = 50 buffer.type = RealParameter @@ -159,13 +154,13 @@ walls-grid.help = Grid resolution for wall surfaces. Use 200 (maximum) if you wa walls-grid.category = Level of Detail zone-geometry = 2 -zone-geometry.type = IntegerParameter -zone-geometry.help = Simplification level of the zone geometry (1 is the lowest). +zone-geometry.type = RealParameter +zone-geometry.help = Simplification level of zone geometry; tolerance distance of points in meters from original. zone-geometry.category = Level of Detail surrounding-geometry = 5 -surrounding-geometry.type = IntegerParameter -surrounding-geometry.help = Simplification level of the surroundings geometry (1 is the lowest). +surrounding-geometry.type = RealParameter +surrounding-geometry.help = Simplification level of surroundings geometry; tolerance distance of points in meters from original. surrounding-geometry.category = Level of Detail consider-floors = true @@ -658,10 +653,10 @@ optimization.type = BooleanParameter optimization.help = True if executing CEA district-energy-supply-system-optimization. Change default parameters under the respective tabs (i.e. supply system part 1 - decentralized and supply system part 2 - centralized, beta version) and Save to Config. [test] -workflow = quick -workflow.type = ChoiceParameter -workflow.choices = quick, medium, slow, unittests -workflow.help = The test workflow to run +type = unittest +type.type = ChoiceParameter +type.choices = unittest, integration +type.help = The test workflow to run [trace-inputlocator] scripts = archetypes-mapper, demand, emissions diff --git a/cea/dev/contributors_template b/cea/dev/contributors_template index 2d8b0438e8..da7d21304b 100644 --- a/cea/dev/contributors_template +++ b/cea/dev/contributors_template @@ -1,36 +1,15 @@ Developers: - * [Amr Elesawy](https://cityenergyanalyst.com/people) - * [Jimeno A. Fonseca](https://cityenergyanalyst.com/people) - * [Gabriel Happle](https://cityenergyanalyst.com/people) - * [Shanshan Hsieh](https://cityenergyanalyst.com/people) - * [Reynold Mok](https://cityenergyanalyst.com/people) - * [Martín Mosteiro Romero](https://cityenergyanalyst.com/people) - * [Mathias Niffeler](https://cityenergyanalyst.com/people) - * [Anastasiya Popova](https://cityenergyanalyst.com/people) - * [Zhongming Shi](https://cityenergyanalyst.com/people) - * [Luis Santos](https://cityenergyanalyst.com/people) - * [Bhargava Krishna Sreepathi](https://cityenergyanalyst.com/people) - * [Daren Thomas](https://cityenergyanalyst.com/people) + * Reynold Mok + * Martín Mosteiro Romero + * Mathias Niffeler + * Zhongming Shi Product lead: - * [Zhongming Shi](https://cityenergyanalyst.com/people) + * Zhongming Shi Product sponsor: * [Arno Schlueter](https://systems.arch.ethz.ch/arno-schlueter) * [Toni Piëch Foundation](https://www.tonipiechfoundation.org/) - Collaborators: - * Jose Bello - * Kian Wee Chen - * Jack Hawthorne - * Fazel Khayatian - * Victor Marty - * Rowan Molony - * Paul Neitzel - * Thuy-An Nguyen - * Bo Lie Ong - * Emanuel Riegelbauer - * Lennart Rogenhofer - * Toivo Säwén - * Sebastian Troiztsch - * Tim Vollrath + Others: + * https://cityenergyanalyst.com/people diff --git a/cea/inputlocator.py b/cea/inputlocator.py index a0adbe7018..aea1a25553 100644 --- a/cea/inputlocator.py +++ b/cea/inputlocator.py @@ -1,7 +1,7 @@ """ inputlocator.py - locate input files by name based on the reference folder structure. """ - +import atexit import os import cea.schemas import shutil @@ -34,10 +34,8 @@ def __init__(self, scenario, plugins=None): self._wrap_locator_methods(plugins) self.plugins = plugins - self._temp_directory = tempfile.TemporaryDirectory() - - def __del__(self): - self._temp_directory.cleanup() + self._temp_directory = tempfile.mkdtemp() + atexit.register(self._cleanup_temp_directory) def __getstate__(self): """Make sure we can pickle an InputLocator...""" @@ -59,6 +57,14 @@ def __setstate__(self, state): self._wrap_locator_methods(self.plugins) self._temp_directory = state["_temp_directory"] + def __del__(self): + self._cleanup_temp_directory() + + def _cleanup_temp_directory(self): + # Cleanup the temporary directory when the object is destroyed + if os.path.exists(self._temp_directory): + shutil.rmtree(self._temp_directory) + def _wrap_locator_methods(self, plugins): """ For each locator method defined in schemas.yml, wrap it in a callable object (preserving the @@ -1151,7 +1157,7 @@ def get_timeseries_plots_file(self, building, category=''): # OTHER def get_temporary_folder(self): """Temporary folder as returned by `tempfile`.""" - return self._temp_directory.name + return self._temp_directory def get_temporary_file(self, filename): """Returns the path to a file in the temporary folder with the name `filename`""" diff --git a/cea/interfaces/dashboard/api/inputs.py b/cea/interfaces/dashboard/api/inputs.py index 66be0bd320..1c640ce586 100644 --- a/cea/interfaces/dashboard/api/inputs.py +++ b/cea/interfaces/dashboard/api/inputs.py @@ -263,6 +263,8 @@ def get_building_properties(): columns[column_name]['regex'] = column['regex'] if 'example' in column: columns[column_name]['example'] = column['example'] + if 'nullable' in column: + columns[column_name]['nullable'] = column['nullable'] columns[column_name]['description'] = column["description"] columns[column_name]['unit'] = column["unit"] store['columns'][db] = columns diff --git a/cea/interfaces/dashboard/server/jobs.py b/cea/interfaces/dashboard/server/jobs.py index 6b0c5cf987..e592a573de 100644 --- a/cea/interfaces/dashboard/server/jobs.py +++ b/cea/interfaces/dashboard/server/jobs.py @@ -101,6 +101,7 @@ def post(self): @api.route("/list") class ListJobs(Resource): + @api.marshal_with(job_info_model, as_list=True) def get(self): return [asdict(job) for job in jobs.values()] diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index be361406f9..2df53a6f91 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -13,6 +13,7 @@ from itertools import repeat import numpy as np +import pandas as pd import py4design.py3dmodel.calculate as calculate import py4design.py3dmodel.construct as construct import py4design.py3dmodel.fetch as fetch @@ -20,7 +21,7 @@ import py4design.py3dmodel.utility as utility from OCC.Core.IntCurvesFace import IntCurvesFace_ShapeIntersector from OCC.Core.gp import gp_Pnt, gp_Lin, gp_Ax1, gp_Dir -from osgeo import osr +from osgeo import osr, gdal from py4design import urbangeom import cea @@ -37,6 +38,9 @@ __email__ = "cea@arch.ethz.ch" __status__ = "Production" +from cea.utilities.standardize_coordinates import (get_lat_lon_projected_shapefile, get_projected_coordinate_system, + crs_to_epsg) + def identify_surfaces_type(occface_list): roof_list = [] @@ -73,17 +77,18 @@ def identify_surfaces_type(occface_list): return facade_list_north, facade_list_west, facade_list_east, facade_list_south, roof_list, footprint_list -def calc_intersection(terrain_intersection_curves, edges_coords, edges_dir): +def calc_intersection(surface, edges_coords, edges_dir, tolerance): """ This script calculates the intersection of the building edges to the terrain, - :param terrain_intersection_curves: - :param edges_coords: - :param edges_dir: - :return: intersecting points, intersecting faces """ - building_line = gp_Lin(gp_Ax1(gp_Pnt(edges_coords[0], edges_coords[1], edges_coords[2]), - gp_Dir(edges_dir[0], edges_dir[1], edges_dir[2]))) - terrain_intersection_curves.PerformNearest(building_line, 0.0, float("+inf")) + point = gp_Pnt(edges_coords[0], edges_coords[1], edges_coords[2]) + direction = gp_Dir(edges_dir[0], edges_dir[1], edges_dir[2]) + line = gp_Lin(gp_Ax1(point, direction)) + + terrain_intersection_curves = IntCurvesFace_ShapeIntersector() + terrain_intersection_curves.Load(surface, tolerance) + terrain_intersection_curves.PerformNearest(line, float("-inf"), float("+inf")) + if terrain_intersection_curves.IsDone(): npts = terrain_intersection_curves.NbPnt() if npts != 0: @@ -114,8 +119,7 @@ def calc_building_solids(buildings_df, geometry_simplification, elevation_map, n nfloor_col_name = "floors_ag" # simplify geometry for buildings of interest - geometries = buildings_df.geometry.map( - lambda geometry: geometry.simplify(geometry_simplification, preserve_topology=True)) + geometries = buildings_df.geometry.simplify(geometry_simplification, preserve_topology=True) height = buildings_df[height_col_name].astype(float) nfloors = buildings_df[nfloor_col_name].astype(int) @@ -139,7 +143,7 @@ def calc_floor_to_floor_height(building_height, number_of_floors): def process_geometries(geometry, elevation_map, range_floors, floor_to_floor_height): elevation_map_for_geometry = elevation_map.get_elevation_map_from_geometry(geometry) # burn buildings footprint into the terrain and return the location of the new face - face_footprint = burn_buildings(geometry, elevation_map_for_geometry) + face_footprint = burn_buildings(geometry, elevation_map_for_geometry, 1e-12) # create floors and form a solid building_solid = calc_solid(face_footprint, range_floors, floor_to_floor_height) @@ -357,7 +361,7 @@ def calc_building_geometry_zone(name, building_solid, all_building_solid_list, a return name -def burn_buildings(geometry, elevation_map): +def burn_buildings(geometry, elevation_map, tolerance): if geometry.has_z: # remove elevation - we'll add it back later by intersecting with the topography point_list_2D = ((a, b) for (a, b, _) in geometry.exterior.coords) @@ -370,14 +374,12 @@ def burn_buildings(geometry, elevation_map): # get the midpt of the face face_midpt = calculate.face_midpt(face) - terrain_tin = elevation_map.generate_tin() + terrain_tin = elevation_map.generate_tin(tolerance) # make shell out of tin_occface_list and create OCC object terrain_shell = construct.make_shell(terrain_tin) - terrain_intersection_curves = IntCurvesFace_ShapeIntersector() - terrain_intersection_curves.Load(terrain_shell, 1e-6) # project the face_midpt to the terrain and get the elevation - inter_pt, inter_face = calc_intersection(terrain_intersection_curves, face_midpt, (0, 0, 1)) + inter_pt, inter_face = calc_intersection(terrain_shell, face_midpt, (0, 0, 1), tolerance) # reconstruct the footprint with the elevation loc_pt = (inter_pt.X(), inter_pt.Y(), inter_pt.Z()) @@ -491,82 +493,100 @@ def calc_intersection_face_solid(potentially_intersecting_solid, point): class ElevationMap(object): - __slots__ = ['elevation_map', 'x_coords', 'y_coords'] + __slots__ = ['elevation_map', 'x_coords', 'y_coords', 'x_size', 'y_size', 'nodata'] - def __init__(self, elevation_map, x_coords, y_coords): + def __init__(self, elevation_map, x_coords, y_coords, x_size, y_size, nodata=None): self.elevation_map = elevation_map self.x_coords = x_coords self.y_coords = y_coords + self.x_size = x_size + self.y_size = y_size + + self.nodata = nodata + @classmethod - def read_raster(cls, raster, raise_above_sea_level=True): + def read_raster(cls, raster): band = raster.GetRasterBand(1) + nodata = band.GetNoDataValue() + a = band.ReadAsArray() - if raise_above_sea_level and (a < 0).any(): - print('Warning: Some heights are below sea level') - # if height is below sea level, the entire case study is lifted to the lowest point is at altitude 0 - print('Adjusting elevation map to above sea level') - a -= a.min() - - (y, x) = np.shape(a) - (upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = raster.GetGeoTransform() + + y, x = np.shape(a) + upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size = raster.GetGeoTransform() + + if x_rotation != 0 or y_rotation != 0: + raise ValueError("Rotation in raster is not supported.") + x_coords = np.arange(start=0, stop=x) * x_size + upper_left_x + (x_size / 2) # add half the cell size y_coords = np.arange(start=0, stop=y) * y_size + upper_left_y + (y_size / 2) # to centre the point - return cls(a, x_coords, y_coords) + return cls(a, x_coords, y_coords, x_size, y_size, nodata) - def get_elevation_map_from_geometry(self, geometry, extra_points=5): + def get_elevation_map_from_geometry(self, geometry, extra_points=3): minx, miny, maxx, maxy = geometry.bounds - x_start = np.where(minx > self.x_coords)[0] - x_end = np.where(maxx < self.x_coords)[0] - y_start = np.where(maxy < self.y_coords)[0] - y_end = np.where(miny > self.y_coords)[0] + # Ensure geometry bounds is within elevation map + if (minx < self.x_coords[0] - self.x_size or maxx > self.x_coords[-1] + self.x_size + or miny < self.y_coords[-1] + self.y_size or maxy > self.y_coords[0] - self.y_size): + raise ValueError("Geometry bounds in not within the elevation map.") + + x_start = np.searchsorted(self.x_coords - self.x_size, minx, side='left') - 1 + x_end = np.searchsorted(self.x_coords + self.x_size, maxx, side='right') + y_start = len(self.y_coords) - np.searchsorted((self.y_coords - self.y_size)[::-1], maxy, side='left') - 1 + y_end = len(self.y_coords) - np.searchsorted((self.y_coords + self.y_size)[::-1], miny, side='right') - x_start = max(x_start[-1] - extra_points, 0) - x_end = min(x_end[0] + extra_points, len(self.x_coords)) - y_start = max(y_start[-1] - extra_points, 0) - y_end = min(y_end[0] + extra_points, len(self.y_coords)) + # Consider extra points + x_start = max(x_start - extra_points, 0) + x_end = min(x_end + extra_points, len(self.x_coords)) + y_start = max(y_start - extra_points, 0) + y_end = min(y_end + extra_points, len(self.y_coords)) new_elevation_map = self.elevation_map[y_start:y_end + 1, x_start:x_end + 1] new_x_coords = self.x_coords[x_start:x_end + 1] new_y_coords = self.y_coords[y_start:y_end + 1] - return ElevationMap(new_elevation_map, new_x_coords, new_y_coords) + return ElevationMap(new_elevation_map, new_x_coords, new_y_coords, self.x_size, self.y_size, self.nodata) - def generate_tin(self): - (y_index, x_index) = np.nonzero(self.elevation_map >= 0) + def generate_tin(self, tolerance=1e-6): + # Ignore no data values from raster + y_index, x_index = np.nonzero(self.elevation_map != self.nodata) _x_coords = self.x_coords[x_index] _y_coords = self.y_coords[y_index] - raster_points = [(x, y, z) for x, y, z in zip(_x_coords, _y_coords, self.elevation_map[y_index, x_index])] + raster_points = ((x, y, z) for x, y, z in zip(_x_coords, _y_coords, self.elevation_map[y_index, x_index])) - tin_occface_list = construct.delaunay3d(raster_points) + tin_occface_list = construct.delaunay3d(raster_points, tolerance=tolerance) return tin_occface_list def standardize_coordinate_systems(zone_df, surroundings_df, terrain_raster): - # Get projection of terrain and apply to zone and surroundings - terrian_projection = terrain_raster.GetProjection() - proj4_str = osr.SpatialReference(wkt=terrian_projection).ExportToProj4() - zone_df = zone_df.to_crs(proj4_str) - surroundings_df = surroundings_df.to_crs(proj4_str) + # Change all to projected cr (to meters) + lat, lon = get_lat_lon_projected_shapefile(zone_df) + crs = get_projected_coordinate_system(lat, lon) - return zone_df, surroundings_df, terrain_raster + reprojected_zone_df = zone_df.to_crs(crs) + reprojected_surroundings_df = surroundings_df.to_crs(crs) + + reprojected_terrain = gdal.Warp( + '', # Empty string as the output file path means in-memory + terrain_raster, + format='VRT', # Use VRT format for in-memory operation + dstSRS=crs + ) + + print(f"Reprojected scene to `EPSG:{crs_to_epsg(crs)}`") + return reprojected_zone_df, reprojected_surroundings_df, reprojected_terrain def check_terrain_bounds(zone_df, surroundings_df, terrain_raster): + total_df = pd.concat([zone_df, surroundings_df]) + # minx, miny, maxx, maxy - zone_bounds = zone_df.geometry.total_bounds - if len(surroundings_df): - surroundings_bounds = surroundings_df.geometry.total_bounds - else: # set bounds to zone if no surroundings - surroundings_bounds = zone_bounds - geometry_bounds = (min(zone_bounds[0], surroundings_bounds[0]), min(zone_bounds[1], surroundings_bounds[1]), - max(zone_bounds[2], surroundings_bounds[2]), max(zone_bounds[3], surroundings_bounds[3])) - - (upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = terrain_raster.GetGeoTransform() + geometry_bounds = total_df.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 diff --git a/cea/schemas.py b/cea/schemas.py index 020c442e17..4aebf0c353 100644 --- a/cea/schemas.py +++ b/cea/schemas.py @@ -5,7 +5,7 @@ """ import os -import pickle +from typing import List, Optional, Dict import pandas as pd import yaml @@ -28,41 +28,23 @@ __schemas = {} -def schemas(plugins): +def schemas(plugins: Optional[List] = None) -> Dict: """Return the contents of the schemas.yml file :parameter plugins: the list of plugins to generate the schemas for. Use ``config.plugins`` for this. :type plugins: List[cea.plugin.CeaPlugin] """ + if plugins is None: + plugins = [] + # loading schemas.yml is quite expensive - try to avoid it as much as possible by # maintaining a cache of the schemas - using the plugin list as a key global __schemas key = ":".join(str(p) for p in plugins) if key not in __schemas: - # load schemas.yml from disk - here again we use a cache: a pickled version is stored in the user folder. schemas_yml = os.path.join(os.path.dirname(__file__), 'schemas.yml') - schemas_pickle = os.path.expanduser("~/schemas.yml.pickle") - - def load_schemas_dict_from_yaml(): - schemas_dict = yaml.load(open(schemas_yml, "rb"), Loader=yaml.CLoader) - # ... so write out a pickle for next time - with open(schemas_pickle, "wb") as schemas_pickle_fp: - pickle.dump(schemas_dict, schemas_pickle_fp) - return schemas_dict - - # compare the dates of the two files - use the pickle if it's newer - schemas_dict = None - if os.path.exists(schemas_pickle) and os.path.getmtime(schemas_pickle) > os.path.getmtime(schemas_yml): - with open(schemas_pickle, "r") as schemas_pickle_fp: - try: - schemas_dict = pickle.load(schemas_pickle_fp) - except Exception: - schemas_dict = None - - if not schemas_dict: - # ok. this will take a while to read... - schemas_dict = load_schemas_dict_from_yaml() - + with open(schemas_yml, "r") as f: + schemas_dict = yaml.load(f, Loader=yaml.CLoader) __schemas[key] = schemas_dict # add the plugins - these don't use caches as their schemas.yml are (probably) much shorter diff --git a/cea/schemas.yml b/cea/schemas.yml index 1e05001e0d..a02f3572c5 100644 --- a/cea/schemas.yml +++ b/cea/schemas.yml @@ -2619,8 +2619,7 @@ get_database_air_conditioning_systems: system type: float unit: '[C]' - values: '{0.0...n}' - min: 0.0 + values: '{n...n}' HOT_WATER: columns: Description: @@ -16084,13 +16083,13 @@ get_zone_geometry: values: alphanumeric nullable: true house_name: - description: House number (if any) + description: House name (if any) type: string unit: 'NA' values: alphanumeric nullable: true resi_type: - description: Residential type, e.g. HDB for Singapore (if any) + description: Residential type, if Singapore's HDB or not (if any) type: string unit: 'NA' values: alphanumeric @@ -16102,7 +16101,7 @@ get_zone_geometry: values: alphanumeric nullable: true country: - description: Country 2-digit code in the address (if any) + description: Country (2-digit) code in the address (if any) type: string unit: 'NA' values: alphanumeric diff --git a/cea/scripts.yml b/cea/scripts.yml index bcf301e22e..4adf4fe44c 100644 --- a/cea/scripts.yml +++ b/cea/scripts.yml @@ -28,7 +28,7 @@ Data management: - name: terrain-helper label: Terrain helper - description: Query topography with a fixed elevation + description: Query topography data from third party sources interfaces: [cli, dashboard] module: cea.datamanagement.terrain_helper parameters: ['general:scenario', terrain-helper] diff --git a/cea/technologies/solar/photovoltaic.py b/cea/technologies/solar/photovoltaic.py index 46314d0f29..d2396adfd5 100644 --- a/cea/technologies/solar/photovoltaic.py +++ b/cea/technologies/solar/photovoltaic.py @@ -6,7 +6,7 @@ import time from itertools import repeat from math import radians, degrees, asin, sin, acos, cos, exp, tan, atan, ceil, log -from multiprocessing import Pool +from multiprocessing.dummy import Pool import numpy as np import pandas as pd @@ -761,21 +761,22 @@ def aggregate_results_func(args): return aggregate_results(args[0], args[1]) -def write_aggregate_results(config, locator, building_names, num_process=1): +def write_aggregate_results(config, locator, building_names): aggregated_hourly_results_df = pd.DataFrame() aggregated_annual_results = pd.DataFrame() panel_type = config.solar.type_PVpanel - pool = Pool(processes=num_process) - args = [(locator, x) for x in np.array_split(building_names, num_process) if x.size != 0] - for i, x in enumerate(pool.map(aggregate_results_func, args)): - hourly_results_df, annual_results = x - if i == 0: - aggregated_hourly_results_df = hourly_results_df - aggregated_annual_results = annual_results - else: - aggregated_hourly_results_df = aggregated_hourly_results_df + hourly_results_df - aggregated_annual_results = pd.concat([aggregated_annual_results, annual_results], axis=1, sort=False) + num_process = 4 + with Pool(processes=num_process) as pool: + args = [(locator, x) for x in np.array_split(building_names, num_process) if x.size != 0] + for i, x in enumerate(pool.map(aggregate_results_func, args)): + hourly_results_df, annual_results = x + if i == 0: + aggregated_hourly_results_df = hourly_results_df + aggregated_annual_results = annual_results + else: + aggregated_hourly_results_df = aggregated_hourly_results_df + hourly_results_df + aggregated_annual_results = pd.concat([aggregated_annual_results, annual_results], axis=1, sort=False) # save hourly results aggregated_hourly_results_df.to_csv(locator.PV_totals(panel_type=panel_type), index=True, float_format='%.2f', na_rep='nan') @@ -821,7 +822,7 @@ def main(config): building_names) # aggregate results from all buildings - write_aggregate_results(config, locator, building_names, num_process) + write_aggregate_results(config, locator, building_names) if __name__ == '__main__': diff --git a/cea/tests/__init__.py b/cea/tests/__init__.py index 1d69ef7af9..7f0b9928ad 100644 --- a/cea/tests/__init__.py +++ b/cea/tests/__init__.py @@ -2,16 +2,11 @@ Run the CEA scripts and unit tests as part of our CI efforts (cf. The Jenkins) """ - - - - import os -import shutil -import tempfile +import unittest + import cea.config import cea.inputlocator -import cea.workflows.workflow __author__ = "Daren Thomas" __copyright__ = "Copyright 2020, Architecture and Building Systems - ETH Zurich" @@ -22,22 +17,25 @@ __email__ = "cea@arch.ethz.ch" __status__ = "Production" +from cea.tests.test_workflow import TestWorkflows + def main(config): - workflow_yml = os.path.join(os.path.dirname(__file__), "workflow_{workflow}.yml".format(workflow=config.test.workflow)) + test_type = config.test.type + + if test_type == "unittest": + test_suite = unittest.defaultTestLoader.discover(os.path.dirname(__file__)) + result = unittest.TextTestRunner().run(test_suite) - default_config = cea.config.Configuration(cea.config.DEFAULT_CONFIG) - default_config.project = os.path.join(tempfile.gettempdir(), "reference-case-open") - default_config.workflow.workflow = workflow_yml - default_config.workflow.resume = False - default_config.workflow.resume_file = os.path.join(tempfile.gettempdir(), "resume.yml") # force resume file to be temporary + if not result.wasSuccessful(): + raise AssertionError("Unittests failed.") - if os.path.exists(default_config.project): - # make sure we're working on a clean slate - shutil.rmtree(default_config.project) + elif test_type == "integration": + TestWorkflows()._test_workflows() - cea.workflows.workflow.main(default_config) + else: + raise Exception(f"Test type '{test_type}' not supported") if __name__ == '__main__': - main(cea.config.Configuration()) \ No newline at end of file + main(cea.config.Configuration()) diff --git a/cea/tests/run_unit_tests.py b/cea/tests/run_unit_tests.py deleted file mode 100644 index 8789fc2a94..0000000000 --- a/cea/tests/run_unit_tests.py +++ /dev/null @@ -1,33 +0,0 @@ -""" -Run all the unit tests in the cea/tests folder -""" - - - - - -import os -import unittest -import cea.config -import cea.workflows.workflow - -__author__ = "Daren Thomas" -__copyright__ = "Copyright 2020, Architecture and Building Systems - ETH Zurich" -__credits__ = ["Daren Thomas"] -__license__ = "MIT" -__version__ = "0.1" -__maintainer__ = "Daren Thomas" -__email__ = "cea@arch.ethz.ch" -__status__ = "Production" - - -def main(_): - test_suite = unittest.defaultTestLoader.discover(os.path.dirname(__file__)) - result = unittest.TextTestRunner(verbosity=1).run(test_suite) - - if not result.wasSuccessful(): - raise AssertionError("Unittests failed.") - - -if __name__ == "__main__": - main(cea.config.Configuration) \ No newline at end of file diff --git a/cea/tests/test_workflow.py b/cea/tests/test_workflow.py new file mode 100644 index 0000000000..55b9a28f36 --- /dev/null +++ b/cea/tests/test_workflow.py @@ -0,0 +1,32 @@ +import glob +import multiprocessing.pool +import os +import unittest +from typing import List + +import cea.config +from cea.workflows.workflow import main + + +class TestWorkflows(unittest.TestCase): + @staticmethod + def get_test_workflows() -> List[str]: + dirname = os.path.join(os.path.realpath(os.path.dirname(__file__)), "workflows") + workflows = [workflow for workflow in glob.glob(os.path.join(dirname, "*.yml"))] + + return [os.path.join(dirname, workflow) for workflow in workflows] + + @staticmethod + def get_workflow_config(workflow: str) -> cea.config.Configuration: + config = cea.config.Configuration() + config.workflow.workflow = workflow + + return config + + def _test_workflows(self): + with multiprocessing.pool.Pool() as p: + p.map(main, [self.get_workflow_config(workflow) for workflow in self.get_test_workflows()]) + + +if __name__ == '__main__': + unittest.main() diff --git a/cea/tests/workflow_medium.yml b/cea/tests/workflow_medium.yml deleted file mode 100644 index a6ce3282b2..0000000000 --- a/cea/tests/workflow_medium.yml +++ /dev/null @@ -1,138 +0,0 @@ -- config: . - "general:multiprocessing": False - "general:project": "${CEA_general_project}/../reference-case-open" - "radiation:daysim-bin-directory": "${CEA_radiation_daysim-bin-directory}" -#- script: run-unit-tests -- script: extract-reference-case - parameters: - destination: "{general:project}/.." - case: open -- config: . - "general:scenario-name": baseline - -# HEATING CASE (CH) -- script: data-initializer - parameters: - databases-path: CH - databases: [archetypes, assemblies, components] -- script: weather-helper - parameters: - weather: Zug-inducity_1990_2010_TMY -- script: archetypes-mapper - parameters: - input-databases: ['comfort', 'architecture', 'air-conditioning', 'internal-loads', 'supply', 'schedules'] - buildings: [] -- script: radiation - parameters: - neglect-adjacent-buildings: false -- script: schedule-maker -- script: demand -- script: emissions -- script: system-costs -- script: water-body-potential -- script: sewage-potential -- script: shallow-geothermal-potential -- script: photovoltaic -- script: solar-collector - parameters: - type-scpanel: FP -- script: solar-collector - parameters: - type-scpanel: ET -- script: photovoltaic-thermal - parameters: - type-scpanel: FP -- script: photovoltaic-thermal - parameters: - type-scpanel: ET -- script: network-layout - parameters: - network-type: DH -- script: thermal-network - parameters: - network-type: DH - network-model: detailed -# stop-t: 744 # run for one month - start-t: 0 - stop-t: 24 # run for one day -- script: thermal-network - parameters: - network-type: DH - network-model: simplified -- script: decentralized -- script: optimization-new - parameters: - network-type: DH - ga-number-of-generations: 2 - ga-population-size: 5 -#- script: multi-criteria-analysis -# parameters: -# generation: 2 -#- script: run-all-plots -# parameters: -# network-type: DH -# network-name: "" - -# COOLING CASE (SG) -- script: data-initializer - parameters: - databases-path: SG - databases: [archetypes, assemblies, components] -- script: archetypes-mapper - parameters: - input-databases: [comfort, architecture, air-conditioning, internal-loads, supply, schedules] -- script: weather-helper - parameters: - weather: Singapore-Changi_1990_2010_TMY -- script: radiation - parameters: - neglect-adjacent-buildings: false -- script: schedule-maker -- script: demand -- script: emissions -- script: system-costs -- script: water-body-potential -- script: sewage-potential -- script: shallow-geothermal-potential -- script: photovoltaic -- script: solar-collector - parameters: - type-scpanel: FP -- script: solar-collector - parameters: - type-scpanel: ET -- script: photovoltaic-thermal - parameters: - type-scpanel: FP -- script: photovoltaic-thermal - parameters: - type-scpanel: ET -- script: network-layout - parameters: - network-type: DC - consider-only-buildings-with-demand: on -- script: thermal-network - parameters: - network-type: DC - network-model: detailed - #stop-t: 744 # run for one month - start-t: 0 - stop-t: 24 # run for one day -- script: thermal-network - parameters: - network-type: DC - network-model: simplified -- script: decentralized -#- script: optimization-new -# parameters: -# network-type: DC -# ga-number-of-generations: 2 -# ga-population-size: 5 -#- script: multi-criteria-analysis -# parameters: -# generation: 2 -#- script: run-all-plots -# parameters: -# plant-node: NODE41 -# network-type: DC -# network-name: "" diff --git a/cea/tests/workflow_quick.yml b/cea/tests/workflow_quick.yml deleted file mode 100644 index 75832b90af..0000000000 --- a/cea/tests/workflow_quick.yml +++ /dev/null @@ -1,28 +0,0 @@ -- config: default - "general:multiprocessing": False - "general:project": "${CEA_general_project}/../reference-case-open" - "radiation:daysim-bin-directory": "${CEA_radiation_daysim-bin-directory}" -- script: run-unit-tests -- script: extract-reference-case - parameters: - destination: "{general:project}/.." - case: open -- config: . - "general:scenario-name": baseline -- script: data-initializer - parameters: - databases-path: CH - databases: [archetypes, assemblies, components] -- script: archetypes-mapper - parameters: - input-databases: [comfort, architecture, air-conditioning, internal-loads, supply, schedules] -- script: weather-helper - parameters: - weather: Zug_inducity_2009 -- script: schedule-maker -- script: demand -- script: emissions -- script: system-costs -- script: water-body-potential -- script: sewage-potential -- script: shallow-geothermal-potential \ No newline at end of file diff --git a/cea/tests/workflow_unittests.yml b/cea/tests/workflow_unittests.yml deleted file mode 100644 index f2c3300729..0000000000 --- a/cea/tests/workflow_unittests.yml +++ /dev/null @@ -1,4 +0,0 @@ -- config: . - "general:multiprocessing": on - "radiation:daysim-bin-directory": "${CEA_radiation_daysim-bin-directory}" -- script: run-unit-tests \ No newline at end of file diff --git a/cea/tests/workflow_slow.yml b/cea/tests/workflows/sg_cooling.yml similarity index 52% rename from cea/tests/workflow_slow.yml rename to cea/tests/workflows/sg_cooling.yml index dcddcc0ca5..dd845f7b40 100644 --- a/cea/tests/workflow_slow.yml +++ b/cea/tests/workflows/sg_cooling.yml @@ -1,8 +1,7 @@ - config: default "general:multiprocessing": False - "general:project": "${CEA_general_project}/../reference-case-open" + "general:project": "${CEA_general_project}/../sg_cooling/reference-case-open" "radiation:daysim-bin-directory": "${CEA_radiation_daysim-bin-directory}" -- script: run-unit-tests - script: extract-reference-case parameters: destination: "{general:project}/.." @@ -10,10 +9,9 @@ - config: . "general:scenario-name": baseline -# HEATING CASE (CH) - script: data-initializer parameters: - databases-path: CH + databases-path: SG databases: [archetypes, assemblies, components] - script: archetypes-mapper parameters: @@ -21,69 +19,6 @@ - script: surroundings-helper - script: streets-helper - script: terrain-helper -- script: weather-helper - parameters: - weather: Zug-inducity_1990_2010_TMY -- script: radiation - parameters: - neglect-adjacent-buildings: false -- script: schedule-maker -- script: demand -- script: emissions -- script: system-costs -- script: water-body-potential -- script: sewage-potential -- script: shallow-geothermal-potential -- script: photovoltaic -- script: solar-collector - parameters: - type-scpanel: FP -- script: solar-collector - parameters: - type-scpanel: ET -- script: photovoltaic-thermal - parameters: - type-scpanel: FP -- script: photovoltaic-thermal - parameters: - type-scpanel: ET -- script: network-layout - parameters: - network-type: DH -- script: thermal-network - parameters: - network-type: DH - network-model: detailed - stop-t: 744 # run for one month -- script: thermal-network - parameters: - network-type: DH - network-model: simplified -- script: decentralized -- script: optimization-new - parameters: - network-type: DH - ga-number-of-generations: 2 - ga-population-size: 5 -#- script: multi-criteria-analysis -# parameters: -# generation: 2 -#- script: run-all-plots -# parameters: -# network-type: DH -# network-name: "" - -## COOLING CASE (SG) -- script: data-initializer - parameters: - databases-path: SG - databases: [archetypes, assemblies, components] -- script: archetypes-mapper - parameters: - input-databases: [comfort, architecture, air-conditioning, internal-loads, supply, schedules] -#- script: surroundings-helper -#- script: streets-helper -#- script: terrain-helper - script: weather-helper parameters: weather: Singapore-Changi_1990_2010_TMY @@ -124,11 +59,11 @@ network-type: DC network-model: simplified - script: decentralized -#- script: optimization-new -# parameters: -# network-type: DC -# ga-number-of-generations: 2 -# ga-population-size: 5 +- script: optimization-new + parameters: + network-type: DC + ga-number-of-generations: 2 + ga-population-size: 5 #- script: multi-criteria-analysis # parameters: # generation: 2 diff --git a/cea/tests/workflows/zug_heating.yml b/cea/tests/workflows/zug_heating.yml new file mode 100644 index 0000000000..549fb757a9 --- /dev/null +++ b/cea/tests/workflows/zug_heating.yml @@ -0,0 +1,72 @@ +- config: default + "general:multiprocessing": False + "general:project": "${CEA_general_project}/../zug_heating/reference-case-open" + "radiation:daysim-bin-directory": "${CEA_radiation_daysim-bin-directory}" +- script: extract-reference-case + parameters: + destination: "{general:project}/.." + case: open +- config: . + "general:scenario-name": baseline + +- script: data-initializer + parameters: + databases-path: CH + databases: [archetypes, assemblies, components] +- script: archetypes-mapper + parameters: + input-databases: [comfort, architecture, air-conditioning, internal-loads, supply, schedules] +- script: surroundings-helper +- script: streets-helper +- script: terrain-helper +- script: weather-helper + parameters: + weather: Zug-inducity_1990_2010_TMY +- script: radiation + parameters: + neglect-adjacent-buildings: false +- script: schedule-maker +- script: demand +- script: emissions +- script: system-costs +- script: water-body-potential +- script: sewage-potential +- script: shallow-geothermal-potential +- script: photovoltaic +- script: solar-collector + parameters: + type-scpanel: FP +- script: solar-collector + parameters: + type-scpanel: ET +- script: photovoltaic-thermal + parameters: + type-scpanel: FP +- script: photovoltaic-thermal + parameters: + type-scpanel: ET +- script: network-layout + parameters: + network-type: DH +- script: thermal-network + parameters: + network-type: DH + network-model: detailed + stop-t: 744 # run for one month +- script: thermal-network + parameters: + network-type: DH + network-model: simplified +- script: decentralized +- script: optimization-new + parameters: + network-type: DH + ga-number-of-generations: 2 + ga-population-size: 5 +#- script: multi-criteria-analysis +# parameters: +# generation: 2 +#- script: run-all-plots +# parameters: +# network-type: DH +# network-name: "" diff --git a/cea/utilities/standardize_coordinates.py b/cea/utilities/standardize_coordinates.py index 721bd3f5a4..eb2d42dbc0 100644 --- a/cea/utilities/standardize_coordinates.py +++ b/cea/utilities/standardize_coordinates.py @@ -12,6 +12,8 @@ __email__ = "cea@arch.ethz.ch" __status__ = "Production" +from pyproj import CRS + def shapefile_to_WSG_and_UTM(shapefile_path): @@ -48,17 +50,18 @@ def raster_to_WSG_and_UTM(raster_path, lat, lon): def get_geographic_coordinate_system(): - return "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" + return CRS.from_epsg(4326).to_wkt() def get_projected_coordinate_system(lat, lon): easting, northing, zone_number, zone_letter = utm.from_latlon(lat, lon) - if zone_letter in "NPQRSTUVWXX": - return "+proj=utm +zone=" + str(zone_number) + " +ellps=WGS84 +datum=WGS84 +units=m +no_defs" - elif zone_letter in "CDEFGHJKLM": - return "+proj=utm +zone=" + str(zone_number) + " +ellps=WGS84 +datum=WGS84 +units=m +no_defs +south" - else: - Exception('The projected coordinate system is unknown, lon{}, lat{}').format(lat, lon) + epsg = f"326{zone_number}" if lon >= 0 else f"327{zone_number}" + + return CRS.from_epsg(int(epsg)).to_wkt() + + +def crs_to_epsg(crs: str) -> int: + return CRS.from_string(crs).to_epsg() def get_lat_lon_projected_shapefile(data): diff --git a/cea/utilities/trace_inputlocator/__init__.py b/cea/utilities/trace_inputlocator/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/cea/tests/trace_inputlocator.output.gv b/cea/utilities/trace_inputlocator/trace_inputlocator.output.gv similarity index 100% rename from cea/tests/trace_inputlocator.output.gv rename to cea/utilities/trace_inputlocator/trace_inputlocator.output.gv diff --git a/cea/tests/trace_inputlocator.output.new.yml b/cea/utilities/trace_inputlocator/trace_inputlocator.output.new.yml similarity index 100% rename from cea/tests/trace_inputlocator.output.new.yml rename to cea/utilities/trace_inputlocator/trace_inputlocator.output.new.yml diff --git a/cea/tests/trace_inputlocator.py b/cea/utilities/trace_inputlocator/trace_inputlocator.py similarity index 100% rename from cea/tests/trace_inputlocator.py rename to cea/utilities/trace_inputlocator/trace_inputlocator.py diff --git a/cea/tests/trace_inputlocator.template.gv b/cea/utilities/trace_inputlocator/trace_inputlocator.template.gv similarity index 100% rename from cea/tests/trace_inputlocator.template.gv rename to cea/utilities/trace_inputlocator/trace_inputlocator.template.gv diff --git a/cea/workflows/workflow.py b/cea/workflows/workflow.py index c94801995b..4c5386fced 100644 --- a/cea/workflows/workflow.py +++ b/cea/workflows/workflow.py @@ -32,7 +32,7 @@ def run(config, script, **kwargs): def run_with_trace(config, script, **kwargs): """Same as run, but use the trace-inputlocator functionality to capture InputLocator calls""" - from cea.tests.trace_inputlocator import create_trace_function, update_trace_data, meta_to_yaml + from cea.utilities.trace_inputlocator.trace_inputlocator import create_trace_function, update_trace_data, meta_to_yaml if "multiprocessing" in kwargs: # we can only trace single processes