diff --git a/.ipynb_checkpoints/main-checkpoint.py b/.ipynb_checkpoints/main-checkpoint.py new file mode 100755 index 0000000..d3e6be0 --- /dev/null +++ b/.ipynb_checkpoints/main-checkpoint.py @@ -0,0 +1,636 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon May 7 14:23:00 2018 + + Copyright (C) <2020> + Michael.Neuhauser@bfw.gv.at + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +""" +# import standard libraries +import os +import sys +import psutil +import numpy as np +from datetime import datetime +from multiprocessing import cpu_count +import multiprocessing as mp +import logging +from xml.etree import ElementTree as ET + +# Flow-Py Libraries +import raster_io as io +import Simulation as Sim +import flow_core as fc + +# Libraries for GUI, PyQt5 +from PyQt5.QtCore import pyqtSlot, QCoreApplication +from PyQt5.QtWidgets import QFileDialog, QMessageBox, QMainWindow, QApplication + +from Flow_GUI import Ui_MainWindow + + +class Flow_Py_EXEC(): + + def __init__(self): + + app = QApplication(sys.argv) + MainWindow = QMainWindow() + self.ui = Ui_MainWindow() + self.ui.setupUi(MainWindow) + + # self.showMaximized() + #self.ui.setWindowTitle("Flow-Py") + #self.ui.setWindowIcon(QIcon('logo.jpg')) + + self.directory = os.getcwd() + + self.ui.alpha_Edit.setText('25') + self.ui.exp_Edit.setText('8') + self.ui.flux_Edit.setText('0.003') + self.ui.z_Edit.setText('8848') + + self.ui.wDir_Button.clicked.connect(self.open_wDir) + self.ui.DEM_Button.clicked.connect(self.open_dhm) + self.ui.Release_Button.clicked.connect(self.open_release) + self.ui.infra_Button.clicked.connect(self.open_infra) + self.ui.forest_Button.clicked.connect(self.open_forest) + #self.ui.process_Box.currentIndexChanged.connect(self.processChanged) + self.ui.calc_Button.clicked.connect(self.calculation) + self.ui.actionSave.triggered.connect(self.save) + self.ui.actionLoad.triggered.connect(self.load) + self.ui.actionQuit.triggered.connect(self.quit) + + self.calc_class = None + self.prot_for_bool = False + self.threads_calc = 0 + self.progress_value = 0 + self.cpu_count = 1 + self.thread_list = [] + self.start_list = [] + self.end_list = [] + for i in range(self.cpu_count): + self.thread_list.append(0) + self.start_list.append(0) + self.end_list.append(0) + + # show the constructed window + MainWindow.show() + sys.exit(app.exec_()) + + def set_gui_bool(self, bool): + self.ui.calc_Button.setEnabled(bool) + self.ui.wDir_lineEdit.setEnabled(bool) + self.ui.DEM_lineEdit.setEnabled(bool) + self.ui.release_lineEdit.setEnabled(bool) + self.ui.infra_lineEdit.setEnabled(bool) + self.ui.forest_lineEdit.setEnabled(bool) + + def save(self): + """Save the input paths""" + name = QFileDialog.getSaveFileName(None, 'Save File', + ".xml")[0] + if len(name) != 0: + + root = ET.Element('root') + wdir = ET.SubElement(root, 'wDir') + dhm = ET.SubElement(root, 'DHM') + release = ET.SubElement(root, 'Release') + infra = ET.SubElement(root, 'Infrastructure') + forest = ET.SubElement(root, 'Forest') + + wdir.set('Directory', 'Working') + dhm.set('Directory', 'DHM') + release.set('Directory', 'Release') + infra.set('Directory', 'Infrastructure') + forest.set('Directory', 'Forest') + + wdir.text = self.ui.wDir_lineEdit.text() + dhm.text = self.ui.DEM_lineEdit.text() + release.text = self.ui.release_lineEdit.text() + infra.text = self.ui.infra_lineEdit.text() + forest.text = self.ui.forest_lineEdit.text() + + tree = ET.ElementTree(root) + tree.write(name) + + def load(self): + xml_file = QFileDialog.getOpenFileNames(None, 'Open xml', + self.directory, + "xml (*.xml);;All Files (*.*)")[0] + + if len(xml_file) != 0: + tree = ET.parse(xml_file[0]) + root = tree.getroot() + + try: + self.ui.wDir_lineEdit.setText(root[0].text) + self.directory = root[0].text + except: + print("No Working Directory Path in File!") + + try: + self.ui.DEM_lineEdit.setText(root[1].text) + except: + print("No DEM Path in File!") + + try: + self.ui.release_lineEdit.setText(root[2].text) + except: + print("No Release Path in File!") + + try: + self.ui.infra_lineEdit.setText(root[3].text) + except: + print("No Infrastructure Path in File!") + + try: + self.ui.forest_lineEdit.setText(root[4].text) + except: + print("No Forest Path in File!") + + def quit(self): + QCoreApplication.quit() + + def open_wDir(self): + """Open the Working Directory, where results are stored""" + directory = QFileDialog.getExistingDirectory(None, 'Open Working Directory', + self.directory, + QFileDialog.ShowDirsOnly) + if len(directory) != 0: + self.directory = directory + self.ui.wDir_lineEdit.setText(self.directory) + + def open_dhm(self): + """Open digital elevation model""" + dem_file = QFileDialog.getOpenFileNames(None, 'Open DEM', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(dem_file[0]) != 0: + dem = dem_file[0] + self.ui.DEM_lineEdit.setText(dem[0]) + + def open_release(self): + """Open release layer""" + release_file = QFileDialog.getOpenFileNames(None, 'Open Release', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(release_file[0]) != 0: + release = release_file[0] + self.ui.release_lineEdit.setText(release[0]) + + def open_infra(self): + """Open infrastructure layer""" + infra_file = QFileDialog.getOpenFileNames(None, 'Open Infrastructure Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(infra_file[0]) != 0: + infra = infra_file[0] + self.ui.infra_lineEdit.setText(infra[0]) + + def open_forest(self): + """Open forest layer""" + forest_file = QFileDialog.getOpenFileNames(None, 'Open Forest Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + forest = forest_file[0] + self.ui.forest_lineEdit.setText(forest[0]) + + def update_progressBar(self, float, thread, start, end): + self.thread_list[thread] = float + self.start_list[thread] = start + self.end_list[thread] = end + + self.progress_value = sum(self.thread_list) / len(self.thread_list) + self.progressBar.setValue(self.progress_value) + for i in range(len(self.thread_list)): + sys.stdout.write( + "Thread {}: Startcell {} of {} = {}%"'\n'.format(i + 1, self.start_list[i], self.end_list[i], + self.thread_list[i])) + sys.stdout.flush() + for i in range(len(self.thread_list)): + sys.stdout.write('\x1b[1A' + '\x1b[2K') # Go 1 line up and erase that line + + @staticmethod + def showdialog(path): + msg = QMessageBox() + msg.setIcon(QMessageBox.Critical) + msg.setText("No " + path + " set") + msg.setWindowTitle("Error") + msg.setStandardButtons(QMessageBox.Ok) + msg.exec_() + + def calculation(self): + self.start = datetime.now().replace(microsecond=0) + self.calc_bool = False + + # Check if input is ok + if self.ui.wDir_lineEdit.text() == '': + self.showdialog('Working Directory') + return + if self.ui.DEM_lineEdit.text() == '': + self.showdialog('DEM Layer') + return + if self.ui.release_lineEdit.text() == '': + self.showdialog('Release Layer') + return + # Disable all input line Edits and Buttons + self.set_gui_bool(False) + + # Create result directory + time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + try: + os.makedirs(self.ui.wDir_lineEdit.text() + '/res_{}/'.format(time_string)) + self.res_dir = ('/res_{}/'.format(time_string)) + except FileExistsError: + self.res_dir = ('/res_{}/'.format(time_string)) + + # Setup logger + + for handler in logging.root.handlers[:]: + logging.root.removeHandler(handler) + + logging.basicConfig(level=logging.INFO, + format='%(asctime)s %(levelname)-8s %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + filename=(self.directory + self.res_dir + 'log_{}.txt').format(time_string), + filemode='w') + + # Start of Calculation + logging.info('Start Calculation') + # Read in raster files + try: + dem, header = io.read_raster(self.ui.DEM_lineEdit.text()) + logging.info('DEM File: {}'.format(self.ui.DEM_lineEdit.text())) + except FileNotFoundError: + print("Wrong filepath or filename") + self.set_gui_bool(True) + return + + try: + release, release_header = io.read_raster(self.ui.release_lineEdit.text()) + logging.info('Release File: {}'.format(self.ui.release_lineEdit.text())) + except FileNotFoundError: + print("Wrong filepath or filename") + self.set_gui_bool(True) + return + + # Check if Layers have same size!!! + if header['ncols'] == release_header['ncols'] and header['nrows'] == release_header['nrows']: + print("DEM and Release Layer ok!") + else: + print("Error: Release Layer doesn't match DEM!") + self.set_gui_bool(True) + return + + try: + infra, infra_header = io.read_raster(self.ui.infra_lineEdit.text()) + if header['ncols'] == infra_header['ncols'] and header['nrows'] == infra_header['nrows']: + print("Infra Layer ok!") + self.calc_bool = True + logging.info('Infrastructure File: {}'.format(self.ui.infra_lineEdit.text())) + else: + print("Error: Infra Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(self.ui.forest_lineEdit.text()) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + self.prot_for_bool = True + logging.info('Forest File: {}'.format(self.forest_lineEdit.text())) + else: + print("Error: Forest Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + forest = np.zeros_like(dem) + + logging.info('Files read in') + + alpha = self.ui.alpha_Edit.text() + exp = self.ui.exp_Edit.text() + flux_threshold = self.ui.flux_Edit.text() + max_z = self.ui.z_Edit.text() + + logging.info('Alpha Angle: {}'.format(alpha)) + logging.info('Exponent: {}'.format(exp)) + logging.info('Flux Threshold: {}'.format(flux_threshold)) + logging.info('Max Z_delta: {}'.format(max_z)) + logging.info + + self.z_delta = np.zeros_like(dem) + self.flux = np.zeros_like(dem) + self.cell_counts = np.zeros_like(dem) + self.z_delta_sum = np.zeros_like(dem) + self.backcalc = np.zeros_like(dem) + self.fp_ta = np.zeros_like(dem) + self.sl_ta = np.zeros_like(dem) + + # Calculation + self.calc_class = Sim.Simulation(dem, header, release, release_header, infra, forest, self.calc_bool, alpha, exp, flux_threshold, max_z) + self.calc_class.value_changed.connect(self.update_progressBar) + self.calc_class.finished.connect(self.thread_finished) + logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) + self.calc_class.start() + + def thread_finished(self, z_delta, flux, count_array, z_delta_sum, backcalc, fp_ta, sl_ta): + logging.info('Calculation finished, getting results.') + for i in range(len(z_delta)): + self.z_delta = np.maximum(self.z_delta, z_delta[i]) + self.flux = np.maximum(self.flux, flux[i]) + self.cell_counts += count_array[i] + self.z_delta_sum += z_delta_sum[i] + self.backcalc = np.maximum(self.backcalc, backcalc[i]) + self.fp_ta = np.maximum(self.fp_ta, fp_ta[i]) + self.sl_ta = np.maximum(self.sl_ta, sl_ta[i]) + self.output() + + def output(self): + # time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + logging.info('Writing Output Files') + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "flux{}".format(self.ui.outputBox.currentText()), + self.flux) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "z_delta{}".format(self.ui.outputBox.currentText()), + self.z_delta) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "FP_travel_angle{}".format(self.ui.outputBox.currentText()), + self.fp_ta) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "SL_travel_angle{}".format(self.ui.outputBox.currentText()), + self.sl_ta) + if not self.calc_bool: + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "cell_counts{}".format(self.ui.outputBox.currentText()), + self.cell_counts) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "z_delta_sum{}".format(self.ui.outputBox.currentText()), + self.z_delta_sum) + if self.calc_bool: + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "backcalculation{}".format(self.ui.outputBox.currentText()), + self.backcalc) + + print("Calculation finished") + end = datetime.now().replace(microsecond=0) + logging.info('Calculation needed: ' + str(end - self.start) + ' seconds') + + # Handle GUI + #self.ui.progressBar.setValue(100) + self.set_gui_bool(True) + + +def main(args, kwargs): + + + alpha = args[0] + exp = args[1] + directory = args[2] + dem_path = args[3] + release_path = args[4] + if 'infra' in kwargs: + infra_path = kwargs.get('infra') + #print(infra_path) + else: + infra_path = None + + if 'forest' in kwargs: + forest_path = kwargs.get('forest') + #print(forest_path) + else: + forest_path = None + + if 'flux' in kwargs: + flux_threshold = float(kwargs.get('flux')) + #print(flux_threshold) + else: + flux_threshold = 3 * 10 ** -4 + + if 'max_z' in kwargs: + max_z = kwargs.get('max_z') + #print(max_z) + else: + max_z = 8848 + # Recomendet values: + # Avalanche = 270 + # Rockfall = 50 + # Soil Slide = 12 + + print("Starting...") + print("...") + + start = datetime.now().replace(microsecond=0) + calc_bool = False + # Create result directory + time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + try: + os.makedirs(directory + '/res_{}/'.format(time_string)) + res_dir = ('/res_{}/'.format(time_string)) + except FileExistsError: + res_dir = ('/res_{}/'.format(time_string)) + + # Setup logger + for handler in logging.root.handlers[:]: + logging.root.removeHandler(handler) + + logging.basicConfig(level=logging.INFO, + format='%(asctime)s %(levelname)-8s %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + filename=(directory + res_dir + 'log_{}.txt').format(time_string), + filemode='w') + + # Start of Calculation + logging.info('Start Calculation') + logging.info('Alpha Angle: {}'.format(alpha)) + logging.info('Exponent: {}'.format(exp)) + logging.info('Flux Threshold: {}'.format(flux_threshold)) + logging.info('Max Z_delta: {}'.format(max_z)) + logging.info + # Read in raster files + try: + dem, header = io.read_raster(dem_path) + logging.info('DEM File: {}'.format(dem_path)) + except FileNotFoundError: + print("DEM: Wrong filepath or filename") + return + + try: + release, release_header = io.read_raster(release_path) + logging.info('Release File: {}'.format(release_path)) + except FileNotFoundError: + print("Wrong filepath or filename") + return + + # Check if Layers have same size!!! + if header['ncols'] == release_header['ncols'] and header['nrows'] == release_header['nrows']: + print("DEM and Release Layer ok!") + else: + print("Error: Release Layer doesn't match DEM!") + return + + try: + infra, infra_header = io.read_raster(infra_path) + if header['ncols'] == infra_header['ncols'] and header['nrows'] == infra_header['nrows']: + print("Infra Layer ok!") + calc_bool = True + logging.info('Infrastructure File: {}'.format(infra_path)) + else: + print("Error: Infra Layer doesn't match DEM!") + return + except: + infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(forest_path) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + prot_for_bool = True + logging.info('Forest File: {}'.format(forest_path)) + else: + print("Error: Forest Layer doesn't match DEM!") + return + except: + forest = np.zeros_like(dem) + + logging.info('Files read in') + + z_delta = np.zeros_like(dem) + flux = np.zeros_like(dem) + cell_counts = np.zeros_like(dem) + z_delta_sum = np.zeros_like(dem) + backcalc = np.zeros_like(dem) + fp_ta = np.zeros_like(dem) + sl_ta = np.zeros_like(dem) + + avaiable_memory = psutil.virtual_memory()[1] + needed_memory = sys.getsizeof(dem) + + max_number_procces = int(avaiable_memory / (needed_memory * 10)) + + print( + "There are {} Bytes of Memory avaiable and {} Bytes needed per process. \nMax. Nr. of Processes = {}".format( + avaiable_memory, needed_memory*10, max_number_procces)) + + # Calculation + logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) + + if calc_bool: + release_list = fc.split_release(release, release_header, min(mp.cpu_count() * 2, max_number_procces)) + + print("{} Processes started.".format(len(release_list))) + pool = mp.Pool(len(release_list)) + results = pool.map(fc.calculation, + [[dem, header, infra, forest, release_pixel, alpha, exp, flux_threshold, max_z] + for release_pixel in release_list]) + pool.close() + pool.join() + else: + release_list = fc.split_release(release, release_header, min(mp.cpu_count() * 4, max_number_procces)) + + print("{} Processes started.".format(len(release_list))) + pool = mp.Pool(mp.cpu_count()) + # results = pool.map(gc.calculation, iterable) + results = pool.map(fc.calculation_effect, + [[dem, header, forest, release_pixel, alpha, exp, flux_threshold, max_z] for + release_pixel in release_list]) + pool.close() + pool.join() + + z_delta_list = [] + flux_list = [] + cc_list = [] + z_delta_sum_list = [] + backcalc_list = [] + fp_ta_list = [] + sl_ta_list = [] + for i in range(len(results)): + res = results[i] + res = list(res) + z_delta_list.append(res[0]) + flux_list.append(res[1]) + cc_list.append(res[2]) + z_delta_sum_list.append(res[3]) + backcalc_list.append(res[4]) + fp_ta_list.append(res[5]) + sl_ta_list.append(res[6]) + + logging.info('Calculation finished, getting results.') + for i in range(len(z_delta_list)): + z_delta = np.maximum(z_delta, z_delta_list[i]) + flux = np.maximum(flux, flux_list[i]) + cell_counts += cc_list[i] + z_delta_sum += z_delta_sum_list[i] + backcalc = np.maximum(backcalc, backcalc_list[i]) + fp_ta = np.maximum(fp_ta, fp_ta_list[i]) + sl_ta = np.maximum(sl_ta, sl_ta_list[i]) + + + # time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + logging.info('Writing Output Files') + output_format = '.tif' + io.output_raster(dem_path, + directory + res_dir + "flux{}".format(output_format), + flux) + io.output_raster(dem_path, + directory + res_dir + "z_delta{}".format(output_format), + z_delta) + io.output_raster(dem_path, + directory + res_dir + "FP_travel_angle{}".format(output_format), + fp_ta) + io.output_raster(dem_path, + directory + res_dir + "SL_travel_angle{}".format(output_format), + sl_ta) + if not calc_bool: # if no infra + io.output_raster(dem_path, + directory + res_dir + "cell_counts{}".format(output_format), + cell_counts) + io.output_raster(dem_path, + directory + res_dir + "z_delta_sum{}".format( output_format), + z_delta_sum) + if calc_bool: # if infra + io.output_raster(dem_path, + directory + res_dir + "backcalculation{}".format(output_format), + backcalc) + + print("Calculation finished") + print("...") + end = datetime.now().replace(microsecond=0) + logging.info('Calculation needed: ' + str(end - start) + ' seconds') + + +if __name__ == '__main__': + #mp.set_start_method('spawn') # used in Windows + argv = sys.argv[1:] + #argv = ["28", "8", "../runs/open_slope/", "../runs/open_slope/parabola.asc", "../runs/open_slope/release.tif", "forest=../runs/open_slope/forest_deep.tif", "flux=0.0003" , "max_z=270"] #"forest=../runs/open_slope/forest.tif", + #argv = ['--gui'] + if len(argv) < 1: + print("Too few input arguments!!!") + sys.exit(1) + if len(argv) == 1 and argv[0] == '--gui': + Flow_Py_EXEC() + else: + args=[arg for arg in argv if arg.find('=')<0] + kwargs={kw[0]:kw[1] for kw in [ar.split('=') for ar in argv if ar.find('=')>0]} + + main(args, kwargs) + +# example dam: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif +# example dam/infra: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif infra=./examples/dam/infra.tif flux=0.003 max_z=270 +# example dam/forest: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif forest=./examples/dam/infra.tif flux=0.003 max_z=270 diff --git a/.ipynb_checkpoints/main_edit-checkpoint.py b/.ipynb_checkpoints/main_edit-checkpoint.py new file mode 100644 index 0000000..a1ce79f --- /dev/null +++ b/.ipynb_checkpoints/main_edit-checkpoint.py @@ -0,0 +1,636 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon May 7 14:23:00 2018 + + Copyright (C) <2020> + Michael.Neuhauser@bfw.gv.at + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +""" +# import standard libraries +import os +import sys +import psutil +import numpy as np +from datetime import datetime +from multiprocessing import cpu_count +import multiprocessing as mp +import logging +from xml.etree import ElementTree as ET + +# Flow-Py Libraries +import raster_io as io +import Simulation as Sim +import flow_core as fc + +# Libraries for GUI, PyQt5 +from PyQt5.QtCore import pyqtSlot, QCoreApplication +from PyQt5.QtWidgets import QFileDialog, QMessageBox, QMainWindow, QApplication + +from Flow_GUI import Ui_MainWindow + + +class Flow_Py_EXEC(): + + def __init__(self): + + app = QApplication(sys.argv) + MainWindow = QMainWindow() + self.ui = Ui_MainWindow() + self.ui.setupUi(MainWindow) + + # self.showMaximized() + #self.ui.setWindowTitle("Flow-Py") + #self.ui.setWindowIcon(QIcon('logo.jpg')) + + self.directory = os.getcwd() + + self.ui.alpha_Edit.setText('25') + self.ui.exp_Edit.setText('8') + self.ui.flux_Edit.setText('0.003') + self.ui.z_Edit.setText('8848') + + self.ui.wDir_Button.clicked.connect(self.open_wDir) + self.ui.DEM_Button.clicked.connect(self.open_dhm) + self.ui.Release_Button.clicked.connect(self.open_release) + self.ui.infra_Button.clicked.connect(self.open_infra) + self.ui.forest_Button.clicked.connect(self.open_forest) + #self.ui.process_Box.currentIndexChanged.connect(self.processChanged) + self.ui.calc_Button.clicked.connect(self.calculation) + self.ui.actionSave.triggered.connect(self.save) + self.ui.actionLoad.triggered.connect(self.load) + self.ui.actionQuit.triggered.connect(self.quit) + + self.calc_class = None + self.prot_for_bool = False + self.threads_calc = 0 + self.progress_value = 0 + self.cpu_count = 1 + self.thread_list = [] + self.start_list = [] + self.end_list = [] + for i in range(self.cpu_count): + self.thread_list.append(0) + self.start_list.append(0) + self.end_list.append(0) + + # show the constructed window + MainWindow.show() + sys.exit(app.exec_()) + + def set_gui_bool(self, bool): + self.ui.calc_Button.setEnabled(bool) + self.ui.wDir_lineEdit.setEnabled(bool) + self.ui.DEM_lineEdit.setEnabled(bool) + self.ui.release_lineEdit.setEnabled(bool) + self.ui.infra_lineEdit.setEnabled(bool) + self.ui.forest_lineEdit.setEnabled(bool) + + def save(self): + """Save the input paths""" + name = QFileDialog.getSaveFileName(None, 'Save File', + ".xml")[0] + if len(name) != 0: + + root = ET.Element('root') + wdir = ET.SubElement(root, 'wDir') + dhm = ET.SubElement(root, 'DHM') + release = ET.SubElement(root, 'Release') + infra = ET.SubElement(root, 'Infrastructure') + forest = ET.SubElement(root, 'Forest') + + wdir.set('Directory', 'Working') + dhm.set('Directory', 'DHM') + release.set('Directory', 'Release') + infra.set('Directory', 'Infrastructure') + forest.set('Directory', 'Forest') + + wdir.text = self.ui.wDir_lineEdit.text() + dhm.text = self.ui.DEM_lineEdit.text() + release.text = self.ui.release_lineEdit.text() + infra.text = self.ui.infra_lineEdit.text() + forest.text = self.ui.forest_lineEdit.text() + + tree = ET.ElementTree(root) + tree.write(name) + + def load(self): + xml_file = QFileDialog.getOpenFileNames(None, 'Open xml', + self.directory, + "xml (*.xml);;All Files (*.*)")[0] + + if len(xml_file) != 0: + tree = ET.parse(xml_file[0]) + root = tree.getroot() + + try: + self.ui.wDir_lineEdit.setText(root[0].text) + self.directory = root[0].text + except: + print("No Working Directory Path in File!") + + try: + self.ui.DEM_lineEdit.setText(root[1].text) + except: + print("No DEM Path in File!") + + try: + self.ui.release_lineEdit.setText(root[2].text) + except: + print("No Release Path in File!") + + try: + self.ui.infra_lineEdit.setText(root[3].text) + except: + print("No Infrastructure Path in File!") + + try: + self.ui.forest_lineEdit.setText(root[4].text) + except: + print("No Forest Path in File!") + + def quit(self): + QCoreApplication.quit() + + def open_wDir(self): + """Open the Working Directory, where results are stored""" + directory = QFileDialog.getExistingDirectory(None, 'Open Working Directory', + self.directory, + QFileDialog.ShowDirsOnly) + if len(directory) != 0: + self.directory = directory + self.ui.wDir_lineEdit.setText(self.directory) + + def open_dhm(self): + """Open digital elevation model""" + dem_file = QFileDialog.getOpenFileNames(None, 'Open DEM', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(dem_file[0]) != 0: + dem = dem_file[0] + self.ui.DEM_lineEdit.setText(dem[0]) + + def open_release(self): + """Open release layer""" + release_file = QFileDialog.getOpenFileNames(None, 'Open Release', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(release_file[0]) != 0: + release = release_file[0] + self.ui.release_lineEdit.setText(release[0]) + + def open_infra(self): + """Open infrastructure layer""" + infra_file = QFileDialog.getOpenFileNames(None, 'Open Infrastructure Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(infra_file[0]) != 0: + infra = infra_file[0] + self.ui.infra_lineEdit.setText(infra[0]) + + def open_forest(self): + """Open forest layer""" + forest_file = QFileDialog.getOpenFileNames(None, 'Open Forest Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + forest = forest_file[0] + self.ui.forest_lineEdit.setText(forest[0]) + + def update_progressBar(self, float, thread, start, end): + self.thread_list[thread] = float + self.start_list[thread] = start + self.end_list[thread] = end + + self.progress_value = sum(self.thread_list) / len(self.thread_list) + self.progressBar.setValue(self.progress_value) + for i in range(len(self.thread_list)): + sys.stdout.write( + "Thread {}: Startcell {} of {} = {}%"'\n'.format(i + 1, self.start_list[i], self.end_list[i], + self.thread_list[i])) + sys.stdout.flush() + for i in range(len(self.thread_list)): + sys.stdout.write('\x1b[1A' + '\x1b[2K') # Go 1 line up and erase that line + + @staticmethod + def showdialog(path): + msg = QMessageBox() + msg.setIcon(QMessageBox.Critical) + msg.setText("No " + path + " set") + msg.setWindowTitle("Error") + msg.setStandardButtons(QMessageBox.Ok) + msg.exec_() + + def calculation(self): + self.start = datetime.now().replace(microsecond=0) + self.calc_bool = False + + # Check if input is ok + if self.ui.wDir_lineEdit.text() == '': + self.showdialog('Working Directory') + return + if self.ui.DEM_lineEdit.text() == '': + self.showdialog('DEM Layer') + return + if self.ui.release_lineEdit.text() == '': + self.showdialog('Release Layer') + return + # Disable all input line Edits and Buttons + self.set_gui_bool(False) + + # Create result directory + time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + try: + os.makedirs(self.ui.wDir_lineEdit.text() + '/result/') + self.res_dir = ('/result/') + except FileExistsError: + self.res_dir = ('/result/') + + # Setup logger + + for handler in logging.root.handlers[:]: + logging.root.removeHandler(handler) + + logging.basicConfig(level=logging.INFO, + format='%(asctime)s %(levelname)-8s %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + filename=(self.directory + self.res_dir + 'log_{}.txt').format(time_string), + filemode='w') + + # Start of Calculation + logging.info('Start Calculation') + # Read in raster files + try: + dem, header = io.read_raster(self.ui.DEM_lineEdit.text()) + logging.info('DEM File: {}'.format(self.ui.DEM_lineEdit.text())) + except FileNotFoundError: + print("Wrong filepath or filename") + self.set_gui_bool(True) + return + + try: + release, release_header = io.read_raster(self.ui.release_lineEdit.text()) + logging.info('Release File: {}'.format(self.ui.release_lineEdit.text())) + except FileNotFoundError: + print("Wrong filepath or filename") + self.set_gui_bool(True) + return + + # Check if Layers have same size!!! + if header['ncols'] == release_header['ncols'] and header['nrows'] == release_header['nrows']: + print("DEM and Release Layer ok!") + else: + print("Error: Release Layer doesn't match DEM!") + self.set_gui_bool(True) + return + + try: + infra, infra_header = io.read_raster(self.ui.infra_lineEdit.text()) + if header['ncols'] == infra_header['ncols'] and header['nrows'] == infra_header['nrows']: + print("Infra Layer ok!") + self.calc_bool = True + logging.info('Infrastructure File: {}'.format(self.ui.infra_lineEdit.text())) + else: + print("Error: Infra Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(self.ui.forest_lineEdit.text()) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + self.prot_for_bool = True + logging.info('Forest File: {}'.format(self.forest_lineEdit.text())) + else: + print("Error: Forest Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + forest = np.zeros_like(dem) + + logging.info('Files read in') + + alpha = self.ui.alpha_Edit.text() + exp = self.ui.exp_Edit.text() + flux_threshold = self.ui.flux_Edit.text() + max_z = self.ui.z_Edit.text() + + logging.info('Alpha Angle: {}'.format(alpha)) + logging.info('Exponent: {}'.format(exp)) + logging.info('Flux Threshold: {}'.format(flux_threshold)) + logging.info('Max Z_delta: {}'.format(max_z)) + logging.info + + self.z_delta = np.zeros_like(dem) + self.flux = np.zeros_like(dem) + self.cell_counts = np.zeros_like(dem) + self.z_delta_sum = np.zeros_like(dem) + self.backcalc = np.zeros_like(dem) + self.fp_ta = np.zeros_like(dem) + self.sl_ta = np.zeros_like(dem) + + # Calculation + self.calc_class = Sim.Simulation(dem, header, release, release_header, infra, forest, self.calc_bool, alpha, exp, flux_threshold, max_z) + self.calc_class.value_changed.connect(self.update_progressBar) + self.calc_class.finished.connect(self.thread_finished) + logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) + self.calc_class.start() + + def thread_finished(self, z_delta, flux, count_array, z_delta_sum, backcalc, fp_ta, sl_ta): + logging.info('Calculation finished, getting results.') + for i in range(len(z_delta)): + self.z_delta = np.maximum(self.z_delta, z_delta[i]) + self.flux = np.maximum(self.flux, flux[i]) + self.cell_counts += count_array[i] + self.z_delta_sum += z_delta_sum[i] + self.backcalc = np.maximum(self.backcalc, backcalc[i]) + self.fp_ta = np.maximum(self.fp_ta, fp_ta[i]) + self.sl_ta = np.maximum(self.sl_ta, sl_ta[i]) + self.output() + + def output(self): + # time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + logging.info('Writing Output Files') + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "flux{}".format(self.ui.outputBox.currentText()), + self.flux) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "z_delta{}".format(self.ui.outputBox.currentText()), + self.z_delta) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "FP_travel_angle{}".format(self.ui.outputBox.currentText()), + self.fp_ta) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "SL_travel_angle{}".format(self.ui.outputBox.currentText()), + self.sl_ta) + if not self.calc_bool: + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "cell_counts{}".format(self.ui.outputBox.currentText()), + self.cell_counts) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "z_delta_sum{}".format(self.ui.outputBox.currentText()), + self.z_delta_sum) + if self.calc_bool: + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "backcalculation{}".format(self.ui.outputBox.currentText()), + self.backcalc) + + print("Calculation finished") + end = datetime.now().replace(microsecond=0) + logging.info('Calculation needed: ' + str(end - self.start) + ' seconds') + + # Handle GUI + #self.ui.progressBar.setValue(100) + self.set_gui_bool(True) + + +def main(args, kwargs): + + + alpha = args[0] + exp = args[1] + directory = args[2] + dem_path = args[3] + release_path = args[4] + if 'infra' in kwargs: + infra_path = kwargs.get('infra') + #print(infra_path) + else: + infra_path = None + + if 'forest' in kwargs: + forest_path = kwargs.get('forest') + #print(forest_path) + else: + forest_path = None + + if 'flux' in kwargs: + flux_threshold = float(kwargs.get('flux')) + #print(flux_threshold) + else: + flux_threshold = 3 * 10 ** -4 + + if 'max_z' in kwargs: + max_z = kwargs.get('max_z') + #print(max_z) + else: + max_z = 8848 + # Recomendet values: + # Avalanche = 270 + # Rockfall = 50 + # Soil Slide = 12 + + print("Starting...") + print("...") + + start = datetime.now().replace(microsecond=0) + calc_bool = False + # Create result directory + time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + try: + os.makedirs(directory + '/result/') + res_dir = ('/result/') + except FileExistsError: + res_dir = ('/result/') + + # Setup logger + for handler in logging.root.handlers[:]: + logging.root.removeHandler(handler) + + logging.basicConfig(level=logging.INFO, + format='%(asctime)s %(levelname)-8s %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + filename=(directory + res_dir + 'log_{}.txt').format(time_string), + filemode='w') + + # Start of Calculation + logging.info('Start Calculation') + logging.info('Alpha Angle: {}'.format(alpha)) + logging.info('Exponent: {}'.format(exp)) + logging.info('Flux Threshold: {}'.format(flux_threshold)) + logging.info('Max Z_delta: {}'.format(max_z)) + logging.info + # Read in raster files + try: + dem, header = io.read_raster(dem_path) + logging.info('DEM File: {}'.format(dem_path)) + except FileNotFoundError: + print("DEM: Wrong filepath or filename") + return + + try: + release, release_header = io.read_raster(release_path) + logging.info('Release File: {}'.format(release_path)) + except FileNotFoundError: + print("Wrong filepath or filename") + return + + # Check if Layers have same size!!! + if header['ncols'] == release_header['ncols'] and header['nrows'] == release_header['nrows']: + print("DEM and Release Layer ok!") + else: + print("Error: Release Layer doesn't match DEM!") + return + + try: + infra, infra_header = io.read_raster(infra_path) + if header['ncols'] == infra_header['ncols'] and header['nrows'] == infra_header['nrows']: + print("Infra Layer ok!") + calc_bool = True + logging.info('Infrastructure File: {}'.format(infra_path)) + else: + print("Error: Infra Layer doesn't match DEM!") + return + except: + infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(forest_path) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + prot_for_bool = True + logging.info('Forest File: {}'.format(forest_path)) + else: + print("Error: Forest Layer doesn't match DEM!") + return + except: + forest = np.zeros_like(dem) + + logging.info('Files read in') + + z_delta = np.zeros_like(dem) + flux = np.zeros_like(dem) + cell_counts = np.zeros_like(dem) + z_delta_sum = np.zeros_like(dem) + backcalc = np.zeros_like(dem) + fp_ta = np.zeros_like(dem) + sl_ta = np.zeros_like(dem) + + avaiable_memory = psutil.virtual_memory()[1] + needed_memory = sys.getsizeof(dem) + + max_number_procces = int(avaiable_memory / (needed_memory * 10)) + + print( + "There are {} Bytes of Memory avaiable and {} Bytes needed per process. \nMax. Nr. of Processes = {}".format( + avaiable_memory, needed_memory*10, max_number_procces)) + + # Calculation + logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) + + if calc_bool: + release_list = fc.split_release(release, release_header, min(mp.cpu_count() * 2, max_number_procces)) + + print("{} Processes started.".format(len(release_list))) + pool = mp.Pool(len(release_list)) + results = pool.map(fc.calculation, + [[dem, header, infra, forest, release_pixel, alpha, exp, flux_threshold, max_z] + for release_pixel in release_list]) + pool.close() + pool.join() + else: + release_list = fc.split_release(release, release_header, min(mp.cpu_count() * 4, max_number_procces)) + + print("{} Processes started.".format(len(release_list))) + pool = mp.Pool(mp.cpu_count()) + # results = pool.map(gc.calculation, iterable) + results = pool.map(fc.calculation_effect, + [[dem, header, forest, release_pixel, alpha, exp, flux_threshold, max_z] for + release_pixel in release_list]) + pool.close() + pool.join() + + z_delta_list = [] + flux_list = [] + cc_list = [] + z_delta_sum_list = [] + backcalc_list = [] + fp_ta_list = [] + sl_ta_list = [] + for i in range(len(results)): + res = results[i] + res = list(res) + z_delta_list.append(res[0]) + flux_list.append(res[1]) + cc_list.append(res[2]) + z_delta_sum_list.append(res[3]) + backcalc_list.append(res[4]) + fp_ta_list.append(res[5]) + sl_ta_list.append(res[6]) + + logging.info('Calculation finished, getting results.') + for i in range(len(z_delta_list)): + z_delta = np.maximum(z_delta, z_delta_list[i]) + flux = np.maximum(flux, flux_list[i]) + cell_counts += cc_list[i] + z_delta_sum += z_delta_sum_list[i] + backcalc = np.maximum(backcalc, backcalc_list[i]) + fp_ta = np.maximum(fp_ta, fp_ta_list[i]) + sl_ta = np.maximum(sl_ta, sl_ta_list[i]) + + + # time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + logging.info('Writing Output Files') + output_format = '.tif' + io.output_raster(dem_path, + directory + res_dir + "flux{}".format(output_format), + flux) + io.output_raster(dem_path, + directory + res_dir + "z_delta{}".format(output_format), + z_delta) + io.output_raster(dem_path, + directory + res_dir + "FP_travel_angle{}".format(output_format), + fp_ta) + io.output_raster(dem_path, + directory + res_dir + "SL_travel_angle{}".format(output_format), + sl_ta) + if not calc_bool: # if no infra + io.output_raster(dem_path, + directory + res_dir + "cell_counts{}".format(output_format), + cell_counts) + io.output_raster(dem_path, + directory + res_dir + "z_delta_sum{}".format( output_format), + z_delta_sum) + if calc_bool: # if infra + io.output_raster(dem_path, + directory + res_dir + "backcalculation{}".format(output_format), + backcalc) + + print("Calculation finished") + print("...") + end = datetime.now().replace(microsecond=0) + logging.info('Calculation needed: ' + str(end - start) + ' seconds') + + +if __name__ == '__main__': + #mp.set_start_method('spawn') # used in Windows + argv = sys.argv[1:] + #argv = ["28", "8", "../runs/open_slope/", "../runs/open_slope/parabola.asc", "../runs/open_slope/release.tif", "forest=../runs/open_slope/forest_deep.tif", "flux=0.0003" , "max_z=270"] #"forest=../runs/open_slope/forest.tif", + #argv = ['--gui'] + if len(argv) < 1: + print("Too few input arguments!!!") + sys.exit(1) + if len(argv) == 1 and argv[0] == '--gui': + Flow_Py_EXEC() + else: + args=[arg for arg in argv if arg.find('=')<0] + kwargs={kw[0]:kw[1] for kw in [ar.split('=') for ar in argv if ar.find('=')>0]} + + main(args, kwargs) + +# example dam: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif +# example dam/infra: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif infra=./examples/dam/infra.tif flux=0.003 max_z=270 +# example dam/forest: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif forest=./examples/dam/infra.tif flux=0.003 max_z=270 diff --git a/Flow_GUI.py b/Flow_GUI.py index d8cf211..4918bc4 100644 --- a/Flow_GUI.py +++ b/Flow_GUI.py @@ -2,71 +2,77 @@ # Form implementation generated from reading ui file '.\Flow_GUI.ui' # -# Created by: PyQt5 UI code generator 5.9.2 +# Created by: PyQt5 UI code generator 5.15.5 # -# WARNING! All changes made in this file will be lost! +# WARNING: Any manual changes made to this file will be lost when pyuic5 is +# run again. Do not edit this file unless you know what you are doing. + from PyQt5 import QtCore, QtGui, QtWidgets + class Ui_MainWindow(object): def setupUi(self, MainWindow): MainWindow.setObjectName("MainWindow") MainWindow.resize(546, 352) icon = QtGui.QIcon() - icon.addPixmap(QtGui.QPixmap("FlowPylogo.png"), QtGui.QIcon.Normal, QtGui.QIcon.Off) + icon.addPixmap(QtGui.QPixmap(".\\FlowPylogo.png"), QtGui.QIcon.Normal, QtGui.QIcon.Off) MainWindow.setWindowIcon(icon) self.centralwidget = QtWidgets.QWidget(MainWindow) self.centralwidget.setObjectName("centralwidget") self.gridLayout = QtWidgets.QGridLayout(self.centralwidget) self.gridLayout.setObjectName("gridLayout") + self.forest_Label = QtWidgets.QLabel(self.centralwidget) + self.forest_Label.setObjectName("forest_Label") + self.gridLayout.addWidget(self.forest_Label, 6, 1, 1, 1) self.label = QtWidgets.QLabel(self.centralwidget) self.label.setObjectName("label") - self.gridLayout.addWidget(self.label, 11, 1, 1, 1) + self.gridLayout.addWidget(self.label, 12, 1, 1, 1) self.line = QtWidgets.QFrame(self.centralwidget) self.line.setFrameShape(QtWidgets.QFrame.HLine) self.line.setFrameShadow(QtWidgets.QFrame.Sunken) self.line.setObjectName("line") self.gridLayout.addWidget(self.line, 4, 1, 1, 1) - self.infra_lineEdit = QtWidgets.QLineEdit(self.centralwidget) - self.infra_lineEdit.setObjectName("infra_lineEdit") - self.gridLayout.addWidget(self.infra_lineEdit, 5, 2, 1, 1) - self.wDir_Button = QtWidgets.QToolButton(self.centralwidget) - self.wDir_Button.setCheckable(False) - self.wDir_Button.setObjectName("wDir_Button") - self.gridLayout.addWidget(self.wDir_Button, 1, 3, 1, 1) + self.line_5 = QtWidgets.QFrame(self.centralwidget) + self.line_5.setFrameShape(QtWidgets.QFrame.HLine) + self.line_5.setFrameShadow(QtWidgets.QFrame.Sunken) + self.line_5.setObjectName("line_5") + self.gridLayout.addWidget(self.line_5, 8, 3, 1, 1) self.line_4 = QtWidgets.QFrame(self.centralwidget) self.line_4.setFrameShape(QtWidgets.QFrame.HLine) self.line_4.setFrameShadow(QtWidgets.QFrame.Sunken) self.line_4.setObjectName("line_4") - self.gridLayout.addWidget(self.line_4, 4, 2, 1, 1) + self.gridLayout.addWidget(self.line_4, 4, 3, 1, 1) self.Release_Button = QtWidgets.QToolButton(self.centralwidget) self.Release_Button.setObjectName("Release_Button") - self.gridLayout.addWidget(self.Release_Button, 3, 3, 1, 1) + self.gridLayout.addWidget(self.Release_Button, 3, 4, 1, 1) self.calc_Button = QtWidgets.QPushButton(self.centralwidget) self.calc_Button.setObjectName("calc_Button") - self.gridLayout.addWidget(self.calc_Button, 12, 2, 1, 1) - self.release_label = QtWidgets.QLabel(self.centralwidget) - self.release_label.setObjectName("release_label") - self.gridLayout.addWidget(self.release_label, 3, 1, 1, 1) - self.wDir_label = QtWidgets.QLabel(self.centralwidget) - self.wDir_label.setObjectName("wDir_label") - self.gridLayout.addWidget(self.wDir_label, 1, 1, 1, 1) + self.gridLayout.addWidget(self.calc_Button, 13, 3, 1, 1) self.infra_Button = QtWidgets.QToolButton(self.centralwidget) self.infra_Button.setObjectName("infra_Button") - self.gridLayout.addWidget(self.infra_Button, 5, 3, 1, 1) + self.gridLayout.addWidget(self.infra_Button, 5, 4, 1, 1) self.line_6 = QtWidgets.QFrame(self.centralwidget) self.line_6.setFrameShape(QtWidgets.QFrame.HLine) self.line_6.setFrameShadow(QtWidgets.QFrame.Sunken) self.line_6.setObjectName("line_6") - self.gridLayout.addWidget(self.line_6, 10, 2, 1, 1) + self.gridLayout.addWidget(self.line_6, 11, 3, 1, 1) + self.infra_lineEdit = QtWidgets.QLineEdit(self.centralwidget) + self.infra_lineEdit.setObjectName("infra_lineEdit") + self.gridLayout.addWidget(self.infra_lineEdit, 5, 3, 1, 1) + self.wDir_Button = QtWidgets.QToolButton(self.centralwidget) + self.wDir_Button.setCheckable(False) + self.wDir_Button.setObjectName("wDir_Button") + self.gridLayout.addWidget(self.wDir_Button, 1, 4, 1, 1) + self.release_label = QtWidgets.QLabel(self.centralwidget) + self.release_label.setObjectName("release_label") + self.gridLayout.addWidget(self.release_label, 3, 1, 1, 1) + self.wDir_label = QtWidgets.QLabel(self.centralwidget) + self.wDir_label.setObjectName("wDir_label") + self.gridLayout.addWidget(self.wDir_label, 1, 1, 1, 1) self.DEM_Button = QtWidgets.QToolButton(self.centralwidget) self.DEM_Button.setObjectName("DEM_Button") - self.gridLayout.addWidget(self.DEM_Button, 2, 3, 1, 1) - self.line_5 = QtWidgets.QFrame(self.centralwidget) - self.line_5.setFrameShape(QtWidgets.QFrame.HLine) - self.line_5.setFrameShadow(QtWidgets.QFrame.Sunken) - self.line_5.setObjectName("line_5") - self.gridLayout.addWidget(self.line_5, 7, 2, 1, 1) + self.gridLayout.addWidget(self.DEM_Button, 2, 4, 1, 1) self.horizontalLayout = QtWidgets.QHBoxLayout() self.horizontalLayout.setObjectName("horizontalLayout") self.label_7 = QtWidgets.QLabel(self.centralwidget) @@ -92,45 +98,32 @@ def setupUi(self, MainWindow): self.horizontalLayout.addWidget(self.exp_Edit) spacerItem1 = QtWidgets.QSpacerItem(40, 20, QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Minimum) self.horizontalLayout.addItem(spacerItem1) - self.gridLayout.addLayout(self.horizontalLayout, 8, 2, 1, 1) + self.gridLayout.addLayout(self.horizontalLayout, 9, 3, 1, 1) + self.DEM_label = QtWidgets.QLabel(self.centralwidget) + self.DEM_label.setObjectName("DEM_label") + self.gridLayout.addWidget(self.DEM_label, 2, 1, 1, 1) + self.DEM_lineEdit = QtWidgets.QLineEdit(self.centralwidget) + self.DEM_lineEdit.setObjectName("DEM_lineEdit") + self.gridLayout.addWidget(self.DEM_lineEdit, 2, 3, 1, 1) + self.infra_label = QtWidgets.QLabel(self.centralwidget) + self.infra_label.setObjectName("infra_label") + self.gridLayout.addWidget(self.infra_label, 5, 1, 1, 1) + self.wDir_lineEdit = QtWidgets.QLineEdit(self.centralwidget) + self.wDir_lineEdit.setObjectName("wDir_lineEdit") + self.gridLayout.addWidget(self.wDir_lineEdit, 1, 3, 1, 1) + self.line_3 = QtWidgets.QFrame(self.centralwidget) + self.line_3.setFrameShape(QtWidgets.QFrame.HLine) + self.line_3.setFrameShadow(QtWidgets.QFrame.Sunken) + self.line_3.setObjectName("line_3") + self.gridLayout.addWidget(self.line_3, 11, 1, 1, 1) self.label_2 = QtWidgets.QLabel(self.centralwidget) self.label_2.setObjectName("label_2") - self.gridLayout.addWidget(self.label_2, 8, 1, 1, 1) + self.gridLayout.addWidget(self.label_2, 9, 1, 1, 1) self.line_2 = QtWidgets.QFrame(self.centralwidget) self.line_2.setFrameShape(QtWidgets.QFrame.HLine) self.line_2.setFrameShadow(QtWidgets.QFrame.Sunken) self.line_2.setObjectName("line_2") - self.gridLayout.addWidget(self.line_2, 7, 1, 1, 1) - self.release_lineEdit = QtWidgets.QLineEdit(self.centralwidget) - self.release_lineEdit.setObjectName("release_lineEdit") - self.gridLayout.addWidget(self.release_lineEdit, 3, 2, 1, 1) - self.line_3 = QtWidgets.QFrame(self.centralwidget) - self.line_3.setFrameShape(QtWidgets.QFrame.HLine) - self.line_3.setFrameShadow(QtWidgets.QFrame.Sunken) - self.line_3.setObjectName("line_3") - self.gridLayout.addWidget(self.line_3, 10, 1, 1, 1) - self.horizontalLayout_2 = QtWidgets.QHBoxLayout() - self.horizontalLayout_2.setObjectName("horizontalLayout_2") - self.outputBox = QtWidgets.QComboBox(self.centralwidget) - self.outputBox.setObjectName("outputBox") - self.outputBox.addItem("") - self.outputBox.addItem("") - self.horizontalLayout_2.addWidget(self.outputBox) - spacerItem2 = QtWidgets.QSpacerItem(40, 20, QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Minimum) - self.horizontalLayout_2.addItem(spacerItem2) - self.gridLayout.addLayout(self.horizontalLayout_2, 11, 2, 1, 1) - self.infra_label = QtWidgets.QLabel(self.centralwidget) - self.infra_label.setObjectName("infra_label") - self.gridLayout.addWidget(self.infra_label, 5, 1, 1, 1) - self.wDir_lineEdit = QtWidgets.QLineEdit(self.centralwidget) - self.wDir_lineEdit.setObjectName("wDir_lineEdit") - self.gridLayout.addWidget(self.wDir_lineEdit, 1, 2, 1, 1) - self.DEM_label = QtWidgets.QLabel(self.centralwidget) - self.DEM_label.setObjectName("DEM_label") - self.gridLayout.addWidget(self.DEM_label, 2, 1, 1, 1) - self.DEM_lineEdit = QtWidgets.QLineEdit(self.centralwidget) - self.DEM_lineEdit.setObjectName("DEM_lineEdit") - self.gridLayout.addWidget(self.DEM_lineEdit, 2, 2, 1, 1) + self.gridLayout.addWidget(self.line_2, 8, 1, 1, 1) self.horizontalLayout_3 = QtWidgets.QHBoxLayout() self.horizontalLayout_3.setObjectName("horizontalLayout_3") self.label_4 = QtWidgets.QLabel(self.centralwidget) @@ -139,20 +132,39 @@ def setupUi(self, MainWindow): self.flux_Edit = QtWidgets.QLineEdit(self.centralwidget) self.flux_Edit.setObjectName("flux_Edit") self.horizontalLayout_3.addWidget(self.flux_Edit) - spacerItem3 = QtWidgets.QSpacerItem(40, 20, QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Minimum) - self.horizontalLayout_3.addItem(spacerItem3) + spacerItem2 = QtWidgets.QSpacerItem(40, 20, QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Minimum) + self.horizontalLayout_3.addItem(spacerItem2) self.label_3 = QtWidgets.QLabel(self.centralwidget) self.label_3.setObjectName("label_3") self.horizontalLayout_3.addWidget(self.label_3) self.z_Edit = QtWidgets.QLineEdit(self.centralwidget) self.z_Edit.setObjectName("z_Edit") self.horizontalLayout_3.addWidget(self.z_Edit) + spacerItem3 = QtWidgets.QSpacerItem(40, 20, QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Minimum) + self.horizontalLayout_3.addItem(spacerItem3) + self.gridLayout.addLayout(self.horizontalLayout_3, 10, 3, 1, 1) + self.release_lineEdit = QtWidgets.QLineEdit(self.centralwidget) + self.release_lineEdit.setObjectName("release_lineEdit") + self.gridLayout.addWidget(self.release_lineEdit, 3, 3, 1, 1) + self.horizontalLayout_2 = QtWidgets.QHBoxLayout() + self.horizontalLayout_2.setObjectName("horizontalLayout_2") + self.outputBox = QtWidgets.QComboBox(self.centralwidget) + self.outputBox.setObjectName("outputBox") + self.outputBox.addItem("") + self.outputBox.addItem("") + self.horizontalLayout_2.addWidget(self.outputBox) spacerItem4 = QtWidgets.QSpacerItem(40, 20, QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Minimum) - self.horizontalLayout_3.addItem(spacerItem4) - self.gridLayout.addLayout(self.horizontalLayout_3, 9, 2, 1, 1) + self.horizontalLayout_2.addItem(spacerItem4) + self.gridLayout.addLayout(self.horizontalLayout_2, 12, 3, 1, 1) + self.forest_lineEdit = QtWidgets.QLineEdit(self.centralwidget) + self.forest_lineEdit.setObjectName("forest_lineEdit") + self.gridLayout.addWidget(self.forest_lineEdit, 6, 3, 1, 1) + self.forest_Button = QtWidgets.QToolButton(self.centralwidget) + self.forest_Button.setObjectName("forest_Button") + self.gridLayout.addWidget(self.forest_Button, 6, 4, 1, 1) MainWindow.setCentralWidget(self.centralwidget) self.menubar = QtWidgets.QMenuBar(MainWindow) - self.menubar.setGeometry(QtCore.QRect(0, 0, 546, 26)) + self.menubar.setGeometry(QtCore.QRect(0, 0, 546, 21)) self.menubar.setDefaultUp(False) self.menubar.setObjectName("menubar") self.menuFile = QtWidgets.QMenu(self.menubar) @@ -181,23 +193,25 @@ def setupUi(self, MainWindow): def retranslateUi(self, MainWindow): _translate = QtCore.QCoreApplication.translate MainWindow.setWindowTitle(_translate("MainWindow", "Flow-Py")) + self.forest_Label.setText(_translate("MainWindow", "Forest")) self.label.setText(_translate("MainWindow", "Output Format")) - self.wDir_Button.setText(_translate("MainWindow", "...")) self.Release_Button.setText(_translate("MainWindow", "...")) self.calc_Button.setText(_translate("MainWindow", "Calculate")) + self.infra_Button.setText(_translate("MainWindow", "...")) + self.wDir_Button.setText(_translate("MainWindow", "...")) self.release_label.setText(_translate("MainWindow", "Release Layer")) self.wDir_label.setText(_translate("MainWindow", "Working Directory")) - self.infra_Button.setText(_translate("MainWindow", "...")) self.DEM_Button.setText(_translate("MainWindow", "...")) self.label_7.setText(_translate("MainWindow", "alpha")) self.label_6.setText(_translate("MainWindow", "exp ")) - self.label_2.setText(_translate("MainWindow", "Parameters")) - self.outputBox.setItemText(0, _translate("MainWindow", ".tif")) - self.outputBox.setItemText(1, _translate("MainWindow", ".asc")) - self.infra_label.setText(_translate("MainWindow", "Infrastructure")) self.DEM_label.setText(_translate("MainWindow", "DEM Layer")) + self.infra_label.setText(_translate("MainWindow", "Infrastructure")) + self.label_2.setText(_translate("MainWindow", "Parameters")) self.label_4.setText(_translate("MainWindow", "flux ")) self.label_3.setText(_translate("MainWindow", "max_z")) + self.outputBox.setItemText(0, _translate("MainWindow", ".tif")) + self.outputBox.setItemText(1, _translate("MainWindow", ".asc")) + self.forest_Button.setText(_translate("MainWindow", "...")) self.menuFile.setTitle(_translate("MainWindow", "File")) self.toolBar.setWindowTitle(_translate("MainWindow", "toolBar")) self.actionSave.setText(_translate("MainWindow", "Save")) @@ -213,4 +227,3 @@ def retranslateUi(self, MainWindow): ui.setupUi(MainWindow) MainWindow.show() sys.exit(app.exec_()) - diff --git a/Flow_GUI.ui b/Flow_GUI.ui index 2d1b027..9eeeed5 100755 --- a/Flow_GUI.ui +++ b/Flow_GUI.ui @@ -19,7 +19,14 @@ - + + + + Forest + + + + Output Format @@ -33,83 +40,83 @@ - - - - - - - ... - - - false + + + + Qt::Horizontal - + Qt::Horizontal - + ... - + Calculate - - + + - Release Layer + ... - - - - Working Directory + + + + Qt::Horizontal - + + + + ... + + false + - - - - Qt::Horizontal + + + + Release Layer - - + + - ... + Working Directory - - - - Qt::Horizontal + + + + ... - + @@ -182,82 +189,48 @@ - - + + - Parameters + DEM Layer - - - - Qt::Horizontal + + + + + + + Infrastructure - - + + - + Qt::Horizontal - - - - - - - .tif - - - - - .asc - - - - - - - - Qt::Horizontal - - - - 40 - 20 - - - - - - - - + + - Infrastructure + Parameters - - - - - - - DEM Layer + + + + Qt::Horizontal - - - - + @@ -307,6 +280,50 @@ + + + + + + + + + + .tif + + + + + .asc + + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + + + + + + ... + + + @@ -315,7 +332,7 @@ 0 0 546 - 26 + 21 diff --git a/Simulation.py b/Simulation.py index 911e1f4..b916f0b 100644 --- a/Simulation.py +++ b/Simulation.py @@ -37,13 +37,14 @@ class Simulation(QThread): value_changed = pyqtSignal(float) finished = pyqtSignal(list, list, list, list, list, list, list) - def __init__(self, dem, header, release, release_header, infra, calc_bool, alpha, exp, flux, max_z): + def __init__(self, dem, header, release, release_header, infra, forest, calc_bool, alpha, exp, flux, max_z): QThread.__init__(self) self.dem = dem self.header = header self.release = release self.release_header = release_header self.infra = infra + self.forest = forest self.numberofprocesses = mp.cpu_count() self.calc_bool = calc_bool self.alpha = alpha @@ -76,7 +77,7 @@ def run(self): print("{} Processes started.".format(len(release_list))) pool = mp.Pool(len(release_list)) - results = pool.map(fc.calculation, [[self.dem, self.header, self.infra, release_pixel, self.alpha, self.exp, self.flux, self.max_z] for release_pixel in release_list]) + results = pool.map(fc.calculation, [[self.dem, self.header, self.infra, self.forest, release_pixel, self.alpha, self.exp, self.flux, self.max_z] for release_pixel in release_list]) pool.close() pool.join() else: @@ -86,7 +87,7 @@ def run(self): pool = mp.Pool(mp.cpu_count()) #results = pool.map(gc.calculation, iterable) results = pool.map(fc.calculation_effect, - [[self.dem, self.header, release_pixel, + [[self.dem, self.header, self.forest, release_pixel, self.alpha, self.exp, self.flux, self.max_z] for release_pixel in release_list]) pool.close() diff --git a/flow_class.py b/flow_class.py index 0a74e4c..9aa5a40 100644 --- a/flow_class.py +++ b/flow_class.py @@ -28,13 +28,14 @@ class Cell: - def __init__(self, rowindex, colindex, dem_ng, cellsize, flux, z_delta, parent, alpha, exp, flux_threshold, max_z_delta, startcell): + def __init__(self, rowindex, colindex, dem_ng, forest, cellsize, flux, z_delta, parent, alpha, exp, flux_threshold, max_z_delta, startcell): '''This class handles the spreading over the DEM! Depending on the process different alpha angles are used for energy dissipation.''' self.rowindex = rowindex self.colindex = colindex self.altitude = dem_ng[1, 1] self.dem_ng = dem_ng + self.forest = forest self.cellsize = cellsize self.tan_beta = np.zeros_like(self.dem_ng) self.dist = np.zeros_like(self.dem_ng) @@ -51,13 +52,21 @@ def __init__(self, rowindex, colindex, dem_ng, cellsize, flux, z_delta, parent, self.max_distance = 0 self.min_gamma = 0 self.max_gamma = 0 - self.sl_gamma = 0 + self.sl_gamma = 0 + # Parameters for Forest Friction, right now use it just for avalanches + #self.alpha_forest = 10 # Max added friction angel + self.max_added_friction_forest = 20 # degrees added to friction angle + self.min_added_friction_forest = 2 # minimium effect forested terrain can have + self.no_friction_effect_v = 270 # velocity shared for friction and detrainment methods + self.max_added_detrainment_forest = 0.0003 # + self.min_added_detrainment_forest = 0.00001 + self.no_detrainmnet_effect_v = 270 # if type(startcell) == bool: # check, if start cell exist (start cell is release point) self.is_start = True # set is_start to True else: self.startcell = startcell # give startcell to cell - self.is_start = False # set is_start to False + self.is_start = False # s et is_start to False self.parent = [] if type(parent) == Cell: @@ -66,6 +75,20 @@ def __init__(self, rowindex, colindex, dem_ng, cellsize, flux, z_delta, parent, def add_os(self, flux): self.flux += flux + def forest_detrainment(self ): + """ + linear decrease of forest effect with regard to alpha increase and kinetic energy height + This is the detrainment routine for forest. It should reduce the routing flux of the avalanche. + self.max_added_detrainment_forest = .001 + self.min_added_detrainment_forest = 0 + self.no_detrainmnet_effect_v = 45 + """ + no_detrainmnet_effect_zdelta = self.no_detrainmnet_effect_v**(2)/(np.sqrt(2) * 9.8) # change veloctiy into kinetic energy line hight (z_delta) 9.8 = gravity. derrived from 1/2mv^2 = mgh + rest = self.max_added_detrainment_forest * self.forest # detrainment effect scalled to forest, should be zero for non-forested area + slope = (rest - self.min_added_detrainment_forest)/(0-no_detrainmnet_effect_zdelta) # rise over run (should be negative slope) + self.detrainment = max(self.min_added_detrainment_forest, slope * self.z_delta + rest) # y = mx + b, shere z_delta is the x + + def add_parent(self, parent): self.parent.append(parent) @@ -90,8 +113,22 @@ def calc_sl_travelangle(self): def calc_z_delta(self): self.z_delta_neighbour = np.zeros((3, 3)) self.z_gamma = self.altitude - self.dem_ng + no_friction_effect_zdelta = self.no_friction_effect_v**(2) / (np.sqrt(2) * 9.8) # change veloctiy into energy line hight (z_delta) 9.8 = gravity. derrived from 1/2mv^2 = mgh ds = np.array([[np.sqrt(2), 1, np.sqrt(2)], [1, 0, 1], [np.sqrt(2), 1, np.sqrt(2)]]) - tan_alpha = np.tan(np.deg2rad(self.alpha)) + ## Calculation for Forest Friction leads to new alpha_calc + if self.forest > 0: + if self.z_delta < no_friction_effect_zdelta: # no min_added forest values becuase of this line + rest = self.max_added_friction_forest * self.forest # friction at rest v=0 would be applied to start cells + slope = (rest - self.min_added_friction_forest) / (0 - no_friction_effect_zdelta) # rise over run + friction = max(self.min_added_friction_forest, + slope * self.z_delta + rest) # y = mx + b, shere z_delta is the x + alpha_calc = self.alpha + max(0, friction) + else: + alpha_calc = self.alpha + self.min_added_friction_forest + else: + alpha_calc = self.alpha + # Normal calculation + tan_alpha = np.tan(np.deg2rad(alpha_calc)) self.z_alpha = ds * self.cellsize * tan_alpha self.z_delta_neighbour = self.z_delta + self.z_gamma - self.z_alpha self.z_delta_neighbour[self.z_delta_neighbour < 0] = 0 @@ -178,7 +215,7 @@ def calc_persistence(self): # pers[1, 1] = 0 # self.persistence += pers * maxweight # ============================================================================= - + def calc_distribution(self): self.calc_z_delta() @@ -186,14 +223,21 @@ def calc_distribution(self): self.persistence *= self.no_flow self.calc_tanbeta() #print(self.persistence) + self.forest_detrainment() if not self.is_start: self.calc_fp_travelangle() self.calc_sl_travelangle() + self.flux = max(0.0003, self.flux - self.detrainment) # here we subtract the detrainment from the flux before moving flux to new cells. + threshold = self.flux_threshold - if np.sum(self.r_t) > 0: + if np.sum(self.r_t) > 0: # if there is routing flux self.dist = (self.persistence * self.r_t) / np.sum(self.persistence * self.r_t) * self.flux + #detrainment = forest_scale(self, max_detrainment, min_detrainment, FSI, max_FE_velocity, z_delta ) + #self.dist = self.dist - self.detrainment # todo max needed to keep it always positive + #todo make sure that I am only removing detrained mass once not every spatial iteration that includes a cell as a neighbor + #print(self.detrainment, "detrainment" , self.dist, "mass") # This lines handle if a distribution to a neighbour cell is lower then the threshold, so we don´t lose # flux. # The flux of this cells will then spread equally to all neighbour cells diff --git a/flow_core.py b/flow_core.py index c0a148f..45f10c1 100644 --- a/flow_core.py +++ b/flow_core.py @@ -137,8 +137,7 @@ def split_release(release, header_release, pieces): release_list.append(c) print("Release Split from {} to {}".format(breakpoint_x, i)) breakpoint_x = i - - + return release_list @@ -149,7 +148,7 @@ def calculation(args): Input parameters: dem The digital elevation model header The header of the elevation model - forest The forest layer + infras The infrastructure layer process Which process to calculate (Avalanche, Rockfall, SoilSlides) release The list of release arrays alpha @@ -158,8 +157,9 @@ def calculation(args): max_z_delta Output parameters: - elh Array like DEM with the max. Energy Line Height for every - pixel + z_delta_array Array like DEM with the max. Energy Line Height for every + pixel + z_delta_sum Array... mass_array Array with max. concentration factor saved count_array Array with the number of hits for every pixel elh_sum Array with the sum of Energy Line Height @@ -169,11 +169,12 @@ def calculation(args): dem = args[0] header = args[1] infra = args[2] - release = args[3] - alpha = args[4] - exp = args[5] - flux_threshold = args[6] - max_z_delta = args[7] + forest = args[3] + release = args[4] + alpha = args[5] + exp = args[6] + flux_threshold = args[7] + max_z_delta = args[8] #print(len(args), max_z_delta) z_delta_array = np.zeros_like(dem) @@ -289,11 +290,12 @@ def calculation_effect(args): dem = args[0] header = args[1] - release = args[2] - alpha = args[3] - exp = args[4] - flux_threshold = args[5] - max_z_delta = args[6] + forest = args[2] + release = args[3] + alpha = args[4] + exp = args[5] + flux_threshold = args[6] + max_z_delta = args[7] z_delta_array = np.zeros_like(dem) z_delta_sum = np.zeros_like(dem) @@ -325,7 +327,8 @@ def calculation_effect(args): startcell_idx += 1 continue - startcell = Cell(row_idx, col_idx, dem_ng, cellsize, 1, 0, None, + startcell = Cell(row_idx, col_idx, dem_ng, forest[row_idx, col_idx], + cellsize, 1, 0, None, alpha, exp, flux_threshold, max_z_delta, True) # If this is a startcell just give a Bool to startcell otherwise the object startcell @@ -358,7 +361,9 @@ def calculation_effect(args): if (nodata in dem_ng) or np.size(dem_ng) < 9: continue cell_list.append( - Cell(row[k], col[k], dem_ng, cellsize, flux[k], z_delta[k], cell, alpha, exp, flux_threshold, max_z_delta, startcell)) + Cell(row[k], col[k], dem_ng, forest[row[k], col[k]], + cellsize, flux[k], z_delta[k], cell, alpha, exp, + flux_threshold, max_z_delta, startcell)) for cell in cell_list: z_delta_array[cell.rowindex, cell.colindex] = max(z_delta_array[cell.rowindex, cell.colindex], cell.z_delta) diff --git a/main.py b/main.py index b532c52..d3e6be0 100755 --- a/main.py +++ b/main.py @@ -66,7 +66,7 @@ def __init__(self): self.ui.DEM_Button.clicked.connect(self.open_dhm) self.ui.Release_Button.clicked.connect(self.open_release) self.ui.infra_Button.clicked.connect(self.open_infra) - #self.ui.forest_Button.clicked.connect(self.open_forest) + self.ui.forest_Button.clicked.connect(self.open_forest) #self.ui.process_Box.currentIndexChanged.connect(self.processChanged) self.ui.calc_Button.clicked.connect(self.calculation) self.ui.actionSave.triggered.connect(self.save) @@ -96,6 +96,7 @@ def set_gui_bool(self, bool): self.ui.DEM_lineEdit.setEnabled(bool) self.ui.release_lineEdit.setEnabled(bool) self.ui.infra_lineEdit.setEnabled(bool) + self.ui.forest_lineEdit.setEnabled(bool) def save(self): """Save the input paths""" @@ -120,7 +121,7 @@ def save(self): dhm.text = self.ui.DEM_lineEdit.text() release.text = self.ui.release_lineEdit.text() infra.text = self.ui.infra_lineEdit.text() - #forest.text = self.ui.forest_lineEdit.text() + forest.text = self.ui.forest_lineEdit.text() tree = ET.ElementTree(root) tree.write(name) @@ -198,6 +199,14 @@ def open_infra(self): if len(infra_file[0]) != 0: infra = infra_file[0] self.ui.infra_lineEdit.setText(infra[0]) + + def open_forest(self): + """Open forest layer""" + forest_file = QFileDialog.getOpenFileNames(None, 'Open Forest Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + forest = forest_file[0] + self.ui.forest_lineEdit.setText(forest[0]) def update_progressBar(self, float, thread, start, end): self.thread_list[thread] = float @@ -298,6 +307,19 @@ def calculation(self): return except: infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(self.ui.forest_lineEdit.text()) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + self.prot_for_bool = True + logging.info('Forest File: {}'.format(self.forest_lineEdit.text())) + else: + print("Error: Forest Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + forest = np.zeros_like(dem) logging.info('Files read in') @@ -321,7 +343,7 @@ def calculation(self): self.sl_ta = np.zeros_like(dem) # Calculation - self.calc_class = Sim.Simulation(dem, header, release, release_header, infra, self.calc_bool, alpha, exp, flux_threshold, max_z) + self.calc_class = Sim.Simulation(dem, header, release, release_header, infra, forest, self.calc_bool, alpha, exp, flux_threshold, max_z) self.calc_class.value_changed.connect(self.update_progressBar) self.calc_class.finished.connect(self.thread_finished) logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) @@ -389,6 +411,12 @@ def main(args, kwargs): else: infra_path = None + if 'forest' in kwargs: + forest_path = kwargs.get('forest') + #print(forest_path) + else: + forest_path = None + if 'flux' in kwargs: flux_threshold = float(kwargs.get('flux')) #print(flux_threshold) @@ -468,6 +496,18 @@ def main(args, kwargs): return except: infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(forest_path) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + prot_for_bool = True + logging.info('Forest File: {}'.format(forest_path)) + else: + print("Error: Forest Layer doesn't match DEM!") + return + except: + forest = np.zeros_like(dem) logging.info('Files read in') @@ -497,7 +537,7 @@ def main(args, kwargs): print("{} Processes started.".format(len(release_list))) pool = mp.Pool(len(release_list)) results = pool.map(fc.calculation, - [[dem, header, infra, release_pixel, alpha, exp, flux_threshold, max_z] + [[dem, header, infra, forest, release_pixel, alpha, exp, flux_threshold, max_z] for release_pixel in release_list]) pool.close() pool.join() @@ -508,7 +548,7 @@ def main(args, kwargs): pool = mp.Pool(mp.cpu_count()) # results = pool.map(gc.calculation, iterable) results = pool.map(fc.calculation_effect, - [[dem, header, release_pixel, alpha, exp, flux_threshold, max_z] for + [[dem, header, forest, release_pixel, alpha, exp, flux_threshold, max_z] for release_pixel in release_list]) pool.close() pool.join() @@ -578,6 +618,7 @@ def main(args, kwargs): if __name__ == '__main__': #mp.set_start_method('spawn') # used in Windows argv = sys.argv[1:] + #argv = ["28", "8", "../runs/open_slope/", "../runs/open_slope/parabola.asc", "../runs/open_slope/release.tif", "forest=../runs/open_slope/forest_deep.tif", "flux=0.0003" , "max_z=270"] #"forest=../runs/open_slope/forest.tif", #argv = ['--gui'] if len(argv) < 1: print("Too few input arguments!!!") @@ -591,4 +632,5 @@ def main(args, kwargs): main(args, kwargs) # example dam: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif -# example dam: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif infra=./examples/dam/infra.tif flux=0.003 max_z=270 +# example dam/infra: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif infra=./examples/dam/infra.tif flux=0.003 max_z=270 +# example dam/forest: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif forest=./examples/dam/infra.tif flux=0.003 max_z=270 diff --git a/main_edit.py b/main_edit.py new file mode 100644 index 0000000..a1ce79f --- /dev/null +++ b/main_edit.py @@ -0,0 +1,636 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon May 7 14:23:00 2018 + + Copyright (C) <2020> + Michael.Neuhauser@bfw.gv.at + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +""" +# import standard libraries +import os +import sys +import psutil +import numpy as np +from datetime import datetime +from multiprocessing import cpu_count +import multiprocessing as mp +import logging +from xml.etree import ElementTree as ET + +# Flow-Py Libraries +import raster_io as io +import Simulation as Sim +import flow_core as fc + +# Libraries for GUI, PyQt5 +from PyQt5.QtCore import pyqtSlot, QCoreApplication +from PyQt5.QtWidgets import QFileDialog, QMessageBox, QMainWindow, QApplication + +from Flow_GUI import Ui_MainWindow + + +class Flow_Py_EXEC(): + + def __init__(self): + + app = QApplication(sys.argv) + MainWindow = QMainWindow() + self.ui = Ui_MainWindow() + self.ui.setupUi(MainWindow) + + # self.showMaximized() + #self.ui.setWindowTitle("Flow-Py") + #self.ui.setWindowIcon(QIcon('logo.jpg')) + + self.directory = os.getcwd() + + self.ui.alpha_Edit.setText('25') + self.ui.exp_Edit.setText('8') + self.ui.flux_Edit.setText('0.003') + self.ui.z_Edit.setText('8848') + + self.ui.wDir_Button.clicked.connect(self.open_wDir) + self.ui.DEM_Button.clicked.connect(self.open_dhm) + self.ui.Release_Button.clicked.connect(self.open_release) + self.ui.infra_Button.clicked.connect(self.open_infra) + self.ui.forest_Button.clicked.connect(self.open_forest) + #self.ui.process_Box.currentIndexChanged.connect(self.processChanged) + self.ui.calc_Button.clicked.connect(self.calculation) + self.ui.actionSave.triggered.connect(self.save) + self.ui.actionLoad.triggered.connect(self.load) + self.ui.actionQuit.triggered.connect(self.quit) + + self.calc_class = None + self.prot_for_bool = False + self.threads_calc = 0 + self.progress_value = 0 + self.cpu_count = 1 + self.thread_list = [] + self.start_list = [] + self.end_list = [] + for i in range(self.cpu_count): + self.thread_list.append(0) + self.start_list.append(0) + self.end_list.append(0) + + # show the constructed window + MainWindow.show() + sys.exit(app.exec_()) + + def set_gui_bool(self, bool): + self.ui.calc_Button.setEnabled(bool) + self.ui.wDir_lineEdit.setEnabled(bool) + self.ui.DEM_lineEdit.setEnabled(bool) + self.ui.release_lineEdit.setEnabled(bool) + self.ui.infra_lineEdit.setEnabled(bool) + self.ui.forest_lineEdit.setEnabled(bool) + + def save(self): + """Save the input paths""" + name = QFileDialog.getSaveFileName(None, 'Save File', + ".xml")[0] + if len(name) != 0: + + root = ET.Element('root') + wdir = ET.SubElement(root, 'wDir') + dhm = ET.SubElement(root, 'DHM') + release = ET.SubElement(root, 'Release') + infra = ET.SubElement(root, 'Infrastructure') + forest = ET.SubElement(root, 'Forest') + + wdir.set('Directory', 'Working') + dhm.set('Directory', 'DHM') + release.set('Directory', 'Release') + infra.set('Directory', 'Infrastructure') + forest.set('Directory', 'Forest') + + wdir.text = self.ui.wDir_lineEdit.text() + dhm.text = self.ui.DEM_lineEdit.text() + release.text = self.ui.release_lineEdit.text() + infra.text = self.ui.infra_lineEdit.text() + forest.text = self.ui.forest_lineEdit.text() + + tree = ET.ElementTree(root) + tree.write(name) + + def load(self): + xml_file = QFileDialog.getOpenFileNames(None, 'Open xml', + self.directory, + "xml (*.xml);;All Files (*.*)")[0] + + if len(xml_file) != 0: + tree = ET.parse(xml_file[0]) + root = tree.getroot() + + try: + self.ui.wDir_lineEdit.setText(root[0].text) + self.directory = root[0].text + except: + print("No Working Directory Path in File!") + + try: + self.ui.DEM_lineEdit.setText(root[1].text) + except: + print("No DEM Path in File!") + + try: + self.ui.release_lineEdit.setText(root[2].text) + except: + print("No Release Path in File!") + + try: + self.ui.infra_lineEdit.setText(root[3].text) + except: + print("No Infrastructure Path in File!") + + try: + self.ui.forest_lineEdit.setText(root[4].text) + except: + print("No Forest Path in File!") + + def quit(self): + QCoreApplication.quit() + + def open_wDir(self): + """Open the Working Directory, where results are stored""" + directory = QFileDialog.getExistingDirectory(None, 'Open Working Directory', + self.directory, + QFileDialog.ShowDirsOnly) + if len(directory) != 0: + self.directory = directory + self.ui.wDir_lineEdit.setText(self.directory) + + def open_dhm(self): + """Open digital elevation model""" + dem_file = QFileDialog.getOpenFileNames(None, 'Open DEM', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(dem_file[0]) != 0: + dem = dem_file[0] + self.ui.DEM_lineEdit.setText(dem[0]) + + def open_release(self): + """Open release layer""" + release_file = QFileDialog.getOpenFileNames(None, 'Open Release', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(release_file[0]) != 0: + release = release_file[0] + self.ui.release_lineEdit.setText(release[0]) + + def open_infra(self): + """Open infrastructure layer""" + infra_file = QFileDialog.getOpenFileNames(None, 'Open Infrastructure Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + if len(infra_file[0]) != 0: + infra = infra_file[0] + self.ui.infra_lineEdit.setText(infra[0]) + + def open_forest(self): + """Open forest layer""" + forest_file = QFileDialog.getOpenFileNames(None, 'Open Forest Layer', + self.directory, + "ascii (*.asc);;tif (*.tif);;All Files (*.*)") + forest = forest_file[0] + self.ui.forest_lineEdit.setText(forest[0]) + + def update_progressBar(self, float, thread, start, end): + self.thread_list[thread] = float + self.start_list[thread] = start + self.end_list[thread] = end + + self.progress_value = sum(self.thread_list) / len(self.thread_list) + self.progressBar.setValue(self.progress_value) + for i in range(len(self.thread_list)): + sys.stdout.write( + "Thread {}: Startcell {} of {} = {}%"'\n'.format(i + 1, self.start_list[i], self.end_list[i], + self.thread_list[i])) + sys.stdout.flush() + for i in range(len(self.thread_list)): + sys.stdout.write('\x1b[1A' + '\x1b[2K') # Go 1 line up and erase that line + + @staticmethod + def showdialog(path): + msg = QMessageBox() + msg.setIcon(QMessageBox.Critical) + msg.setText("No " + path + " set") + msg.setWindowTitle("Error") + msg.setStandardButtons(QMessageBox.Ok) + msg.exec_() + + def calculation(self): + self.start = datetime.now().replace(microsecond=0) + self.calc_bool = False + + # Check if input is ok + if self.ui.wDir_lineEdit.text() == '': + self.showdialog('Working Directory') + return + if self.ui.DEM_lineEdit.text() == '': + self.showdialog('DEM Layer') + return + if self.ui.release_lineEdit.text() == '': + self.showdialog('Release Layer') + return + # Disable all input line Edits and Buttons + self.set_gui_bool(False) + + # Create result directory + time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + try: + os.makedirs(self.ui.wDir_lineEdit.text() + '/result/') + self.res_dir = ('/result/') + except FileExistsError: + self.res_dir = ('/result/') + + # Setup logger + + for handler in logging.root.handlers[:]: + logging.root.removeHandler(handler) + + logging.basicConfig(level=logging.INFO, + format='%(asctime)s %(levelname)-8s %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + filename=(self.directory + self.res_dir + 'log_{}.txt').format(time_string), + filemode='w') + + # Start of Calculation + logging.info('Start Calculation') + # Read in raster files + try: + dem, header = io.read_raster(self.ui.DEM_lineEdit.text()) + logging.info('DEM File: {}'.format(self.ui.DEM_lineEdit.text())) + except FileNotFoundError: + print("Wrong filepath or filename") + self.set_gui_bool(True) + return + + try: + release, release_header = io.read_raster(self.ui.release_lineEdit.text()) + logging.info('Release File: {}'.format(self.ui.release_lineEdit.text())) + except FileNotFoundError: + print("Wrong filepath or filename") + self.set_gui_bool(True) + return + + # Check if Layers have same size!!! + if header['ncols'] == release_header['ncols'] and header['nrows'] == release_header['nrows']: + print("DEM and Release Layer ok!") + else: + print("Error: Release Layer doesn't match DEM!") + self.set_gui_bool(True) + return + + try: + infra, infra_header = io.read_raster(self.ui.infra_lineEdit.text()) + if header['ncols'] == infra_header['ncols'] and header['nrows'] == infra_header['nrows']: + print("Infra Layer ok!") + self.calc_bool = True + logging.info('Infrastructure File: {}'.format(self.ui.infra_lineEdit.text())) + else: + print("Error: Infra Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(self.ui.forest_lineEdit.text()) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + self.prot_for_bool = True + logging.info('Forest File: {}'.format(self.forest_lineEdit.text())) + else: + print("Error: Forest Layer doesn't match DEM!") + self.set_gui_bool(True) + return + except: + forest = np.zeros_like(dem) + + logging.info('Files read in') + + alpha = self.ui.alpha_Edit.text() + exp = self.ui.exp_Edit.text() + flux_threshold = self.ui.flux_Edit.text() + max_z = self.ui.z_Edit.text() + + logging.info('Alpha Angle: {}'.format(alpha)) + logging.info('Exponent: {}'.format(exp)) + logging.info('Flux Threshold: {}'.format(flux_threshold)) + logging.info('Max Z_delta: {}'.format(max_z)) + logging.info + + self.z_delta = np.zeros_like(dem) + self.flux = np.zeros_like(dem) + self.cell_counts = np.zeros_like(dem) + self.z_delta_sum = np.zeros_like(dem) + self.backcalc = np.zeros_like(dem) + self.fp_ta = np.zeros_like(dem) + self.sl_ta = np.zeros_like(dem) + + # Calculation + self.calc_class = Sim.Simulation(dem, header, release, release_header, infra, forest, self.calc_bool, alpha, exp, flux_threshold, max_z) + self.calc_class.value_changed.connect(self.update_progressBar) + self.calc_class.finished.connect(self.thread_finished) + logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) + self.calc_class.start() + + def thread_finished(self, z_delta, flux, count_array, z_delta_sum, backcalc, fp_ta, sl_ta): + logging.info('Calculation finished, getting results.') + for i in range(len(z_delta)): + self.z_delta = np.maximum(self.z_delta, z_delta[i]) + self.flux = np.maximum(self.flux, flux[i]) + self.cell_counts += count_array[i] + self.z_delta_sum += z_delta_sum[i] + self.backcalc = np.maximum(self.backcalc, backcalc[i]) + self.fp_ta = np.maximum(self.fp_ta, fp_ta[i]) + self.sl_ta = np.maximum(self.sl_ta, sl_ta[i]) + self.output() + + def output(self): + # time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + logging.info('Writing Output Files') + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "flux{}".format(self.ui.outputBox.currentText()), + self.flux) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "z_delta{}".format(self.ui.outputBox.currentText()), + self.z_delta) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "FP_travel_angle{}".format(self.ui.outputBox.currentText()), + self.fp_ta) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "SL_travel_angle{}".format(self.ui.outputBox.currentText()), + self.sl_ta) + if not self.calc_bool: + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "cell_counts{}".format(self.ui.outputBox.currentText()), + self.cell_counts) + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "z_delta_sum{}".format(self.ui.outputBox.currentText()), + self.z_delta_sum) + if self.calc_bool: + io.output_raster(self.ui.DEM_lineEdit.text(), + self.directory + self.res_dir + "backcalculation{}".format(self.ui.outputBox.currentText()), + self.backcalc) + + print("Calculation finished") + end = datetime.now().replace(microsecond=0) + logging.info('Calculation needed: ' + str(end - self.start) + ' seconds') + + # Handle GUI + #self.ui.progressBar.setValue(100) + self.set_gui_bool(True) + + +def main(args, kwargs): + + + alpha = args[0] + exp = args[1] + directory = args[2] + dem_path = args[3] + release_path = args[4] + if 'infra' in kwargs: + infra_path = kwargs.get('infra') + #print(infra_path) + else: + infra_path = None + + if 'forest' in kwargs: + forest_path = kwargs.get('forest') + #print(forest_path) + else: + forest_path = None + + if 'flux' in kwargs: + flux_threshold = float(kwargs.get('flux')) + #print(flux_threshold) + else: + flux_threshold = 3 * 10 ** -4 + + if 'max_z' in kwargs: + max_z = kwargs.get('max_z') + #print(max_z) + else: + max_z = 8848 + # Recomendet values: + # Avalanche = 270 + # Rockfall = 50 + # Soil Slide = 12 + + print("Starting...") + print("...") + + start = datetime.now().replace(microsecond=0) + calc_bool = False + # Create result directory + time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + try: + os.makedirs(directory + '/result/') + res_dir = ('/result/') + except FileExistsError: + res_dir = ('/result/') + + # Setup logger + for handler in logging.root.handlers[:]: + logging.root.removeHandler(handler) + + logging.basicConfig(level=logging.INFO, + format='%(asctime)s %(levelname)-8s %(message)s', + datefmt='%Y-%m-%d %H:%M:%S', + filename=(directory + res_dir + 'log_{}.txt').format(time_string), + filemode='w') + + # Start of Calculation + logging.info('Start Calculation') + logging.info('Alpha Angle: {}'.format(alpha)) + logging.info('Exponent: {}'.format(exp)) + logging.info('Flux Threshold: {}'.format(flux_threshold)) + logging.info('Max Z_delta: {}'.format(max_z)) + logging.info + # Read in raster files + try: + dem, header = io.read_raster(dem_path) + logging.info('DEM File: {}'.format(dem_path)) + except FileNotFoundError: + print("DEM: Wrong filepath or filename") + return + + try: + release, release_header = io.read_raster(release_path) + logging.info('Release File: {}'.format(release_path)) + except FileNotFoundError: + print("Wrong filepath or filename") + return + + # Check if Layers have same size!!! + if header['ncols'] == release_header['ncols'] and header['nrows'] == release_header['nrows']: + print("DEM and Release Layer ok!") + else: + print("Error: Release Layer doesn't match DEM!") + return + + try: + infra, infra_header = io.read_raster(infra_path) + if header['ncols'] == infra_header['ncols'] and header['nrows'] == infra_header['nrows']: + print("Infra Layer ok!") + calc_bool = True + logging.info('Infrastructure File: {}'.format(infra_path)) + else: + print("Error: Infra Layer doesn't match DEM!") + return + except: + infra = np.zeros_like(dem) + + try: + forest, forest_header = io.read_raster(forest_path) + if header['ncols'] == forest_header['ncols'] and header['nrows'] == forest_header['nrows']: + print("Forest Layer ok!") + prot_for_bool = True + logging.info('Forest File: {}'.format(forest_path)) + else: + print("Error: Forest Layer doesn't match DEM!") + return + except: + forest = np.zeros_like(dem) + + logging.info('Files read in') + + z_delta = np.zeros_like(dem) + flux = np.zeros_like(dem) + cell_counts = np.zeros_like(dem) + z_delta_sum = np.zeros_like(dem) + backcalc = np.zeros_like(dem) + fp_ta = np.zeros_like(dem) + sl_ta = np.zeros_like(dem) + + avaiable_memory = psutil.virtual_memory()[1] + needed_memory = sys.getsizeof(dem) + + max_number_procces = int(avaiable_memory / (needed_memory * 10)) + + print( + "There are {} Bytes of Memory avaiable and {} Bytes needed per process. \nMax. Nr. of Processes = {}".format( + avaiable_memory, needed_memory*10, max_number_procces)) + + # Calculation + logging.info('Multiprocessing starts, used cores: {}'.format(cpu_count())) + + if calc_bool: + release_list = fc.split_release(release, release_header, min(mp.cpu_count() * 2, max_number_procces)) + + print("{} Processes started.".format(len(release_list))) + pool = mp.Pool(len(release_list)) + results = pool.map(fc.calculation, + [[dem, header, infra, forest, release_pixel, alpha, exp, flux_threshold, max_z] + for release_pixel in release_list]) + pool.close() + pool.join() + else: + release_list = fc.split_release(release, release_header, min(mp.cpu_count() * 4, max_number_procces)) + + print("{} Processes started.".format(len(release_list))) + pool = mp.Pool(mp.cpu_count()) + # results = pool.map(gc.calculation, iterable) + results = pool.map(fc.calculation_effect, + [[dem, header, forest, release_pixel, alpha, exp, flux_threshold, max_z] for + release_pixel in release_list]) + pool.close() + pool.join() + + z_delta_list = [] + flux_list = [] + cc_list = [] + z_delta_sum_list = [] + backcalc_list = [] + fp_ta_list = [] + sl_ta_list = [] + for i in range(len(results)): + res = results[i] + res = list(res) + z_delta_list.append(res[0]) + flux_list.append(res[1]) + cc_list.append(res[2]) + z_delta_sum_list.append(res[3]) + backcalc_list.append(res[4]) + fp_ta_list.append(res[5]) + sl_ta_list.append(res[6]) + + logging.info('Calculation finished, getting results.') + for i in range(len(z_delta_list)): + z_delta = np.maximum(z_delta, z_delta_list[i]) + flux = np.maximum(flux, flux_list[i]) + cell_counts += cc_list[i] + z_delta_sum += z_delta_sum_list[i] + backcalc = np.maximum(backcalc, backcalc_list[i]) + fp_ta = np.maximum(fp_ta, fp_ta_list[i]) + sl_ta = np.maximum(sl_ta, sl_ta_list[i]) + + + # time_string = datetime.now().strftime("%Y%m%d_%H%M%S") + logging.info('Writing Output Files') + output_format = '.tif' + io.output_raster(dem_path, + directory + res_dir + "flux{}".format(output_format), + flux) + io.output_raster(dem_path, + directory + res_dir + "z_delta{}".format(output_format), + z_delta) + io.output_raster(dem_path, + directory + res_dir + "FP_travel_angle{}".format(output_format), + fp_ta) + io.output_raster(dem_path, + directory + res_dir + "SL_travel_angle{}".format(output_format), + sl_ta) + if not calc_bool: # if no infra + io.output_raster(dem_path, + directory + res_dir + "cell_counts{}".format(output_format), + cell_counts) + io.output_raster(dem_path, + directory + res_dir + "z_delta_sum{}".format( output_format), + z_delta_sum) + if calc_bool: # if infra + io.output_raster(dem_path, + directory + res_dir + "backcalculation{}".format(output_format), + backcalc) + + print("Calculation finished") + print("...") + end = datetime.now().replace(microsecond=0) + logging.info('Calculation needed: ' + str(end - start) + ' seconds') + + +if __name__ == '__main__': + #mp.set_start_method('spawn') # used in Windows + argv = sys.argv[1:] + #argv = ["28", "8", "../runs/open_slope/", "../runs/open_slope/parabola.asc", "../runs/open_slope/release.tif", "forest=../runs/open_slope/forest_deep.tif", "flux=0.0003" , "max_z=270"] #"forest=../runs/open_slope/forest.tif", + #argv = ['--gui'] + if len(argv) < 1: + print("Too few input arguments!!!") + sys.exit(1) + if len(argv) == 1 and argv[0] == '--gui': + Flow_Py_EXEC() + else: + args=[arg for arg in argv if arg.find('=')<0] + kwargs={kw[0]:kw[1] for kw in [ar.split('=') for ar in argv if ar.find('=')>0]} + + main(args, kwargs) + +# example dam: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif +# example dam/infra: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif infra=./examples/dam/infra.tif flux=0.003 max_z=270 +# example dam/forest: python3 main.py 25 8 ./examples/dam/ ./examples/dam/dam_010m_standard_cr100_sw250_f2500.20.6_n0.asc ./examples/dam/release_dam.tif forest=./examples/dam/infra.tif flux=0.003 max_z=270