From 2202df0d062bfb03bb0cad54daff1857026348f3 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Wed, 22 May 2024 18:36:39 +0200 Subject: [PATCH 01/42] Add nullable info for column schema --- cea/interfaces/dashboard/api/inputs.py | 2 ++ 1 file changed, 2 insertions(+) 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 From da0ab1a87c0f405bf095ca665ce1d812b8a98e76 Mon Sep 17 00:00:00 2001 From: Zhongming Shi Date: Wed, 22 May 2024 18:48:22 +0200 Subject: [PATCH 02/42] Update schemas.yml --- cea/schemas.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cea/schemas.yml b/cea/schemas.yml index 1e05001e0d..b36e848ccf 100644 --- a/cea/schemas.yml +++ b/cea/schemas.yml @@ -16084,13 +16084,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 +16102,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 From 93d2d98b6f47b255310c92b998883315314c0210 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Wed, 22 May 2024 23:35:28 +0200 Subject: [PATCH 03/42] Remove min value --- cea/schemas.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/cea/schemas.yml b/cea/schemas.yml index b36e848ccf..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: From b659a4db5a3e9b2c5b61ae8f6f3e797db7447839 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Wed, 29 May 2024 08:38:01 +0200 Subject: [PATCH 04/42] Create release --- CHANGELOG.md | 2 ++ CREDITS.md | 19 +++++++++++++++++++ cea/__init__.py | 2 +- cea/dev/contributors_template | 35 +++++++---------------------------- 4 files changed, 29 insertions(+), 29 deletions(-) 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..d65b0e8ed3 100644 --- a/CREDITS.md +++ b/CREDITS.md @@ -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/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/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 From d2e0a23841a4a2e58f2d60a377e2865f410fa803 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Wed, 29 May 2024 08:58:09 +0200 Subject: [PATCH 05/42] Fetch micromamba for both mac arch --- .github/workflows/setup_build.yml | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) 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: From 5bdb18e172dc44355a55306bde13305f0083d61d Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Wed, 29 May 2024 13:30:20 +0200 Subject: [PATCH 06/42] Update zenodo links --- CREDITS.md | 2 +- README.rst | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CREDITS.md b/CREDITS.md index d65b0e8ed3..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 ============ diff --git a/README.rst b/README.rst index eab8021677..621a907424 100644 --- a/README.rst +++ b/README.rst @@ -6,8 +6,8 @@ :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 +33,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 From 55b59480ff6688d2164e922eecfea7503abe6902 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Thu, 30 May 2024 13:35:17 +0200 Subject: [PATCH 07/42] Move trace_inputlocator out of test folder --- bin/create_trace_graphviz.py | 2 +- cea/utilities/trace_inputlocator/__init__.py | 0 .../trace_inputlocator}/trace_inputlocator.output.gv | 0 .../trace_inputlocator}/trace_inputlocator.output.new.yml | 0 .../trace_inputlocator}/trace_inputlocator.py | 0 .../trace_inputlocator}/trace_inputlocator.template.gv | 0 cea/workflows/workflow.py | 2 +- 7 files changed, 2 insertions(+), 2 deletions(-) create mode 100644 cea/utilities/trace_inputlocator/__init__.py rename cea/{tests => utilities/trace_inputlocator}/trace_inputlocator.output.gv (100%) rename cea/{tests => utilities/trace_inputlocator}/trace_inputlocator.output.new.yml (100%) rename cea/{tests => utilities/trace_inputlocator}/trace_inputlocator.py (100%) rename cea/{tests => utilities/trace_inputlocator}/trace_inputlocator.template.gv (100%) 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/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 From 7e2d4bed9fbbcd74110e77d1ce99bbb1250adf68 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 14:24:34 +0200 Subject: [PATCH 08/42] Remove invalid lines of code badge --- README.rst | 2 -- 1 file changed, 2 deletions(-) diff --git a/README.rst b/README.rst index 621a907424..061fe4b30d 100644 --- a/README.rst +++ b/README.rst @@ -4,8 +4,6 @@ :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.11382063.svg :target: https://doi.org/10.5281/zenodo.11382063 From 9469e31c469a88ea37100ed2fe452811ee660773 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 14:35:23 +0200 Subject: [PATCH 09/42] Separate tests into unit and integration tests --- cea/default.config | 8 +- cea/tests/__init__.py | 34 ++--- cea/tests/test_workflow.py | 32 ++++ cea/tests/workflow_medium.yml | 138 ------------------ cea/tests/workflow_quick.yml | 28 ---- cea/tests/workflow_unittests.yml | 4 - .../sg_cooling.yml} | 79 +--------- cea/tests/workflows/zug_heating.yml | 72 +++++++++ 8 files changed, 131 insertions(+), 264 deletions(-) create mode 100644 cea/tests/test_workflow.py delete mode 100644 cea/tests/workflow_medium.yml delete mode 100644 cea/tests/workflow_quick.yml delete mode 100644 cea/tests/workflow_unittests.yml rename cea/tests/{workflow_slow.yml => workflows/sg_cooling.yml} (52%) create mode 100644 cea/tests/workflows/zug_heating.yml diff --git a/cea/default.config b/cea/default.config index a592a5c738..7d7b54633a 100644 --- a/cea/default.config +++ b/cea/default.config @@ -658,10 +658,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 = unittests +type.type = ChoiceParameter +type.choices = unittests, integration +type.help = The test workflow to run [trace-inputlocator] scripts = archetypes-mapper, demand, emissions diff --git a/cea/tests/__init__.py b/cea/tests/__init__.py index 1d69ef7af9..709a6144da 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 == "unittests": + 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/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: "" From ab713813262f11d169cd04b0b8371dc48d3ccefc Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 14:38:21 +0200 Subject: [PATCH 10/42] Use threads instead of processes This is to prevent errors due testing, since python does not allow child processes to spawn more processes. --- cea/technologies/solar/photovoltaic.py | 27 +++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) 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__': From fd12a1df9fb409718b0e89c29d123244aed45f62 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 14:38:51 +0200 Subject: [PATCH 11/42] Ensure temporary directory is deleted at exit --- cea/inputlocator.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/cea/inputlocator.py b/cea/inputlocator.py index a0adbe7018..af1aca0bb0 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(lambda: shutil.rmtree(self._temp_directory)) def __getstate__(self): """Make sure we can pickle an InputLocator...""" @@ -1151,7 +1149,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`""" From dccaeab42bc38ed97c327fdcfaf00b9c30009e77 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 20:47:39 +0200 Subject: [PATCH 12/42] Use singular instead of plural --- cea/default.config | 4 ++-- cea/tests/__init__.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cea/default.config b/cea/default.config index 7d7b54633a..17a1f00927 100644 --- a/cea/default.config +++ b/cea/default.config @@ -658,9 +658,9 @@ 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] -type = unittests +type = unittest type.type = ChoiceParameter -type.choices = unittests, integration +type.choices = unittest, integration type.help = The test workflow to run [trace-inputlocator] diff --git a/cea/tests/__init__.py b/cea/tests/__init__.py index 709a6144da..7f0b9928ad 100644 --- a/cea/tests/__init__.py +++ b/cea/tests/__init__.py @@ -23,7 +23,7 @@ def main(config): test_type = config.test.type - if test_type == "unittests": + if test_type == "unittest": test_suite = unittest.defaultTestLoader.discover(os.path.dirname(__file__)) result = unittest.TextTestRunner().run(test_suite) From 494c607cad537c961a1ac1d3f0dcc60995e07260 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 20:54:34 +0200 Subject: [PATCH 13/42] Remove the need for pickling schemas --- cea/schemas.py | 31 ++++++------------------------- 1 file changed, 6 insertions(+), 25 deletions(-) diff --git a/cea/schemas.py b/cea/schemas.py index 020c442e17..bfd01a1ba8 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,22 @@ __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() - + schemas_dict = yaml.load(open(schemas_yml, "rb"), Loader=yaml.CLoader) __schemas[key] = schemas_dict # add the plugins - these don't use caches as their schemas.yml are (probably) much shorter From df773122375dfbf6a052e6d754c4593e763b8da0 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 1 Jun 2024 21:01:40 +0200 Subject: [PATCH 14/42] Use context manager to read schema file --- cea/schemas.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cea/schemas.py b/cea/schemas.py index bfd01a1ba8..4aebf0c353 100644 --- a/cea/schemas.py +++ b/cea/schemas.py @@ -43,7 +43,8 @@ def schemas(plugins: Optional[List] = None) -> Dict: if key not in __schemas: schemas_yml = os.path.join(os.path.dirname(__file__), 'schemas.yml') - schemas_dict = yaml.load(open(schemas_yml, "rb"), Loader=yaml.CLoader) + 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 From 5fa56a8f5dd612e483652e14ffb5b7e1ba926269 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 3 Jun 2024 14:58:25 +0200 Subject: [PATCH 15/42] Update cea test command --- .github/workflows/main.yml | 4 +++- .github/workflows/pull.request.yml | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) 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 From 50f24b72b3e730eb7a29a6b788cd40f931759a3c Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 3 Jun 2024 16:35:22 +0200 Subject: [PATCH 16/42] Remove unused test script --- cea/tests/run_unit_tests.py | 33 --------------------------------- 1 file changed, 33 deletions(-) delete mode 100644 cea/tests/run_unit_tests.py 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 From 7c9736ce9959fd53b91b3e64287c003d79e26e6f Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 3 Jun 2024 18:06:23 +0200 Subject: [PATCH 17/42] Fix some warnings --- cea/datamanagement/surroundings_helper.py | 13 ++++++++----- cea/datamanagement/zone_helper.py | 11 +++++------ 2 files changed, 13 insertions(+), 11 deletions(-) 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/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 From f514a516ad7663fdc811aa51e6ce12cb6981e38c Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 3 Jun 2024 18:06:58 +0200 Subject: [PATCH 18/42] Update coordinate system functions --- cea/utilities/standardize_coordinates.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/cea/utilities/standardize_coordinates.py b/cea/utilities/standardize_coordinates.py index 721bd3f5a4..37a0607ebd 100644 --- a/cea/utilities/standardize_coordinates.py +++ b/cea/utilities/standardize_coordinates.py @@ -48,17 +48,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 "EPSG:4326" 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" + + if lon >= 0: + hemisphere = 'N' else: - Exception('The projected coordinate system is unknown, lon{}, lat{}').format(lat, lon) + hemisphere = 'S' + + return f"EPSG:326{zone_number}" if hemisphere == 'N' else f"EPSG:327{zone_number}" def get_lat_lon_projected_shapefile(data): From 46efc6e4c7652f12bdf1f19b5f9269e6fa1f8dc6 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 3 Jun 2024 18:30:16 +0200 Subject: [PATCH 19/42] Use WKT for CRS string format --- cea/utilities/standardize_coordinates.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/cea/utilities/standardize_coordinates.py b/cea/utilities/standardize_coordinates.py index 37a0607ebd..6cb3089cf5 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,18 +50,14 @@ def raster_to_WSG_and_UTM(raster_path, lat, lon): def get_geographic_coordinate_system(): - return "EPSG:4326" + 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) + epsg = f"326{zone_number}" if lon >= 0 else f"327{zone_number}" - if lon >= 0: - hemisphere = 'N' - else: - hemisphere = 'S' - - return f"EPSG:326{zone_number}" if hemisphere == 'N' else f"EPSG:327{zone_number}" + return CRS.from_epsg(int(epsg)).to_wkt() def get_lat_lon_projected_shapefile(data): From 805be6a77d15cb33cb6069b0c946626acaf45b21 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Tue, 4 Jun 2024 20:27:45 +0200 Subject: [PATCH 20/42] Add api marshalling as list Fixes datetime not JSON serializable error --- cea/interfaces/dashboard/server/jobs.py | 1 + 1 file changed, 1 insertion(+) 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()] From c5115f69f09808093fb3a18a2f8790adcc6521ae Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Wed, 5 Jun 2024 11:29:38 +0200 Subject: [PATCH 21/42] Cleanup temporary directory when object is destroyed --- cea/inputlocator.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/cea/inputlocator.py b/cea/inputlocator.py index af1aca0bb0..aea1a25553 100644 --- a/cea/inputlocator.py +++ b/cea/inputlocator.py @@ -35,7 +35,7 @@ def __init__(self, scenario, plugins=None): self.plugins = plugins self._temp_directory = tempfile.mkdtemp() - atexit.register(lambda: shutil.rmtree(self._temp_directory)) + atexit.register(self._cleanup_temp_directory) def __getstate__(self): """Make sure we can pickle an InputLocator...""" @@ -57,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 From 84f5896460940acccdeed330b9b3f2e7d050ad14 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 01:21:44 +0200 Subject: [PATCH 22/42] Simplify getting building bounds --- cea/resources/radiation/geometry_generator.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index be361406f9..38e93268fa 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 @@ -557,16 +558,12 @@ def standardize_coordinate_systems(zone_df, surroundings_df, terrain_raster): 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 From 24e59f8f90c5a2a5cc5915211c75e114ac2fdbac Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 01:22:54 +0200 Subject: [PATCH 23/42] Use generator for points --- cea/resources/radiation/geometry_generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 38e93268fa..58b292bcd8 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -540,7 +540,7 @@ def generate_tin(self): _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) From 2e7ee37ef5d4d8fa1cb8487a706b39bc8aaffcc7 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 01:28:22 +0200 Subject: [PATCH 24/42] Store nodata value for elevation map This allows us to ignore invalid values in the elevation map. This could also possibly allow rendering of terrain values below zero, since we can check for intersections below zero. --- cea/resources/radiation/geometry_generator.py | 38 +++++++++---------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 58b292bcd8..f8fd129fc7 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -74,17 +74,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: @@ -374,11 +375,9 @@ def burn_buildings(geometry, elevation_map): terrain_tin = elevation_map.generate_tin() # 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()) @@ -492,29 +491,27 @@ 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', 'nodata'] - def __init__(self, elevation_map, x_coords, y_coords): + def __init__(self, elevation_map, x_coords, y_coords, nodata=None): self.elevation_map = elevation_map self.x_coords = x_coords self.y_coords = y_coords + self.nodata = nodata @classmethod def read_raster(cls, raster, raise_above_sea_level=True): 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() 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, nodata) def get_elevation_map_from_geometry(self, geometry, extra_points=5): minx, miny, maxx, maxy = geometry.bounds @@ -536,7 +533,8 @@ def get_elevation_map_from_geometry(self, geometry, extra_points=5): return ElevationMap(new_elevation_map, new_x_coords, new_y_coords) def generate_tin(self): - (y_index, x_index) = np.nonzero(self.elevation_map >= 0) + # 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] From a69bbf51390010f31a3a7669acc0859cfccd8926 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 01:29:54 +0200 Subject: [PATCH 25/42] Add tolerance parameter --- cea/resources/radiation/geometry_generator.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index f8fd129fc7..f0d1e29bd0 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -141,7 +141,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) @@ -359,7 +359,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) @@ -372,7 +372,7 @@ 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) @@ -532,7 +532,7 @@ def get_elevation_map_from_geometry(self, geometry, extra_points=5): return ElevationMap(new_elevation_map, new_x_coords, new_y_coords) - def generate_tin(self): + 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] @@ -540,7 +540,7 @@ def generate_tin(self): 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 From 185fa043e89bc6f892d4d0f2f7cd4ac82465ba4b Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 01:39:37 +0200 Subject: [PATCH 26/42] Remove unused parameter --- cea/resources/radiation/geometry_generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index f0d1e29bd0..c980a0013e 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -500,7 +500,7 @@ def __init__(self, elevation_map, x_coords, y_coords, nodata=None): 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() From 8cd2557f02f03f9e8f907a60f0060215719a8adc Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 04:01:38 +0200 Subject: [PATCH 27/42] Take into account the grid size when checking bounds Since the stored coordinates in the elevation map are the mid points of each grid, we need to consider the left/right/top/bottom side of the grid when comparing the bounds of buildings. --- cea/resources/radiation/geometry_generator.py | 28 ++++++++++++------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index c980a0013e..97fcab05c3 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -491,12 +491,16 @@ def calc_intersection_face_solid(potentially_intersecting_solid, point): class ElevationMap(object): - __slots__ = ['elevation_map', 'x_coords', 'y_coords', 'nodata'] + __slots__ = ['elevation_map', 'x_coords', 'y_coords', 'x_size', 'y_size', 'nodata'] - def __init__(self, elevation_map, x_coords, y_coords, nodata=None): + 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 @@ -506,20 +510,24 @@ def read_raster(cls, raster): a = band.ReadAsArray() - (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, nodata) + return cls(a, x_coords, y_coords, x_size, y_size, nodata) def get_elevation_map_from_geometry(self, geometry, extra_points=5): 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] + x_start = np.where(minx > self.x_coords - self.x_size)[0] + x_end = np.where(maxx < self.x_coords + self.x_size)[0] + y_start = np.where(maxy < self.y_coords - self.y_size)[0] + y_end = np.where(miny > self.y_coords + self.y_size)[0] x_start = max(x_start[-1] - extra_points, 0) x_end = min(x_end[0] + extra_points, len(self.x_coords)) @@ -530,7 +538,7 @@ def get_elevation_map_from_geometry(self, geometry, extra_points=5): 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, tolerance=1e-6): # Ignore no data values from raster From fb9e1d1dbd82ab928e59a8aca889e192bfb87635 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 04:20:43 +0200 Subject: [PATCH 28/42] Use np.searchsorted to find relevant index Should provide better performance for higher resolution terrain --- cea/resources/radiation/geometry_generator.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 97fcab05c3..3ad3456b18 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -524,15 +524,16 @@ def read_raster(cls, raster): def get_elevation_map_from_geometry(self, geometry, extra_points=5): minx, miny, maxx, maxy = geometry.bounds - x_start = np.where(minx > self.x_coords - self.x_size)[0] - x_end = np.where(maxx < self.x_coords + self.x_size)[0] - y_start = np.where(maxy < self.y_coords - self.y_size)[0] - y_end = np.where(miny > self.y_coords + self.y_size)[0] - - 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)) + 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') + + # 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] From 65b4b3afe350c98540871588b17067f8d69be03d Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 04:21:20 +0200 Subject: [PATCH 29/42] Check geometry bounds with elevation map --- cea/resources/radiation/geometry_generator.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 3ad3456b18..43675ba33c 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -524,6 +524,11 @@ def read_raster(cls, raster): def get_elevation_map_from_geometry(self, geometry, extra_points=5): minx, miny, maxx, maxy = geometry.bounds + # 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 From eb7ab673907cc7f6cf240b624cfe88463f7b1574 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 04:22:31 +0200 Subject: [PATCH 30/42] Reduce number of extra points considered --- cea/resources/radiation/geometry_generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 43675ba33c..669fec39a3 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -521,7 +521,7 @@ def read_raster(cls, raster): 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 # Ensure geometry bounds is within elevation map From 584899ce40d2809c84216733564c6b25a7641424 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 14:25:14 +0200 Subject: [PATCH 31/42] Make sure that the scene is in projected crs --- cea/resources/radiation/geometry_generator.py | 24 +++++++++++++------ 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 669fec39a3..0cd2fb00eb 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -21,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 @@ -38,6 +38,8 @@ __email__ = "cea@arch.ethz.ch" __status__ = "Production" +from cea.utilities.standardize_coordinates import get_lat_lon_projected_shapefile, get_projected_coordinate_system + def identify_surfaces_type(occface_list): roof_list = [] @@ -560,13 +562,21 @@ def generate_tin(self, tolerance=1e-6): 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) + + 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 + ) - return zone_df, surroundings_df, terrain_raster + return reprojected_zone_df, reprojected_surroundings_df, reprojected_terrain def check_terrain_bounds(zone_df, surroundings_df, terrain_raster): From 605f32885a265b4f2dc8bc00d6361e5fb807d704 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 15:49:49 +0200 Subject: [PATCH 32/42] Add more information regarding simplification level --- cea/default.config | 8 ++++---- cea/resources/radiation/geometry_generator.py | 3 +-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/cea/default.config b/cea/default.config index 17a1f00927..95da8f926a 100644 --- a/cea/default.config +++ b/cea/default.config @@ -159,13 +159,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 diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 0cd2fb00eb..5338180d42 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -118,8 +118,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) From c45b047b900936787c27c10eee34fc2388b6b394 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Sat, 8 Jun 2024 15:58:05 +0200 Subject: [PATCH 33/42] Show reprojected crs --- cea/resources/radiation/geometry_generator.py | 4 +++- cea/utilities/standardize_coordinates.py | 4 ++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/cea/resources/radiation/geometry_generator.py b/cea/resources/radiation/geometry_generator.py index 5338180d42..2df53a6f91 100644 --- a/cea/resources/radiation/geometry_generator.py +++ b/cea/resources/radiation/geometry_generator.py @@ -38,7 +38,8 @@ __email__ = "cea@arch.ethz.ch" __status__ = "Production" -from cea.utilities.standardize_coordinates import get_lat_lon_projected_shapefile, get_projected_coordinate_system +from cea.utilities.standardize_coordinates import (get_lat_lon_projected_shapefile, get_projected_coordinate_system, + crs_to_epsg) def identify_surfaces_type(occface_list): @@ -575,6 +576,7 @@ def standardize_coordinate_systems(zone_df, surroundings_df, terrain_raster): dstSRS=crs ) + print(f"Reprojected scene to `EPSG:{crs_to_epsg(crs)}`") return reprojected_zone_df, reprojected_surroundings_df, reprojected_terrain diff --git a/cea/utilities/standardize_coordinates.py b/cea/utilities/standardize_coordinates.py index 6cb3089cf5..eb2d42dbc0 100644 --- a/cea/utilities/standardize_coordinates.py +++ b/cea/utilities/standardize_coordinates.py @@ -60,6 +60,10 @@ def get_projected_coordinate_system(lat, lon): 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): data = data.to_crs(get_geographic_coordinate_system()) lon = data.geometry[0].centroid.coords.xy[0][0] From e4c46c1c0bf82269a06ad206c03b46f930244ea1 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 10 Jun 2024 14:59:11 +0200 Subject: [PATCH 34/42] Rewrite terrain helper to fetch from terrain tiles --- cea/datamanagement/terrain_helper.py | 276 +++++++++++++++------------ 1 file changed, 149 insertions(+), 127 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index 9dfc93c197..0b7217d430 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -1,145 +1,167 @@ -""" -This script extracts terrain elevation from NASA - SRTM -https://www2.jpl.nasa.gov/srtm/ -""" - - - - +import math import os +from typing import Iterable, Tuple +import geopandas as gpd import numpy as np import pandas as pd -import requests -from geopandas import GeoDataFrame as Gdf -from osgeo import gdal, ogr, osr -from shapely.geometry import Polygon +import rasterio +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) - """ - # 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 +URL_FORMAT = "https://s3.amazonaws.com/elevation-tiles-prod/geotiff/{zoom}/{x}/{y}.tif" -def main(config): +def get_tile_number(lat: float, lon: float, zoom: int) -> Tuple[int, int]: """ - Create the terrain.tif file + Converts latitude and longitude to slippy map tile coordinates. + Reference: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames + """ + lat_rad = math.radians(lat) + + 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 - :param config: - :type config: cea.config.Configuration - :return: + +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. """ - assert os.path.exists(config.scenario), 'Scenario not found: %s' % config.scenario + # 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]): + rasters = [] + try: + for tile_url in tile_urls: + rasters.append(rasterio.open(tile_url)) + + # 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 + + +def reproject_raster_array(src_array, src_transform, meta, dst_crs): + # 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 + ) + + 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, min_y, max_x, max_y, zoom=12): + tile_numbers = get_all_tile_numbers(min_x, min_y, max_x, max_y, zoom) + + # Get merged raster array + dest, transform, meta = merge_raster_tiles(URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers) + + # Reproject raster array to bounds crs + dest, transform, meta = reproject_raster_array(dest, transform, meta, 'EPSG:4326') + + with MemoryFile() as memfile: + with memfile.open(**meta) as dataset: + dataset.write(dest) + + bounding_box = box(min_x, min_y, max_x, max_y) + + out_dest, out_transform = mask(dataset, [bounding_box], 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 + + +def main(config): 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()) + total_df = pd.concat([zone_df, surroundings_df]) + total_bounds = total_df.to_crs("EPSG:4326").total_bounds + + dest, transform, meta = fetch_tiff(*total_bounds) + + # Write to disk + os.makedirs(os.path.dirname(locator.get_terrain()), exist_ok=True) + with rasterio.open(locator.get_terrain(), "w", **meta) as f: + f.write(dest) if __name__ == '__main__': From fb3dc4c08511f1b60b86a84192b5fa6b64022e32 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 10 Jun 2024 15:14:13 +0200 Subject: [PATCH 35/42] Update docstring and typing --- cea/datamanagement/terrain_helper.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index 0b7217d430..8f7856d9ac 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -1,11 +1,12 @@ import math import os -from typing import Iterable, Tuple +from typing import Iterable, Tuple, Dict, Union import geopandas as gpd import numpy as np import pandas as pd import rasterio +from pyproj import CRS from rasterio import MemoryFile from rasterio.mask import mask from rasterio.merge import merge @@ -66,7 +67,10 @@ def get_all_tile_numbers(min_x: float, min_y: float, max_x: float, max_y: float, return all_tile_numbers -def merge_raster_tiles(tile_urls: Iterable[str]): +def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + """ + Merge raster files from tile urls into a single raster. + """ rasters = [] try: for tile_url in tile_urls: @@ -88,7 +92,10 @@ def merge_raster_tiles(tile_urls: Iterable[str]): return dest, transform, meta -def reproject_raster_array(src_array, src_transform, meta, dst_crs): +def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, dst_crs: Union[CRS, dict]) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + """ + Reproject raster array to specified CRS. + """ # Get bounds from transform minx, miny = src_transform * (0, src_array.shape[2]) maxx, maxy = src_transform * (src_array.shape[1], 0) @@ -121,14 +128,14 @@ def reproject_raster_array(src_array, src_transform, meta, dst_crs): return dst_array, transform, new_meta -def fetch_tiff(min_x, min_y, max_x, max_y, zoom=12): +def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12) -> Tuple[np.ndarray, rasterio.Affine, Dict]: tile_numbers = get_all_tile_numbers(min_x, min_y, max_x, max_y, zoom) # Get merged raster array dest, transform, meta = merge_raster_tiles(URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers) # Reproject raster array to bounds crs - dest, transform, meta = reproject_raster_array(dest, transform, meta, 'EPSG:4326') + dest, transform, meta = reproject_raster_array(dest, transform, meta, CRS.from_string('EPSG:4326')) with MemoryFile() as memfile: with memfile.open(**meta) as dataset: From 62f0ed33910193b5bc6e110feb92ca7f51d83870 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 10 Jun 2024 16:18:02 +0200 Subject: [PATCH 36/42] Standardise CRS --- cea/datamanagement/terrain_helper.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index 8f7856d9ac..b17acd28b5 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -17,6 +17,8 @@ import cea.inputlocator 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) def get_tile_number(lat: float, lon: float, zoom: int) -> Tuple[int, int]: @@ -32,7 +34,8 @@ def get_tile_number(lat: float, lon: float, zoom: int) -> Tuple[int, int]: return x_tile, y_tile -def get_all_tile_numbers(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int) -> Iterable[Tuple[int, int]]: +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. """ @@ -92,7 +95,8 @@ def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.A return dest, transform, meta -def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, dst_crs: Union[CRS, dict]) -> Tuple[np.ndarray, rasterio.Affine, Dict]: +def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, + dst_crs: Union[CRS, dict]) -> Tuple[np.ndarray, rasterio.Affine, Dict]: """ Reproject raster array to specified CRS. """ @@ -129,14 +133,16 @@ def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, dst def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + tile_numbers = get_all_tile_numbers(min_x, min_y, max_x, max_y, zoom) # Get merged raster array dest, transform, meta = merge_raster_tiles(URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers) # Reproject raster array to bounds crs - dest, transform, meta = reproject_raster_array(dest, transform, meta, CRS.from_string('EPSG:4326')) + dest, transform, meta = reproject_raster_array(dest, transform, meta, DEFAULT_CRS) + # Crop raster based on bounds with MemoryFile() as memfile: with memfile.open(**meta) as dataset: dataset.write(dest) @@ -159,9 +165,9 @@ def main(config): # Get total bounds zone_df = gpd.read_file(locator.get_zone_geometry()) - surroundings_df = gpd.read_file(locator.get_surroundings_geometry()) + surroundings_df = gpd.read_file(locator.get_surroundings_geometry()).to_crs(zone_df.crs) total_df = pd.concat([zone_df, surroundings_df]) - total_bounds = total_df.to_crs("EPSG:4326").total_bounds + total_bounds = total_df.to_crs(DEFAULT_CRS).total_bounds dest, transform, meta = fetch_tiff(*total_bounds) From 43111245498b079f71482ba8b983c4ba0a7911cd Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Mon, 10 Jun 2024 17:37:37 +0200 Subject: [PATCH 37/42] Add buffer to bounds to ensure terrain overlap --- cea/datamanagement/terrain_helper.py | 35 ++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index b17acd28b5..e99296978f 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -132,24 +132,31 @@ def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, return dst_array, transform, new_meta -def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12) -> Tuple[np.ndarray, rasterio.Affine, Dict]: +def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12, + src_crs: Union[CRS, dict] = DEFAULT_CRS) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + """ + 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. + """ + + bounding_box = gpd.GeoDataFrame(geometry=[box(min_x, min_y, max_x, max_y)], crs=src_crs) - tile_numbers = get_all_tile_numbers(min_x, min_y, max_x, max_y, zoom) + # 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 dest, transform, meta = merge_raster_tiles(URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers) # Reproject raster array to bounds crs - dest, transform, meta = reproject_raster_array(dest, transform, meta, DEFAULT_CRS) + dest, transform, meta = reproject_raster_array(dest, transform, meta, src_crs) # Crop raster based on bounds with MemoryFile() as memfile: with memfile.open(**meta) as dataset: dataset.write(dest) - bounding_box = box(min_x, min_y, max_x, max_y) - - out_dest, out_transform = mask(dataset, [bounding_box], crop=True) + out_dest, out_transform = mask(dataset, bounding_box.geometry, crop=True) out_meta = dataset.meta.copy() out_meta.update({ "height": out_dest.shape[1], @@ -166,10 +173,18 @@ def main(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 = pd.concat([zone_df, surroundings_df]) - total_bounds = total_df.to_crs(DEFAULT_CRS).total_bounds - - dest, transform, meta = fetch_tiff(*total_bounds) + 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 = fetch_tiff(min_x, min_y, max_x, max_y, src_crs=buffer_df.crs) # Write to disk os.makedirs(os.path.dirname(locator.get_terrain()), exist_ok=True) From a1284f5bed1c354f6eaffddc542d76597c05115b Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Tue, 11 Jun 2024 11:55:42 +0200 Subject: [PATCH 38/42] Add attribution to data --- cea/datamanagement/terrain_helper.py | 52 ++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 3 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index e99296978f..1a18a357c7 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -1,3 +1,4 @@ +import datetime import math import os from typing import Iterable, Tuple, Dict, Union @@ -6,6 +7,7 @@ import numpy as np import pandas as pd import rasterio +import requests from pyproj import CRS from rasterio import MemoryFile from rasterio.mask import mask @@ -20,6 +22,30 @@ 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]: """ @@ -77,7 +103,8 @@ def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.A rasters = [] try: for tile_url in tile_urls: - rasters.append(rasterio.open(tile_url)) + raster = rasterio.open(tile_url) + rasters.append(raster) # Merge rasters dest, transform = merge(rasters) @@ -140,13 +167,15 @@ def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int """ 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 - dest, transform, meta = merge_raster_tiles(URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers) + tile_urls = [URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers] + dest, transform, meta = merge_raster_tiles(tile_urls) # Reproject raster array to bounds crs dest, transform, meta = reproject_raster_array(dest, transform, meta, src_crs) @@ -186,8 +215,25 @@ def main(config): # Fetch tiff data dest, transform, meta = fetch_tiff(min_x, min_y, max_x, max_y, src_crs=buffer_df.crs) - # Write to disk 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") + content = (f"\nCitation:\n" + f"Terrain Tiles was accessed on {datetime.datetime.now().date()} " + f"from https://registry.opendata.aws/terrain-tiles.\n" + f"\nAcknowledgement 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" + f"") + + 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) From 3027a3c2478a8de39d4f894b93574ca9c7e4084f Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Tue, 11 Jun 2024 11:57:28 +0200 Subject: [PATCH 39/42] Remove unused import --- cea/datamanagement/terrain_helper.py | 1 - 1 file changed, 1 deletion(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index 1a18a357c7..62b8ef142a 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -7,7 +7,6 @@ import numpy as np import pandas as pd import rasterio -import requests from pyproj import CRS from rasterio import MemoryFile from rasterio.mask import mask From ec7c390e7492caa64e0509a19b1fce129ba4ae45 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Tue, 11 Jun 2024 18:10:12 +0200 Subject: [PATCH 40/42] Add tile information to references --- cea/datamanagement/terrain_helper.py | 46 ++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 12 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index 62b8ef142a..ea6d184f96 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -1,12 +1,14 @@ import datetime +import io import math import os -from typing import Iterable, Tuple, Dict, Union +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 pyproj import CRS from rasterio import MemoryFile from rasterio.mask import mask @@ -95,16 +97,28 @@ def get_all_tile_numbers(min_x: float, min_y: float, max_x: float, max_y: float, return all_tile_numbers -def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.Affine, Dict]: +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: - raster = rasterio.open(tile_url) + 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() @@ -118,7 +132,7 @@ def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.A for raster in rasters: raster.close() - return dest, transform, meta + return dest, transform, meta, tile_info def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, @@ -159,10 +173,11 @@ def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12, - src_crs: Union[CRS, dict] = DEFAULT_CRS) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + src_crs: Union[CRS, dict] = DEFAULT_CRS) -> Tuple[np.ndarray, rasterio.Affine, Dict, List[Dict]]: """ 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) @@ -174,7 +189,7 @@ def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int # Get merged raster array tile_urls = [URL_FORMAT.format(zoom=zoom, x=x, y=y) for x, y in tile_numbers] - dest, transform, meta = merge_raster_tiles(tile_urls) + 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) @@ -192,7 +207,7 @@ def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int "transform": out_transform }) - return out_dest, out_transform, out_meta + return out_dest, out_transform, out_meta, tile_info def main(config): @@ -212,20 +227,27 @@ def main(config): min_x, min_y, max_x, max_y = buffer_df.total_bounds # Fetch tiff data - dest, transform, meta = fetch_tiff(min_x, min_y, max_x, max_y, src_crs=buffer_df.crs) + dest, transform, meta, tile_info = fetch_tiff(min_x, min_y, max_x, max_y, 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") - content = (f"\nCitation:\n" + 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"\nAcknowledgement of Data Sources:" + 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" - f"") + f"Attribution: {ATTRIBUTION_URL}\n") with open(reference_file, "w") as f: f.write(content) From 320d6bd466d0da30aa9389274bd59f13278983e0 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Tue, 11 Jun 2024 18:21:54 +0200 Subject: [PATCH 41/42] Add grid size as parametter --- cea/datamanagement/terrain_helper.py | 13 ++++++++----- cea/default.config | 5 ----- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/cea/datamanagement/terrain_helper.py b/cea/datamanagement/terrain_helper.py index ea6d184f96..ee72855162 100644 --- a/cea/datamanagement/terrain_helper.py +++ b/cea/datamanagement/terrain_helper.py @@ -136,7 +136,7 @@ def merge_raster_tiles(tile_urls: Iterable[str]) -> Tuple[np.ndarray, rasterio.A def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, - dst_crs: Union[CRS, dict]) -> Tuple[np.ndarray, rasterio.Affine, Dict]: + dst_crs: Union[CRS, dict], grid_size: int) -> Tuple[np.ndarray, rasterio.Affine, Dict]: """ Reproject raster array to specified CRS. """ @@ -146,7 +146,8 @@ def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, transform, width, height = calculate_default_transform( meta["crs"], dst_crs, src_array.shape[2], src_array.shape[1], - minx, miny, maxx, maxy + minx, miny, maxx, maxy, + resolution=(grid_size, grid_size) ) new_meta = meta.copy() @@ -172,7 +173,7 @@ def reproject_raster_array(src_array: np.ndarray, src_transform, meta: Dict, return dst_array, transform, new_meta -def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int = 12, +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]]: """ Fetch raster data array based on bounds in the given source CRS (Default: lat, lon). @@ -192,7 +193,7 @@ def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int 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) + dest, transform, meta = reproject_raster_array(dest, transform, meta, src_crs, grid_size) # Crop raster based on bounds with MemoryFile() as memfile: @@ -211,6 +212,7 @@ def fetch_tiff(min_x: float, min_y: float, max_x: float, max_y: float, zoom: int def main(config): + grid_size = config.terrain_helper.grid_size locator = cea.inputlocator.InputLocator(config.scenario) # Get total bounds @@ -227,7 +229,8 @@ def main(config): 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, src_crs=buffer_df.crs) + 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 diff --git a/cea/default.config b/cea/default.config index 95da8f926a..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 From 8369a1410349e40c1d07b94583b9d25abf026c19 Mon Sep 17 00:00:00 2001 From: Reynold Mok <34395415+reyery@users.noreply.github.com> Date: Tue, 11 Jun 2024 18:26:43 +0200 Subject: [PATCH 42/42] Update description --- cea/scripts.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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]