From 62aeb9a04a8fcc236ed4db0f4c8f204a7789ddbf Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:00:19 +0200 Subject: [PATCH 01/28] Up gdal version (whl) --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 0d7f3f55..a02946d3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -78,7 +78,7 @@ all = [ "setuptools>=61.0.0", ] build = [ - "GDAL @ https://github.com/cgohlke/geospatial-wheels/releases/download/v2023.4.22/GDAL-3.6.4-cp311-cp311-win_amd64.whl", + "GDAL @ https://github.com/cgohlke/geospatial-wheels/releases/download/v2023.7.16/GDAL-3.7.1-cp311-cp311-win_amd64.whl", "pyinstaller", ] dev = [ From e4923e02528fc720f81527b75d32f54a282aefdf Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:01:53 +0200 Subject: [PATCH 02/28] Added gdal cache --- src/delft_fiat/cfg.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/delft_fiat/cfg.py b/src/delft_fiat/cfg.py index d5fde06f..d45713cd 100644 --- a/src/delft_fiat/cfg.py +++ b/src/delft_fiat/cfg.py @@ -9,6 +9,7 @@ import os import tomli +from osgeo import gdal class ConfigReader(dict): @@ -38,6 +39,13 @@ def __init__( # Create the hidden temporary folder self._create_temp_dir() + # Set the cache size per GDAL object + _cache_size = self.get("global.gdal_cache") + if _cache_size is not None: + gdal.SetCacheMax(_cache_size * 1024**2) + else: + gdal.SetCacheMax(50 * 1024**2) + # Do some checking concerning the file paths in the settings file for key, item in self.items(): if key.endswith(("file", "csv")) or key.rsplit(".", 1)[1].startswith( From 5c27bb3660b2de10eeb3b17271e190fa0a0de0b9 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:02:32 +0200 Subject: [PATCH 03/28] Message update --- src/delft_fiat/cli/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/delft_fiat/cli/main.py b/src/delft_fiat/cli/main.py index ae814d06..d9b20395 100644 --- a/src/delft_fiat/cli/main.py +++ b/src/delft_fiat/cli/main.py @@ -88,7 +88,7 @@ def run( dst=cfg.get("output.path"), ) sys.stdout.write(fiat_start_str) - logger.info(f"Delft-Fiat version: {__version__}") + logger.info(f"Delft-FIAT version: {__version__}") obj = FIAT(cfg) obj.run() From 62c405be56ac9449f5ad74f30c085cdd79eaabce Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:03:08 +0200 Subject: [PATCH 04/28] Added Grid structure for bands --- src/delft_fiat/gis/grid.py | 20 ++++ src/delft_fiat/io.py | 208 ++++++++++++++++++++++++++++++++----- 2 files changed, 204 insertions(+), 24 deletions(-) diff --git a/src/delft_fiat/gis/grid.py b/src/delft_fiat/gis/grid.py index 904a2d1d..6a92e8c1 100644 --- a/src/delft_fiat/gis/grid.py +++ b/src/delft_fiat/gis/grid.py @@ -6,6 +6,26 @@ from pathlib import Path +def clip( + band: gdal.Band, + gtf: tuple, + idx: tuple, +): + """_summary_ + + Parameters + ---------- + band : gdal.Band + _description_ + gtf : tuple + _description_ + idx : tuple + _description_ + """ + + pass + + def reproject( gs: "GridSource", crs: str, diff --git a/src/delft_fiat/io.py b/src/delft_fiat/io.py index 1562401c..3a6a8fa9 100644 --- a/src/delft_fiat/io.py +++ b/src/delft_fiat/io.py @@ -1,10 +1,10 @@ -from typing import Any from delft_fiat.error import DriverNotFoundError from delft_fiat.util import ( GEOM_DRIVER_MAP, GRID_DRIVER_MAP, Path, deter_type, + replace_empty, _dtypes_from_string, _dtypes_reversed, _pat, @@ -16,13 +16,16 @@ import atexit import gc import os +import time import weakref from abc import ABCMeta, abstractmethod from io import BufferedReader, BufferedWriter, FileIO, TextIOWrapper +from itertools import product from math import nan, floor, log10 -from numpy import arange, array, column_stack, interp, ndarray +from numpy import arange, array, column_stack, interp, ndarray, where from osgeo import gdal, ogr from osgeo import osr +from typing import Any _IOS = weakref.WeakValueDictionary() _IOS_COUNT = 1 @@ -508,6 +511,122 @@ def read_exp(self): ## Structs +class Grid(_BaseStruct): + def __init__( + self, + band: gdal.Band, + chunk: tuple = None, + ): + """_summary_""" + + self.src = band + self._x = band.XSize + self._y = band.YSize + self._l = 0 + self._u = 0 + self.nodata = band.GetNoDataValue() + self.dtype = band.DataType + self.dtype_name = gdal.GetDataTypeName(self.dtype) + self.dtype_size = gdal.GetDataTypeSize(self.dtype) + + self._last_chunk = None + + if chunk is None: + self._chunk = self.shape + elif len(chunk) == 2: + self._chunk = chunk + else: + ValueError("") + + self._recreate_windows() + + def __del__(self): + self.flush() + self.src = None + + def __iter__(self): + self.flush() + # self._recreate_windows() + self._reset_chunking() + return self + + def __next__(self): + if self._u > self._y: + self.flush() + raise StopIteration + + w = min(self._chunk[0], self._x - self._l) + h = min(self._chunk[1], self._y - self._u) + + window = ( + self._l, + self._u, + w, + h, + ) + chunk = self[window] + + self._l += self._chunk[0] + if self._l > self._x: + self._l = 0 + self._u += self._chunk[1] + + return window, chunk + + # for l, u in self._windows: + # s = time.time() + + # if self._last_chunk is not None: + # self._last_chunk = None + + # w = min(self._chunk[0], self._x - l) + # h = min(self._chunk[1], self._y - u) + + # chunk = self[l, u, w, h] + # self._last_chunk = chunk + + # e = time.time() - s + # print(f"{e} seconds! (iter)") + + # return (l, u, w, h), chunk + + def __getitem__( + self, + window: tuple, + ): + chunk = self.src.ReadAsArray(*window) + return chunk + + def _recreate_windows(self): + self._windows = product( + range(0, self._x, self._chunk[0]), + range(0, self._y, self._chunk[1]), + ) + + def _reset_chunking(self): + self._l = 0 + self._u = 0 + + def flush(self): + self.src.FlushCache() + + @property + def chunking(self): + return self._chunk + + @property + def shape(self): + return self._x, self._y + + def set_chunking( + self, + chunk: tuple, + ): + """_summary_""" + + self._chunk = chunk + + class GeomSource(_BaseIO, _BaseStruct): _type_map = { "int": ogr.OFTInteger64, @@ -784,6 +903,7 @@ class GridSource(_BaseIO, _BaseStruct): def __new__( cls, file: str, + chunk: tuple = None, subset: str = None, var_as_band: bool = False, mode: str = "r", @@ -791,13 +911,14 @@ def __new__( """_summary_""" obj = object.__new__(cls) - obj.__init__(file, subset, var_as_band, mode) + obj.__init__(file, chunk, subset, var_as_band, mode) return obj def __init__( self, file: str, + chunk: tuple = None, subset: str = None, var_as_band: bool = False, mode: str = "r", @@ -808,11 +929,14 @@ def __init__( ---------- file : str _description_ - - Returns - ------- - object - _description_ + chunk : tuple, optional + _description_, by default None + subset : str, optional + _description_, by default None + var_as_band : bool, optional + _description_, by default False + mode : str, optional + _description_, by default "r" Raises ------ @@ -849,6 +973,8 @@ def __init__( self._driver = gdal.GetDriverByName(driver) self.src = None + self._chunk = None + self._dtype = None self.subset_dict = None self.count = 0 self._cur_index = 1 @@ -857,6 +983,13 @@ def __init__( self.src = gdal.OpenEx(str(self._path), open_options=_open_options) self.count = self.src.RasterCount + if chunk is None: + self._chunk = self.shape + elif len(chunk) == 2: + self._chunk = chunk + else: + ValueError("") + if self.count == 0: self.subset_dict = _read_gridsrouce_layers( self.src, @@ -868,7 +1001,7 @@ def __iter__(self): def __next__(self): if self._cur_index < self.count + 1: - r = self.src.GetRasterBand(self._cur_index) + r = self[self._cur_index] self._cur_index += 1 return r else: @@ -878,11 +1011,15 @@ def __getitem__( self, oid: int, ): - return self.src.GetRasterBand(oid) + return Grid( + self.src.GetRasterBand(oid), + chunk=self._chunk, + ) def __reduce__(self): return self.__class__, ( self.path, + self._chunk, self.subset, self._var_as_band, self._mode_str, @@ -907,38 +1044,59 @@ def reopen(self): return GridSource.__new__( GridSource, self.path, + self._chunk, self.subset, self._var_as_band, ) - @_BaseIO._check_mode + @property @_BaseIO._check_state - def create_band( - self, - ): - """_summary_""" + def dtype(self): + if not self._dtype: + _b = self[1] + self._dtype = _b.dtype + del _b + return self._dtype - self.src.AddBand() - self.count += 1 + @property + @_BaseIO._check_state + def shape(self): + return ( + self.src.RasterXSize, + self.src.RasterYSize, + ) @_BaseIO._check_mode @_BaseIO._check_state - def create_source( + def create( self, shape: tuple, nb: int, type: int, - srs: osr.SpatialReference, + options: list = [], ): """_summary_""" self.src = self._driver.Create( - str(self.path), shape[0], shape[1], nb, GridSource._type_map[type] + str(self.path), + *shape, + nb, + type, + options=options, ) - self.src.SetSpatialRef(srs) self.count = nb + @_BaseIO._check_mode + @_BaseIO._check_state + def create_band( + self, + ): + """_summary_""" + + self.src.AddBand() + self.count += 1 + @_BaseIO._check_state def deter_band_names( self, @@ -959,8 +1117,8 @@ def deter_band_names( def get_band_name(self, n: int): """_summary_""" - _desc = self[n].GetDescription() - _meta = self[n].GetMetadata() + _desc = self[n].src.GetDescription() + _meta = self[n].src.GetMetadata() if _desc: return _desc @@ -1211,7 +1369,7 @@ def _build_from_stream( cols.remove(index_col) for c in cols: - _f.append([dtypes[c](item.decode()) for item in _d[c::ncol]]) + _f.append([dtypes[c](item) for item in replace_empty(_d[c::ncol])]) self.data = column_stack((*_f,)) @@ -1558,6 +1716,7 @@ def open_geom( def open_grid( file: str, + chunk: tuple = None, subset: str = None, var_as_band: bool = False, mode: str = "r", @@ -1566,6 +1725,7 @@ def open_grid( return GridSource( file, + chunk, subset, var_as_band, mode, From 72ec505008fcf3e78097b6cc382517db24153ef2 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:03:33 +0200 Subject: [PATCH 05/28] private to public variables --- src/delft_fiat/models/base.py | 16 ++++++++-------- src/delft_fiat/models/geom.py | 34 +++++++++++++++++----------------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/delft_fiat/models/base.py b/src/delft_fiat/models/base.py index 0ea6dfbe..e1a7443c 100644 --- a/src/delft_fiat/models/base.py +++ b/src/delft_fiat/models/base.py @@ -31,11 +31,11 @@ def __init__( # Declarations self.srs = None - self._exposure_data = None - self._exposure_geoms = None - self._exposure_grid = None - self._hazard_grid = None - self._vulnerability_data = None + self.exposure_data = None + self.exposure_geoms = None + self.exposure_grid = None + self.hazard_grid = None + self.vulnerability_data = None self._vul_step_size = 0.01 self._rounding = 2 self._cfg["vulnerability.round"] = self._rounding @@ -66,7 +66,7 @@ def _read_hazard_grid(self): path = self._cfg.get("hazard.file") logger.info(f"Reading hazard data ('{path.name}')") - kw = self._cfg.generate_kwargs("hazard.multiband") + kw = self._cfg.generate_kwargs("hazard.settings") data = open_grid(path, **kw) ## checks logger.info("Executing hazard checks...") @@ -123,7 +123,7 @@ def _read_hazard_grid(self): self._cfg["hazard.band_names"] = ns # When all is done, add it - self._hazard_grid = data + self.hazard_grid = data def _read_vulnerability_data(self): path = self._cfg.get("vulnerability.file") @@ -148,7 +148,7 @@ def _read_vulnerability_data(self): ) data.upscale(self._vul_step_size, inplace=True) # When all is done, add it - self._vulnerability_data = data + self.vulnerability_data = data def _set_model_srs(self): """_summary_""" diff --git a/src/delft_fiat/models/geom.py b/src/delft_fiat/models/geom.py index a786fc8f..b5029455 100644 --- a/src/delft_fiat/models/geom.py +++ b/src/delft_fiat/models/geom.py @@ -71,7 +71,7 @@ def _read_exposure_data(self): self._cfg["output.new_columns"] = cols ## When all is done, add it - self._exposure_data = data + self.exposure_data = data def _read_exposure_geoms(self): """_summary_""" @@ -104,15 +104,15 @@ def _read_exposure_geoms(self): # Add to the dict _d[file.rsplit(".", 1)[1]] = data # When all is done, add it - self._exposure_geoms = _d + self.exposure_geoms = _d - def _patch_up( + def resolve( self, ): """_summary_""" - _exp = self._exposure_data - _gm = self._exposure_geoms + _exp = self.exposure_data + _gm = self.exposure_geoms _risk = self._cfg.get("hazard.risk") _rp_coef = self._cfg.get("hazard.rp_coefficients") _new_cols = self._cfg["output.new_columns"] @@ -199,23 +199,23 @@ def run( _nms = self._cfg.get("hazard.band_names") logger.info("Starting the calculations") - if self._hazard_grid.count > 1: - pcount = min(os.cpu_count(), self._hazard_grid.count) + if self.hazard_grid.count > 1: + pcount = min(os.cpu_count(), self.hazard_grid.count) futures = [] with ProcessPoolExecutor(max_workers=pcount) as Pool: _s = time.time() - for idx in range(self._hazard_grid.count): + for idx in range(self.hazard_grid.count): logger.info( f"Submitting a job for the calculations in regards to band: '{_nms[idx]}'" ) fs = Pool.submit( geom_worker, self._cfg, - self._hazard_grid, + self.hazard_grid, idx + 1, - self._vulnerability_data, - self._exposure_data, - self._exposure_geoms, + self.vulnerability_data, + self.exposure_data, + self.exposure_geoms, ) futures.append(fs) logger.info("Busy...") @@ -229,11 +229,11 @@ def run( target=geom_worker, args=( self._cfg, - self._hazard_grid, + self.hazard_grid, 1, - self._vulnerability_data, - self._exposure_data, - self._exposure_geoms, + self.vulnerability_data, + self.exposure_data, + self.exposure_geoms, ), ) p.start() @@ -243,7 +243,7 @@ def run( logger.info(f"Calculations time: {round(_e, 2)} seconds") logger.info("Producing model output from temporary files") - self._patch_up() + self.resolve() logger.info(f"Output generated in: '{self._cfg['output.path']}'") if not self._keep_temp: From 52ec58e8b47785afcced30e45f7eb30c3953305b Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:03:51 +0200 Subject: [PATCH 06/28] Raster progress --- src/delft_fiat/models/grid.py | 36 ++++++++++++++++------ src/delft_fiat/models/util.py | 57 +++++++++++++++++++++++++++++++++-- 2 files changed, 81 insertions(+), 12 deletions(-) diff --git a/src/delft_fiat/models/grid.py b/src/delft_fiat/models/grid.py index a6aa310f..0cd4aca6 100644 --- a/src/delft_fiat/models/grid.py +++ b/src/delft_fiat/models/grid.py @@ -1,7 +1,9 @@ +from delft_fiat.check import * from delft_fiat.io import BufferTextHandler, open_grid from delft_fiat.log import spawn_logger from delft_fiat.models.base import BaseModel -from delft_fiat.util import _pat, replace_empty +from delft_fiat.models.util import grid_worker_exact +from delft_fiat.util import NEWLINE_CHAR from concurrent.futures import ProcessPoolExecutor from multiprocessing import Process, get_context @@ -9,31 +11,45 @@ logger = spawn_logger("fiat.model.grid") -def worker(haz, vul, exp): - pass - - class GridModel(BaseModel): def __init__( self, cfg: "ConfigReader", ): """_summary_""" - pass + + super().__init__(cfg) + + self._read_exposure_grid() def __del__(self): BaseModel.__del__(self) + def _clean_up(self): + pass + def _read_exposure_grid(self): """_summary_""" - path = self._cfg.get("exposure.grid.file") - data = open_grid(path) + file = self._cfg.get("exposure.grid.file") + logger.info(f"Reading exposure grid ('{file.name}')") + kw = self._cfg.generate_kwargs("exposure.grid.settings") + data = open_grid(file, **kw) ## checks + logger.info("Executing exposure data checks...") - self._exposure_grid = data + self.exposure_grid = data + + def resolve(): + pass def run(self): """_summary_""" - worker() + grid_worker_exact( + self._cfg, + self.hazard_grid, + 1, + self.vulnerability_data, + self.exposure_grid, + ) diff --git a/src/delft_fiat/models/util.py b/src/delft_fiat/models/util.py index be13dced..8f88e01f 100644 --- a/src/delft_fiat/models/util.py +++ b/src/delft_fiat/models/util.py @@ -1,9 +1,11 @@ -from delft_fiat.io import BufferTextHandler +from delft_fiat.io import BufferTextHandler, GridSource from delft_fiat.gis import geom, overlay from delft_fiat.models.calc import calc_haz from delft_fiat.util import NEWLINE_CHAR, _pat, replace_empty from math import isnan +from numpy import full, ravel, unravel_index, where +from osgeo import gdal, osr from pathlib import Path @@ -83,9 +85,60 @@ def geom_worker( writer = None -def grid_worker( +def grid_worker_exact( cfg: "ConfigReader", + haz: "GridSource", + idx: int, + vul: "Table", + exp: "GridSource", ): """_summary_""" + exp_nd = exp[1].nodata + + out_src = GridSource( + "C:\\temp\\output.nc", + mode="w", + ) + + out_src.create( + exp.shape, + 1, + exp.dtype, + options=["FORMAT=NC4", "COMPRESS=DEFLATE"], + ) + out_src.set_srs(exp.get_srs()) + out_src.set_geotransform(exp.get_geotransform()) + + write_band = out_src[1] + write_band.src.SetNoDataValue(exp_nd) + + for (_, h_ch), (_w, e_ch) in zip(haz[idx], exp[1]): + out_ch = full(e_ch.shape, exp_nd) + e_ch = ravel(e_ch) + _coords = where(e_ch != exp_nd)[0] + if len(_coords) == 0: + write_band.src.WriteArray(out_ch, *_w[:2]) + continue + + e_ch = e_ch[_coords] + h_ch = ravel(h_ch) + h_ch = h_ch[_coords] + _hcoords = where(h_ch != haz[idx].nodata)[0] + _coords = _coords[_hcoords] + e_ch = e_ch[_hcoords] + h_ch = h_ch[_hcoords] + + pass + + write_band.flush() + write_band = None + out_src = None + + pass + + +def grid_worker_loose(): + """_summary_""" + pass From cb70b1f6d1f3c3de57fe2faf2652b84d04d8d294 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 30 Aug 2023 12:04:20 +0200 Subject: [PATCH 07/28] Switched to grid struct --- src/delft_fiat/gis/overlay.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/delft_fiat/gis/overlay.py b/src/delft_fiat/gis/overlay.py index aac08edf..fd770a0e 100644 --- a/src/delft_fiat/gis/overlay.py +++ b/src/delft_fiat/gis/overlay.py @@ -1,10 +1,11 @@ from delft_fiat.gis.util import world2pixel, pixel2world +from delft_fiat.io import Grid from osgeo import gdal, ogr def clip( - band: gdal.Band, + band: Grid, srs: "osr.SpatialReference", gtf: tuple, ft: ogr.Feature, @@ -36,7 +37,7 @@ def clip( pxWidth = int(lrX - ulX) + 1 pxHeight = int(lrY - ulY) + 1 - clip = band.ReadAsArray(ulX, ulY, pxWidth, pxHeight) + clip = band[ulX, ulY, pxWidth, pxHeight] # m = mask.ReadAsArray(ulX,ulY,pxWidth,pxHeight) # pts = geom.GetGeometryRef(0) @@ -140,7 +141,7 @@ def mask( def pin( - band: gdal.Band, + band: Grid, gtf: tuple, point: tuple, ) -> "numpy.array": @@ -163,6 +164,6 @@ def pin( X, Y = world2pixel(gtf, *point) - value = band.ReadAsArray(X, Y, 1, 1) + value = band[X, Y, 1, 1] return value[0] From 14df10eafd1aafe0880e15645e3eb76bb143bd63 Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 1 Sep 2023 14:23:11 +0200 Subject: [PATCH 08/28] Bypass missing info for now... --- src/delft_fiat/models/util.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/delft_fiat/models/util.py b/src/delft_fiat/models/util.py index bdaf12d9..f9ede144 100644 --- a/src/delft_fiat/models/util.py +++ b/src/delft_fiat/models/util.py @@ -46,6 +46,8 @@ def geom_worker( row = b"" ft_info_raw = exp[ft.GetField(0)] + if ft_info_raw is None: + continue ft_info = replace_empty(_pat.split(ft_info_raw)) ft_info = [x(y) for x, y in zip(exp.dtypes, ft_info)] row += f"{ft_info[exp.index_col]}".encode() From 06de1c7f3abd0f158bbce93ba20a415693b3c08d Mon Sep 17 00:00:00 2001 From: Brendan Date: Fri, 8 Sep 2023 16:49:10 +0200 Subject: [PATCH 09/28] Working raster model for multiple exposure bands --- src/delft_fiat/io.py | 96 +++++++++++++++++++++++++++-------- src/delft_fiat/models/util.py | 63 ++++++++++++++--------- src/delft_fiat/util.py | 9 ++++ 3 files changed, 123 insertions(+), 45 deletions(-) diff --git a/src/delft_fiat/io.py b/src/delft_fiat/io.py index 190f7b3d..396271ec 100644 --- a/src/delft_fiat/io.py +++ b/src/delft_fiat/io.py @@ -1,5 +1,6 @@ from delft_fiat.error import DriverNotFoundError from delft_fiat.util import ( + DoNotCall, GEOM_DRIVER_MAP, GRID_DRIVER_MAP, Path, @@ -56,20 +57,25 @@ class _BaseIO(metaclass=ABCMeta): "w": 1, } + _closed = False + _path = None + path = None + src = None + def __init__( self, - file: str, + file: str = None, mode: str = "r", ): """_summary_""" + if file is not None: + self.path = Path(file) + self._path = Path(file) + if not mode in _BaseIO._mode_map: raise ValueError("") - self.path = Path(file) - self._path = Path(file) - - self._closed = False self._mode = _BaseIO._mode_map[mode] self._mode_str = mode @@ -511,14 +517,20 @@ def read_exp(self): ## Structs -class Grid(_BaseStruct): +class Grid( + _BaseIO, + _BaseStruct, +): def __init__( self, band: gdal.Band, chunk: tuple = None, + mode: str = "r", ): """_summary_""" + _BaseIO.__init__(self, mode=mode) + self.src = band self._x = band.XSize self._y = band.YSize @@ -540,13 +552,8 @@ def __init__( self._recreate_windows() - def __del__(self): - self.flush() - self.src = None - def __iter__(self): self.flush() - # self._recreate_windows() self._reset_chunking() return self @@ -607,8 +614,14 @@ def _reset_chunking(self): self._l = 0 self._u = 0 + def close(self): + _BaseIO.close(self) + self.src = None + gc.collect() + def flush(self): - self.src.FlushCache() + if self.src is not None: + self.src.FlushCache() @property def chunking(self): @@ -618,6 +631,14 @@ def chunking(self): def shape(self): return self._x, self._y + def get_metadata_item( + self, + entry: str, + ): + """_summary_""" + + return self.src.GetMetadataItem(entry) + def set_chunking( self, chunk: tuple, @@ -626,6 +647,37 @@ def set_chunking( self._chunk = chunk + @_BaseIO._check_mode + def write( + self, + data: array, + ): + """_summary_ + + Parameters + ---------- + data : array + _description_ + """ + + pass + + @_BaseIO._check_mode + def write_chunk( + self, + chunk: array, + upper_left: tuple | list, + ): + """_summary_ + + Parameters + ---------- + chunk : array + _description_ + """ + + self.src.WriteArray(chunk, *upper_left) + class GeomSource(_BaseIO, _BaseStruct): _type_map = { @@ -642,10 +694,6 @@ def __new__( """_summary_""" obj = object.__new__(cls) - obj.__init__( - file, - mode, - ) return obj @@ -741,14 +789,17 @@ def close(self): # return self.layer.GetFeatureCount() def flush(self): - self.src.FlushCache() + if self.src is not None: + self.src.FlushCache() def reopen(self): """_summary_""" if not self._closed: return self - return GeomSource.__new__(GeomSource, self.path) + obj = GeomSource.__new__(GeomSource, self.path) + obj.__init__(self.path) + return obj @_BaseIO._check_mode @_BaseIO._check_state @@ -911,7 +962,6 @@ def __new__( """_summary_""" obj = object.__new__(cls) - obj.__init__(file, chunk, subset, var_as_band, mode) return obj @@ -1014,6 +1064,7 @@ def __getitem__( return Grid( self.src.GetRasterBand(oid), chunk=self._chunk, + mode=self._mode_str, ) def __reduce__(self): @@ -1034,20 +1085,23 @@ def close(self): gc.collect() def flush(self): - self.src.FlushCache() + if self.src is not None: + self.src.FlushCache() def reopen(self): """_summary_""" if not self._closed: return self - return GridSource.__new__( + obj = GridSource.__new__( GridSource, self.path, self._chunk, self.subset, self._var_as_band, ) + obj.__init__(self.path, self._chunk, self.subset, self._var_as_band) + return obj @property @_BaseIO._check_state diff --git a/src/delft_fiat/models/util.py b/src/delft_fiat/models/util.py index f9ede144..15f72aac 100644 --- a/src/delft_fiat/models/util.py +++ b/src/delft_fiat/models/util.py @@ -98,6 +98,7 @@ def grid_worker_exact( """_summary_""" exp_nd = exp[1].nodata + dmf = exp[1].get_metadata_item("damage_function") out_src = GridSource( "C:\\temp\\output.nc", @@ -106,36 +107,50 @@ def grid_worker_exact( out_src.create( exp.shape, - 1, + exp.count, exp.dtype, options=["FORMAT=NC4", "COMPRESS=DEFLATE"], ) out_src.set_srs(exp.get_srs()) out_src.set_geotransform(exp.get_geotransform()) - write_band = out_src[1] - write_band.src.SetNoDataValue(exp_nd) - - for (_, h_ch), (_w, e_ch) in zip(haz[idx], exp[1]): - out_ch = full(e_ch.shape, exp_nd) - e_ch = ravel(e_ch) - _coords = where(e_ch != exp_nd)[0] - if len(_coords) == 0: - write_band.src.WriteArray(out_ch, *_w[:2]) - continue - - e_ch = e_ch[_coords] - h_ch = ravel(h_ch) - h_ch = h_ch[_coords] - _hcoords = where(h_ch != haz[idx].nodata)[0] - _coords = _coords[_hcoords] - e_ch = e_ch[_hcoords] - h_ch = h_ch[_hcoords] - - pass - - write_band.flush() - write_band = None + for _bandn in range(exp.count): + write_band = out_src[_bandn + 1] + write_band.src.SetNoDataValue(exp_nd) + + for (_, h_ch), (_w, e_ch) in zip(haz[idx], exp[_bandn + 1]): + out_ch = full(e_ch.shape, exp_nd) + e_ch = ravel(e_ch) + _coords = where(e_ch != exp_nd)[0] + if len(_coords) == 0: + write_band.src.WriteArray(out_ch, *_w[:2]) + continue + + e_ch = e_ch[_coords] + h_ch = ravel(h_ch) + h_ch = h_ch[_coords] + _hcoords = where(h_ch != haz[idx].nodata)[0] + + if len(_hcoords) == 0: + write_band.src.WriteArray(out_ch, *_w[:2]) + continue + + _coords = _coords[_hcoords] + e_ch = e_ch[_hcoords] + h_ch = h_ch[_hcoords] + h_ch.clip(min(vul.index), max(vul.index)) + + dmm = [vul[round(float(n), 2), dmf] for n in h_ch] + e_ch = e_ch * dmm + + idx2d = unravel_index(_coords, *[exp._chunk]) + out_ch[idx2d] = e_ch + + write_band.write_chunk(out_ch, _w[:2]) + + write_band.flush() + write_band = None + out_src.flush() out_src = None pass diff --git a/src/delft_fiat/util.py b/src/delft_fiat/util.py index 08b412e3..78511a30 100644 --- a/src/delft_fiat/util.py +++ b/src/delft_fiat/util.py @@ -174,6 +174,15 @@ def _text_chunk_gen( yield _nlines, sd +class DoNotCall(type): + def __call__( + self, + *args, + **kwargs, + ): + raise AttributeError("Cannot initialize directly, needs a contructor") + + def replace_empty(l: list): """_summary_""" From 519e596c965d65186c778957b1fc408df1c7b1ab Mon Sep 17 00:00:00 2001 From: Brendan Date: Mon, 18 Sep 2023 18:12:19 +0200 Subject: [PATCH 10/28] Integrator progress --- src/delft_fiat/cfg.py | 18 +- src/delft_fiat/check.py | 458 ++++++---------------------------- src/delft_fiat/main.py | 11 +- src/delft_fiat/models/geom.py | 4 +- src/delft_fiat/models/grid.py | 1 + 5 files changed, 98 insertions(+), 394 deletions(-) diff --git a/src/delft_fiat/cfg.py b/src/delft_fiat/cfg.py index d45713cd..767204d7 100644 --- a/src/delft_fiat/cfg.py +++ b/src/delft_fiat/cfg.py @@ -1,4 +1,8 @@ -from delft_fiat.check import check_config_entries +from delft_fiat.check import ( + check_config_entries, + check_config_geom, + check_config_grid, +) from delft_fiat.util import ( Path, create_hidden_folder, @@ -89,10 +93,14 @@ def get_model_type( ): """_Summary_""" - if "exposure.geom_file" in self: - return 0 - else: - return 1 + _models = [False, False] + + if check_config_geom(self): + _models[0] = True + if check_config_grid(self): + _models[1] = True + + return _models def get_path( self, diff --git a/src/delft_fiat/check.py b/src/delft_fiat/check.py index 030ebd93..9f441237 100644 --- a/src/delft_fiat/check.py +++ b/src/delft_fiat/check.py @@ -38,6 +38,55 @@ def check_config_entries( sys.exit() +def check_config_geom( + cfg: 'ConfigReader', +): + """_summary_""" + + _req_cols = [ + "exposure.geom.crs", + "exposure.geom.csv", + "exposure.geom.file1", + ] + _all_geom = [ + item for item in cfg if item.startswith("exposure.geom") + ] + if len(_all_geom) == 0: + return False + + _check = [item in _all_geom for item in _req_cols] + if not all(_check): + _missing = [item for item, b in zip(_req_cols, _check) if not b] + logger.warning(f"Info for the geometry model was found, but not all. {_missing} was/ were missing") + return False + + return True + + +def check_config_grid( + cfg: 'ConfigReader', +): + """_summary_""" + + _req_cols = [ + "exposure.grid.crs", + "exposure.grid.file", + ] + _all_grid = [ + item for item in cfg if item.startswith("exposure.grid") + ] + if len(_all_grid) == 0: + return False + + _check = [item in _all_grid for item in _req_cols] + if not all(_check): + _missing = [item for item, b in zip(_req_cols, _check) if not b] + logger.warning(f"Info for the grid (raster) model was found, but not all. {_missing} was/ were missing") + return False + + return True + + def check_global_crs( srs: osr.SpatialReference, fname: str, @@ -52,6 +101,31 @@ def check_global_crs( ## GIS +def check_grid_exact( + haz, + exp, +): + """_summary_""" + + if not check_vs_srs( + haz.get_srs(), + exp.get_srs(), + ): + logger.error("") + sys.exit() + + gtf1 = [round(_n,2) for _n in haz.get_geotransform()] + gtf2 = [round(_n,2) for _n in exp.get_geotransform()] + + if gtf1 != gtf2: + logger.error() + sys.exit() + + if haz.shape != exp.shape: + logger.error("") + sys.exit() + + def check_internal_srs( source_srs: osr.SpatialReference, fname: str, @@ -155,387 +229,3 @@ def check_hazard_subsets( ## Vulnerability - - -def check_config_data( - cfg: "ConfigReader", - logger: "Log", -): - """General checking of the configurations file""" - for key, item in cfg.values().items(): - if not cfg[key]: - logging.warning( - f"{key} is a missing value. Check the configuration file '{config_path}'.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_crs(crs, crs_type, config_path): - # Check if all Coordinate Reference Systems are correctly assigned. - if isinstance(crs, list): - for c in crs: - try: - assert CRS.from_user_input(c) - except pyproj.exceptions.CRSError as e: - logging.warning(e) - logging.warning( - f"The Coordinate Reference System (CRS) of the {crs_type} is not correctly assigned. Please add a valid CRS as EPSG/ESRI code or a coordinate system in Well-Known Text (WKT) in the configuration file '{config_path}'.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - elif isinstance(crs, str): - try: - assert CRS.from_user_input(crs) - except pyproj.exceptions.CRSError as e: - logging.warning(e) - logging.warning( - f"The Coordinate Reference System (CRS) of the {crs_type} is not correctly assigned. Please add a valid CRS as EPSG/ESRI code or a coordinate system in Well-Known Text (WKT) in the configuration file '{config_path}'.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_required_columns(df_exposure, exposure_path): - list_required_columns = [ - "Object ID", - "Object Name", - "Primary Object Type", - "Secondary Object Type", - "X Coordinate", - "Y Coordinate", - "Extraction Method", - "Damage Function: Structure", - "Damage Function: Content", - "Damage Function: Other", - "First Floor Elevation", - "Ground Elevation", - "Max Potential Damage: Structure", - "Max Potential Damage: Content", - "Max Potential Damage: Other", - "Object-Location Shapefile Path", - "Object-Location Join ID", - "Join Attribute Field", - ] - list_missing = [] - for c in list_required_columns: - if c not in df_exposure.columns: - list_missing.append(c) - - if len(list_missing) > 0: - logging.warning( - f"The following columns are missing from the exposure data: {list_missing}. Please add those to the exposure data CSV file: {exposure_path}.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_damage_function_compliance_and_assign_damage_factor( - damage_function, inundation_depth, object_id -): - obj_id = -999 - - # Raise a warning if the inundation depth exceeds the range of the damage function. - try: - assert inundation_depth * 100 >= damage_function[0] - assert inundation_depth * 100 <= damage_function[1] - - except AssertionError: - # The inundation depth exceeded the limits of the damage function. - obj_id = object_id - - if inundation_depth * 100 < damage_function[0]: - damage_factor = damage_function[2][0] - elif inundation_depth * 100 > damage_function[1]: - damage_factor = damage_function[2][-1] - - else: - index = [i for i in range(damage_function[0], damage_function[1] + 1)].index( - int(round(Decimal(inundation_depth), 2) * 100) - ) - # if damage_function[0] < 0: - # # The damage fraction was not found because of negative water depths in the damage function. - # # Subtract the index range for the negative damage function from the previously calculated index (damage_function[0] is negative) - # index = int(index - (damage_function[0] / (damage_function[0] - damage_function[1]) * len(damage_function[2]))) - - try: - damage_factor = damage_function[2][index] - except IndexError: - logging.warning( - f"Cannot find an appropriate damage fraction for a water depth of {round(Decimal(inundation_depth), 2)} for Object ID {obj_id}." - ) - - return damage_factor, obj_id - - -def report_object_ids_outside_df( - obj_outside_df_structure, obj_outside_df_content, obj_outside_df_other -): - list_object_ids_outside_df = list() - list_object_ids_outside_df.extend( - list(obj_outside_df_structure) - + list(obj_outside_df_content) - + list(obj_outside_df_other) - ) - list_object_ids_outside_df = list(set(list_object_ids_outside_df)) - list_object_ids_outside_df.remove(-999) - list_object_ids_outside_df = str(list_object_ids_outside_df) - - logging.info( - "The inundation depth exceeded the limits of the damage function for the objects with the following Object IDs: {}".format( - list_object_ids_outside_df - ) - ) - - -def check_damage_function_ids( - damage_functions_exposure, damage_functions_config_file, config_path, exposure_path -): - all_df_exposure = list( - damage_functions_exposure["Damage Function: Structure"].unique() - ) - all_df_exposure.extend( - list(damage_functions_exposure["Damage Function: Content"].unique()) - ) - all_df_exposure.extend( - list(damage_functions_exposure["Damage Function: Other"].unique()) - ) - set_df_exposure = set([df for df in all_df_exposure if df == df]) - - if not len(set_df_exposure) == len( - set_df_exposure & set(damage_functions_config_file) - ): - logging.warning( - f"Not all damage functions that are in the exposure input file '{exposure_path}' are defined in the configuration file '{config_path}'. Please add the Damage Functions {set_df_exposure - set(damage_functions_config_file)} in the configuration file.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_uniqueness_object_ids(df_exposure, exposure_path): - if len(df_exposure.index) != len(df_exposure["Object ID"].unique()): - logging.warning( - f"The Object IDs must be unique. Check the exposure input file: {exposure_path}.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - if not all(df_exposure["Object-Location Join ID"].isna()): - for idx, shp_name in ( - df_exposure.loc[~df_exposure["Object-Location Join ID"].isna()] - .groupby("Object-Location Shapefile Path") - .size() - .to_frame() - .reset_index() - .iterrows() - ): - length_index_shp = len( - df_exposure.loc[ - df_exposure["Object-Location Shapefile Path"] - == shp_name["Object-Location Shapefile Path"] - ].index - ) - length_unique_shp_id = len( - [ - n - for n in df_exposure.loc[ - df_exposure["Object-Location Shapefile Path"] - == shp_name["Object-Location Shapefile Path"], - "Object-Location Join ID", - ].unique() - if n == n - ] - ) - if length_index_shp != length_unique_shp_id: - logging.warning( - f"The Object-Location Join IDs must be unique per linked Shapefile. Check the exposure input file: {exposure_path}.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_id_column_shapefile(gdf, col_name, shp_path): - if col_name not in gdf.columns: - logging.warning( - f"The Join Attribute Field {col_name} is not found in shapefile: {shp_path}.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_gdf_z_coord(gdf, name): - if any(gdf.geometry.values.has_z): - logging.warning( - f"The exposure shapefile {name} contains a Z-coordinate. Please supply a Shapefile with only XY coordinates.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_hazard_projection(hazard_paths, coordinate_systems, config_path): - if all(str(hp).endswith(".tif") for hp in hazard_paths): - # The hazard input are geotiff files. - for h, cs in zip(hazard_paths, coordinate_systems): - projection_tiff = gdal.Open(str(h)).GetProjection() - projection_tiff = CRS.from_user_input(projection_tiff) - if projection_tiff.name != cs.name: - logging.warning( - f"The hazard data {h.name} is in Coordinate System {projection_tiff.name} while the following was specified for this dataset: {cs.name}. The hazard data Coordinate System must be correctly defined in the configuration file {config_path}." - ) - - -def check_damages_equals_damage_functions(exposure, exposure_path): - if any( - exposure["Max Potential Damage: Structure"].isna() - != exposure["Damage Function: Structure"].isna() - ): - logging.warning( - f"There are objects that contain a structure damage function but not a max potential structure damage, or vice versa. Please check your exposure input file: {exposure_path}" - ) - if any( - exposure["Max Potential Damage: Content"].isna() - != exposure["Damage Function: Content"].isna() - ): - logging.warning( - f"There are objects that contain a content damage function but not a max potential content damage, or vice versa. Please check your exposure input file: {exposure_path}" - ) - if any( - exposure["Max Potential Damage: Other"].isna() - != exposure["Damage Function: Other"].isna() - ): - logging.warning( - f"There are objects that contain an 'other' damage function but not a max potential 'other' damage, or vice versa. Please check your exposure input file: {exposure_path}" - ) - - -def check_hazard_extent_resolution(list_hazards): - if len(list_hazards) == 1: - return True - check_hazard_extent = [ - gdal.Open(str(haz)).GetGeoTransform() for haz in list_hazards - ] - if len(set(check_hazard_extent)) == 1: - # All hazard have the exact same extents and resolution - return True - else: - return False - - -def check_exposure_modification(exposure_modification_df, exposure_modification_path): - if exposure_modification_df.empty: - logging.warning( - f"An exposure modification file is submitted but contains no data (except for the header). Please check your exposure modification data ({exposure_modification_path}).\n--------------------The simulation has been stopped.--------------------" - ) - exit() - - -def check_shp_paths_and_make_absolute(df_exposure, fiat_input_path): - for shp_path in ( - df_exposure.loc[:, "Object-Location Shapefile Path"].dropna().unique().tolist() - ): - if not Path(shp_path).is_file(): - df_exposure.loc[ - df_exposure["Object-Location Shapefile Path"] == shp_path, - "Object-Location Shapefile Path", - ] = str(Path(fiat_input_path).joinpath(shp_path)) - else: - df_exposure.loc[ - df_exposure["Object-Location Shapefile Path"] == shp_path, - "Object-Location Shapefile Path", - ] = shp_path - return df_exposure - - -def check_geographic_reference(exposure_df, exposure_path): - # Conduct a check to guarantee that each object is assigned with either an X and Y coordinate or an object-location file reference. - x = len(list(exposure_df["X Coordinate"].dropna())) - y = len(list(exposure_df["Y Coordinate"].dropna())) - shps = len(list(exposure_df["Object-Location Shapefile Path"].dropna())) - - if x != y: - logging.warning( - f"The number of X and Y coordinates are not aligned. Please check the X and Y coordinates in the exposure input file: {exposure_path}.\n--------------------The simulation has been stopped.--------------------" - ) - exit() - if x + shps != len(exposure_df.index): - logging.warning( - f"Not all objects have X-Y coordinates. Please check the X and Y coordinates in the exposure input file: {exposure_path}.\n--------------------The simulation has been stopped.--------------------" - ) - exit() - - -def check_input_data(input_data, exposure_path): - # Check if all required fields are in the exposure data. - list_required_fields = [ - "Object ID", - "Extraction Method", - "Damage Function: Structure", - "Damage Function: Content", - "First Floor Elevation", - "Ground Elevation", - "Max Potential Damage: Structure", - "Max Potential Damage: Content", - ] - for fld in list_required_fields: - try: - assert fld in input_data["df_exposure"].columns - except AssertionError: - logging.warning( - f"Column {fld} is required and missing from the exposure input file '{exposure_path}'.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - # Check if all objects have an extraction method of 'centroid' when not providing a shapefile. If not, set it to centroid. - df = input_data["df_exposure"].loc[ - (input_data["df_exposure"]["Extraction Method"] == "AREA"), - ["Object-Location Shapefile Path"], - ] - if df["Object-Location Shapefile Path"].isna().any(): - logging.warning( - "All objects with an 'area' Extraction Method must be linked to a shapefile in the Object-Location Shapefile Path column. The Extraction Method will be set from 'area' to 'centroid' for those objects without an Object-Location Shapefile Path." - ) - idx_area = ( - input_data["df_exposure"] - .loc[input_data["df_exposure"]["Extraction Method"] == "AREA"] - .index - ) - idx_shp = ( - input_data["df_exposure"] - .loc[input_data["df_exposure"]["Object-Location Shapefile Path"].isna()] - .index - ) - input_data["df_exposure"].loc[ - input_data["df_exposure"].index.isin(list(set(idx_area) & set(idx_shp))), - "Extraction Method", - ] = "CENTROID" - - return input_data - - -def check_lower_case_colnames(data, lowercase_columns, correct_name): - if correct_name not in data.columns: - if correct_name.lower() in lowercase_columns: - idx = lowercase_columns.index(correct_name.lower()) - data.rename(columns={data.columns[idx]: correct_name}, inplace=True) - return data - - -def check_correct_columns_names(exposure): - exposure.columns = [c.title() if "ID" not in c else c for c in exposure.columns] - lowercase_columns = [c.lower() for c in exposure.columns] - exposure = check_lower_case_colnames(exposure, lowercase_columns, "Object ID") - exposure = check_lower_case_colnames( - exposure, lowercase_columns, "Object-Location Join ID" - ) - exposure = check_lower_case_colnames(exposure, lowercase_columns, "Buyout (1=yes)") - exposure = check_lower_case_colnames( - exposure, lowercase_columns, "Aggregation Variable: SVI" - ) - return exposure - - -def check_dem_datum(config_data): - if len(set(config_data["inundation_references"])) != 1: - logging.warning( - f"There may only be one type of 'Inundation Reference' for the flood maps. Please check your configuration file '{config_data['config_path']}'.\n--------------------The simulation has been stopped.--------------------" - ) - sys.exit() - - -def check_exposure_within_extent(gdf_exposure): - if gdf_exposure.empty: - logging.warning( - "No exposure objects are within the flood map extent.\n--------------------The simulation has been stopped.--------------------" - ) - exit() diff --git a/src/delft_fiat/main.py b/src/delft_fiat/main.py index 56c58309..b29b06a1 100644 --- a/src/delft_fiat/main.py +++ b/src/delft_fiat/main.py @@ -34,9 +34,14 @@ def from_path( def run(self): """_summary_""" - model = GeomModel(self.cfg) - model.run() - + _models = self.cfg.get_model_type() + if _models[0]: + model = GeomModel(self.cfg) + model.run() + if _models[1]: + model = GridModel(self.cfg) + model.run() + if __name__ == "__main__": pass diff --git a/src/delft_fiat/models/geom.py b/src/delft_fiat/models/geom.py index 863c93a9..a433c50b 100644 --- a/src/delft_fiat/models/geom.py +++ b/src/delft_fiat/models/geom.py @@ -242,8 +242,8 @@ def run( # If there are more than a hazard band in the dataset # Use a pool to execute the calculations - if self._hazard_grid.count > 1: - pcount = min(self.max_threads, self._hazard_grid.count) + if self.hazard_grid.count > 1: + pcount = min(self.max_threads, self.hazard_grid.count) futures = [] with ProcessPoolExecutor(max_workers=pcount) as Pool: _s = time.time() diff --git a/src/delft_fiat/models/grid.py b/src/delft_fiat/models/grid.py index 0cd4aca6..5fe8be15 100644 --- a/src/delft_fiat/models/grid.py +++ b/src/delft_fiat/models/grid.py @@ -37,6 +37,7 @@ def _read_exposure_grid(self): data = open_grid(file, **kw) ## checks logger.info("Executing exposure data checks...") + check_grid_exact(self.hazard_grid, data) self.exposure_grid = data From 65b06ed2079593813f37413b5303fcf476af0fdd Mon Sep 17 00:00:00 2001 From: Brendan Date: Tue, 19 Sep 2023 09:35:48 +0200 Subject: [PATCH 11/28] Black stuff... --- src/delft_fiat/cfg.py | 2 +- src/delft_fiat/check.py | 36 ++++++++++++++++++------------------ src/delft_fiat/main.py | 8 ++++---- 3 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/delft_fiat/cfg.py b/src/delft_fiat/cfg.py index 767204d7..4ec4b2fe 100644 --- a/src/delft_fiat/cfg.py +++ b/src/delft_fiat/cfg.py @@ -1,5 +1,5 @@ from delft_fiat.check import ( - check_config_entries, + check_config_entries, check_config_geom, check_config_grid, ) diff --git a/src/delft_fiat/check.py b/src/delft_fiat/check.py index 9f441237..7e9e6cca 100644 --- a/src/delft_fiat/check.py +++ b/src/delft_fiat/check.py @@ -39,32 +39,32 @@ def check_config_entries( def check_config_geom( - cfg: 'ConfigReader', + cfg: "ConfigReader", ): """_summary_""" _req_cols = [ "exposure.geom.crs", - "exposure.geom.csv", + "exposure.geom.csv", "exposure.geom.file1", ] - _all_geom = [ - item for item in cfg if item.startswith("exposure.geom") - ] + _all_geom = [item for item in cfg if item.startswith("exposure.geom")] if len(_all_geom) == 0: return False - + _check = [item in _all_geom for item in _req_cols] if not all(_check): _missing = [item for item, b in zip(_req_cols, _check) if not b] - logger.warning(f"Info for the geometry model was found, but not all. {_missing} was/ were missing") + logger.warning( + f"Info for the geometry model was found, but not all. {_missing} was/ were missing" + ) return False - + return True def check_config_grid( - cfg: 'ConfigReader', + cfg: "ConfigReader", ): """_summary_""" @@ -72,16 +72,16 @@ def check_config_grid( "exposure.grid.crs", "exposure.grid.file", ] - _all_grid = [ - item for item in cfg if item.startswith("exposure.grid") - ] + _all_grid = [item for item in cfg if item.startswith("exposure.grid")] if len(_all_grid) == 0: return False - + _check = [item in _all_grid for item in _req_cols] if not all(_check): _missing = [item for item, b in zip(_req_cols, _check) if not b] - logger.warning(f"Info for the grid (raster) model was found, but not all. {_missing} was/ were missing") + logger.warning( + f"Info for the grid (raster) model was found, but not all. {_missing} was/ were missing" + ) return False return True @@ -103,19 +103,19 @@ def check_global_crs( ## GIS def check_grid_exact( haz, - exp, + exp, ): """_summary_""" if not check_vs_srs( - haz.get_srs(), + haz.get_srs(), exp.get_srs(), ): logger.error("") sys.exit() - gtf1 = [round(_n,2) for _n in haz.get_geotransform()] - gtf2 = [round(_n,2) for _n in exp.get_geotransform()] + gtf1 = [round(_n, 2) for _n in haz.get_geotransform()] + gtf2 = [round(_n, 2) for _n in exp.get_geotransform()] if gtf1 != gtf2: logger.error() diff --git a/src/delft_fiat/main.py b/src/delft_fiat/main.py index b29b06a1..bd92c9e0 100644 --- a/src/delft_fiat/main.py +++ b/src/delft_fiat/main.py @@ -38,10 +38,10 @@ def run(self): if _models[0]: model = GeomModel(self.cfg) model.run() - if _models[1]: - model = GridModel(self.cfg) - model.run() - + # if _models[1]: + # model = GridModel(self.cfg) + # model.run() + if __name__ == "__main__": pass From 280cd0363bb2a02180ac00e4738da925add0ad4b Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 21 Sep 2023 15:33:13 +0200 Subject: [PATCH 12/28] Added extra check for exposure data; predefine band in worker --- src/delft_fiat/check.py | 37 +++++++++++++++++++++++++++++++++++ src/delft_fiat/main.py | 6 +++--- src/delft_fiat/models/geom.py | 7 ++++++- src/delft_fiat/models/util.py | 6 ++++-- 4 files changed, 50 insertions(+), 6 deletions(-) diff --git a/src/delft_fiat/check.py b/src/delft_fiat/check.py index 7e9e6cca..43ae0168 100644 --- a/src/delft_fiat/check.py +++ b/src/delft_fiat/check.py @@ -1,6 +1,7 @@ from delft_fiat.log import spawn_logger, setup_default_log from delft_fiat.util import NEWLINE_CHAR, deter_type, generic_path_check +import fnmatch import sys from osgeo import gdal from osgeo import osr @@ -226,6 +227,42 @@ def check_hazard_subsets( ## Exposure +def check_exp_columns( + columns: tuple | list, +): + """_summary_""" + + _man_columns = [ + "Object ID", + "Ground Elevation", + "Ground Floor Height", + ] + + _check = [item in columns for item in _man_columns] + if not all(_check): + _missing = [item for item, b in zip(_man_columns, _check) if not b] + logger.error(f"Missing mandatory exposure columns: {_missing}") + sys.exit() + dmg = fnmatch.filter(columns, "Damage Function: *") + dmg_suffix = [item.split(":")[1].strip() for item in dmg] + mpd = fnmatch.filter(columns, "Max Potential Damage: *") + mpd_suffix = [item.split(":")[1].strip() for item in mpd] + + if not dmg: + logger.error("No damage function were given in ") + sys.exit() + + if not mpd: + logger.error("No maximum potential damages were given in ") + sys.exit() + + _check = [item in mpd_suffix for item in dmg_suffix] + if not any(_check): + logger.error("Damage function and maximum potential damage do not have a single match") + sys.exit() + if not all(_check): + _missing = [item for item, b in zip(dmg_suffix, _check) if not b] + logger.warning(f"No every damage function has a corresponding maximum potential damage: {_missing}") ## Vulnerability diff --git a/src/delft_fiat/main.py b/src/delft_fiat/main.py index bd92c9e0..fc273afc 100644 --- a/src/delft_fiat/main.py +++ b/src/delft_fiat/main.py @@ -38,9 +38,9 @@ def run(self): if _models[0]: model = GeomModel(self.cfg) model.run() - # if _models[1]: - # model = GridModel(self.cfg) - # model.run() + if _models[1]: + model = GridModel(self.cfg) + model.run() if __name__ == "__main__": diff --git a/src/delft_fiat/models/geom.py b/src/delft_fiat/models/geom.py index 1bf5bb9e..63171a33 100644 --- a/src/delft_fiat/models/geom.py +++ b/src/delft_fiat/models/geom.py @@ -1,4 +1,8 @@ -from delft_fiat.check import check_internal_srs, check_vs_srs +from delft_fiat.check import ( + check_exp_columns, + check_internal_srs, + check_vs_srs, +) from delft_fiat.gis import geom, overlay from delft_fiat.gis.crs import get_srs_repr from delft_fiat.io import ( @@ -61,6 +65,7 @@ def _read_exposure_data(self): data = open_exp(path, index="Object ID") ##checks logger.info("Executing exposure data checks...") + check_exp_columns(data.columns) ## Information for output _ex = None diff --git a/src/delft_fiat/models/util.py b/src/delft_fiat/models/util.py index 75f4ec1e..53579d3a 100644 --- a/src/delft_fiat/models/util.py +++ b/src/delft_fiat/models/util.py @@ -21,6 +21,8 @@ def geom_worker( ): """_summary_""" + # Extract the hazard band as an object + band = haz[idx] # Setup some metadata _band_name = cfg["hazard.band_names"][idx - 1] _ref = cfg.get("hazard.elevation_reference") @@ -65,10 +67,10 @@ def geom_worker( # Get the hazard data from the exposure geometrie if ft_info[exp._columns["Extraction Method"]].lower() == "area": - res = overlay.clip(haz[idx], haz.get_srs(), haz.get_geotransform(), ft) + res = overlay.clip(band, haz.get_srs(), haz.get_geotransform(), ft) else: res = overlay.pin( - haz[idx], haz.get_geotransform(), geom.point_in_geom(ft) + band, haz.get_geotransform(), geom.point_in_geom(ft) ) # Calculate the inundation From 4b438d660fa0382b327c589ca513379a2b5f95f4 Mon Sep 17 00:00:00 2001 From: Brendan Date: Thu, 21 Sep 2023 15:33:32 +0200 Subject: [PATCH 13/28] black stuff --- src/delft_fiat/check.py | 9 +++++++-- src/delft_fiat/models/geom.py | 2 +- src/delft_fiat/models/util.py | 4 +--- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/delft_fiat/check.py b/src/delft_fiat/check.py index 43ae0168..1159a6fc 100644 --- a/src/delft_fiat/check.py +++ b/src/delft_fiat/check.py @@ -259,10 +259,15 @@ def check_exp_columns( _check = [item in mpd_suffix for item in dmg_suffix] if not any(_check): - logger.error("Damage function and maximum potential damage do not have a single match") + logger.error( + "Damage function and maximum potential damage do not have a single match" + ) sys.exit() if not all(_check): _missing = [item for item, b in zip(dmg_suffix, _check) if not b] - logger.warning(f"No every damage function has a corresponding maximum potential damage: {_missing}") + logger.warning( + f"No every damage function has a corresponding maximum potential damage: {_missing}" + ) + ## Vulnerability diff --git a/src/delft_fiat/models/geom.py b/src/delft_fiat/models/geom.py index 63171a33..736fa426 100644 --- a/src/delft_fiat/models/geom.py +++ b/src/delft_fiat/models/geom.py @@ -1,6 +1,6 @@ from delft_fiat.check import ( check_exp_columns, - check_internal_srs, + check_internal_srs, check_vs_srs, ) from delft_fiat.gis import geom, overlay diff --git a/src/delft_fiat/models/util.py b/src/delft_fiat/models/util.py index 53579d3a..7421eca7 100644 --- a/src/delft_fiat/models/util.py +++ b/src/delft_fiat/models/util.py @@ -69,9 +69,7 @@ def geom_worker( if ft_info[exp._columns["Extraction Method"]].lower() == "area": res = overlay.clip(band, haz.get_srs(), haz.get_geotransform(), ft) else: - res = overlay.pin( - band, haz.get_geotransform(), geom.point_in_geom(ft) - ) + res = overlay.pin(band, haz.get_geotransform(), geom.point_in_geom(ft)) # Calculate the inundation inun, redf = calc_haz( From 9cda8d94e30c214fb4666d8c3d312ea606c26540 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 27 Sep 2023 16:12:29 +0200 Subject: [PATCH 14/28] Proper reducing... --- src/delft_fiat/cfg.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/delft_fiat/cfg.py b/src/delft_fiat/cfg.py index 4ec4b2fe..34b2c62e 100644 --- a/src/delft_fiat/cfg.py +++ b/src/delft_fiat/cfg.py @@ -14,13 +14,23 @@ import os import tomli from osgeo import gdal +from typing import Any class ConfigReader(dict): def __init__( self, file: str, + extra: dict = None, ): + """_summary_""" + + # container for extra + self._build = True + self._extra = {} + if extra is not None: + self._extra.update(extra) + # Set the root directory self.filepath = Path(file) self.path = self.filepath.parent @@ -64,9 +74,27 @@ def __init__( if isinstance(item, str): self[key] = item.lower() + self._build = False + + # (Re)set the extra values + self.update(self._extra) + def __repr__(self): return f"" + def __reduce__(self): + """_summary_""" + + return self.__class__, ( + self.filepath, + self._extra, + ) + + def __setitem__(self, __key: Any, __value: Any): + if not self._build: + self._extra[__key] = __value + super().__setitem__(__key, __value) + def _create_output_dir( self, path: Path | str, From ce6162c972cd5756de519926916ea2efce1bad50 Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 27 Sep 2023 16:12:48 +0200 Subject: [PATCH 15/28] Integration checks, exp column check (geom model) --- src/delft_fiat/check.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/delft_fiat/check.py b/src/delft_fiat/check.py index e9784a98..e42c262b 100644 --- a/src/delft_fiat/check.py +++ b/src/delft_fiat/check.py @@ -119,7 +119,7 @@ def check_grid_exact( gtf2 = [round(_n, 2) for _n in exp.get_geotransform()] if gtf1 != gtf2: - logger.error() + logger.error("") sys.exit() if haz.shape != exp.shape: From b6b1b868f05fd148820036f16a229e08d4eea6bf Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 27 Sep 2023 16:13:04 +0200 Subject: [PATCH 16/28] Extra grid methods (bounds etc..) --- src/delft_fiat/io.py | 67 +++++++++++++++++++++++++++----------------- 1 file changed, 42 insertions(+), 25 deletions(-) diff --git a/src/delft_fiat/io.py b/src/delft_fiat/io.py index d0fb62c7..411085fb 100644 --- a/src/delft_fiat/io.py +++ b/src/delft_fiat/io.py @@ -550,7 +550,7 @@ def __init__( else: ValueError("") - self._recreate_windows() + self.create_windows() def __iter__(self): self.flush() @@ -604,12 +604,6 @@ def __getitem__( chunk = self.src.ReadAsArray(*window) return chunk - def _recreate_windows(self): - self._windows = product( - range(0, self._x, self._chunk[0]), - range(0, self._y, self._chunk[1]), - ) - def _reset_chunking(self): self._l = 0 self._u = 0 @@ -624,13 +618,30 @@ def flush(self): self.src.FlushCache() @property - def chunking(self): + def chunk(self): return self._chunk @property def shape(self): return self._x, self._y + def create_windows(self): + _lu = tuple( + product( + range(0, self._x, self._chunk[0]), + range(0, self._y, self._chunk[1]), + ), + ) + for _l, _u in _lu: + w = min(self._chunk[0], self._x - _l) + h = min(self._chunk[1], self._y - _u) + yield ( + _l, + _u, + w, + h, + ) + def get_metadata_item( self, entry: str, @@ -639,7 +650,7 @@ def get_metadata_item( return self.src.GetMetadataItem(entry) - def set_chunking( + def set_chunk_size( self, chunk: tuple, ): @@ -1070,14 +1081,14 @@ def __getitem__( ): return Grid( self.src.GetRasterBand(oid), - chunk=self._chunk, + chunk=self.chunk, mode=self._mode_str, ) def __reduce__(self): return self.__class__, ( self.path, - self._chunk, + self.chunk, self.subset, self._var_as_band, self._mode_str, @@ -1103,13 +1114,32 @@ def reopen(self): obj = GridSource.__new__( GridSource, self.path, - self._chunk, + self.chunk, self.subset, self._var_as_band, ) obj.__init__(self.path, self._chunk, self.subset, self._var_as_band) return obj + @property + @_BaseIO._check_state + def bounds(self): + """_summary_""" + + _gtf = self.src.GetGeoTransform() + return ( + _gtf[0], + _gtf[0] + _gtf[1] * self.src.RasterXSize, + _gtf[3] + _gtf[5] * self.src.RasterYSize, + _gtf[3], + ) + + @property + def chunk(self): + """_summary_""" + + return self._chunk + @property @_BaseIO._check_state def dtype(self): @@ -1148,19 +1178,6 @@ def create( self.count = nb - @property - @_BaseIO._check_state - def bounds(self): - """_summary_""" - - _gtf = self.src.GetGeoTransform() - return ( - _gtf[0], - _gtf[0] + _gtf[1] * self.src.RasterXSize, - _gtf[3] + _gtf[5] * self.src.RasterYSize, - _gtf[3], - ) - @_BaseIO._check_mode @_BaseIO._check_state def create_band( From 233154a23eb193f032871bb63a93604f7c52517c Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 27 Sep 2023 16:14:11 +0200 Subject: [PATCH 17/28] Extra keywords from settings file --- src/delft_fiat/models/base.py | 73 +++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 30 deletions(-) diff --git a/src/delft_fiat/models/base.py b/src/delft_fiat/models/base.py index 1976b52b..b8fbd767 100644 --- a/src/delft_fiat/models/base.py +++ b/src/delft_fiat/models/base.py @@ -17,6 +17,7 @@ from os import cpu_count from multiprocessing import Manager from osgeo import osr +from pathlib import Path logger = spawn_logger("fiat.model") @@ -28,8 +29,8 @@ def __init__( ): """_summary_""" - self._cfg = cfg - logger.info(f"Using settings from '{self._cfg.filepath}'") + self.cfg = cfg + logger.info(f"Using settings from '{self.cfg.filepath}'") # Declarations self.srs = None @@ -40,7 +41,7 @@ def __init__( self.vulnerability_data = None self._vul_step_size = 0.01 self._rounding = 2 - self._cfg["vulnerability.round"] = self._rounding + self.cfg["vulnerability.round"] = self._rounding self._outhandler = None self._keep_temp = False self._mp_manager = Manager() @@ -51,8 +52,8 @@ def __init__( self._read_hazard_grid() self._read_vulnerability_data() - if "global.keep_temp_files" in self._cfg: - self._keep_temp = self._cfg.get("global.keep_temp_files") + if "global.keep_temp_files" in self.cfg: + self._keep_temp = self.cfg.get("global.keep_temp_files") @abstractmethod def __del__(self): @@ -69,9 +70,16 @@ def _clean_up(self): def _read_hazard_grid(self): """_summary_""" - path = self._cfg.get("hazard.file") + path = self.cfg.get("hazard.file") logger.info(f"Reading hazard data ('{path.name}')") - kw = self._cfg.generate_kwargs("hazard.settings") + # Set the extra arguments from the settings file + kw = {} + kw.update( + self.cfg.generate_kwargs("global.grid"), + ) + kw.update( + self.cfg.generate_kwargs("hazard.settings"), + ) data = open_grid(path, **kw) ## checks logger.info("Executing hazard checks...") @@ -90,7 +98,7 @@ def _read_hazard_grid(self): if _int_srs is not None: logger.info( f"Setting spatial reference of '{path.name}' \ -from '{self._cfg.filepath.name}' ('{get_srs_repr(_int_srs)}')" +from '{self.cfg.filepath.name}' ('{get_srs_repr(_int_srs)}')" ) raise ValueError("") @@ -102,50 +110,55 @@ def _read_hazard_grid(self): ) logger.info(f"Reprojecting '{path.name}' to '{get_srs_repr(self.srs)}'") _resalg = 0 - if "hazard.resampling_method" in self._cfg: - _resalg = self._cfg.get("hazard.resampling_method") + if "hazard.resampling_method" in self.cfg: + _resalg = self.cfg.get("hazard.resampling_method") data = grid.reproject(data, self.srs.ExportToWkt(), _resalg) # check risk return periods - if self._cfg["hazard.risk"]: + if self.cfg["hazard.risk"]: rp = check_hazard_rp_iden( data.get_band_names(), - self._cfg.get("hazard.return_periods"), + self.cfg.get("hazard.return_periods"), path, ) - self._cfg["hazard.return_periods"] = rp + self.cfg["hazard.return_periods"] = rp # Directly calculate the coefficients rp_coef = calc_rp_coef(rp) - self._cfg["hazard.rp_coefficients"] = rp_coef + self.cfg["hazard.rp_coefficients"] = rp_coef + # Set the risk folder for raster calculations + self.cfg["output.path.risk"] = Path( + self.cfg["output.path"], + "rp_damages", + ) # Information for output ns = check_hazard_band_names( data.deter_band_names(), - self._cfg.get("hazard.risk"), - self._cfg.get("hazard.return_periods"), + self.cfg.get("hazard.risk"), + self.cfg.get("hazard.return_periods"), data.count, ) - self._cfg["hazard.band_names"] = ns + self.cfg["hazard.band_names"] = ns # When all is done, add it self.hazard_grid = data def _read_vulnerability_data(self): - path = self._cfg.get("vulnerability.file") + path = self.cfg.get("vulnerability.file") logger.info(f"Reading vulnerability curves ('{path.name}')") index = "water depth" - if "vulnerability.index" in self._cfg: - index = self._cfg.get("vulnerability.index") + if "vulnerability.index" in self.cfg: + index = self.cfg.get("vulnerability.index") data = open_csv(str(path), index=index) ## checks logger.info("Executing vulnerability checks...") # upscale the data (can be done after the checks) - if "vulnerability.step_size" in self._cfg: - self._vul_step_size = self._cfg.get("vulnerability.step_size") + if "vulnerability.step_size" in self.cfg: + self._vul_step_size = self.cfg.get("vulnerability.step_size") self._rounding = deter_dec(self._vul_step_size) - self._cfg["vulnerability.round"] = self._rounding + self.cfg["vulnerability.round"] = self._rounding logger.info( f"Upscaling vulnerability curves, \ @@ -159,7 +172,7 @@ def _set_max_threads(self): """_summary_""" self.max_threads = cpu_count() - _max_threads = self._cfg.get("global.max_threads") + _max_threads = self.cfg.get("global.max_threads") if _max_threads is not None: self.max_threads = min(self.max_threads, _max_threads) @@ -168,14 +181,14 @@ def _set_max_threads(self): def _set_model_srs(self): """_summary_""" - _srs = self._cfg.get("global.crs") - path = self._cfg.get("hazard.file") + _srs = self.cfg.get("global.crs") + path = self.cfg.get("hazard.file") if _srs is not None: self.srs = osr.SpatialReference() self.srs.SetFromUserInput(_srs) else: # Inferring by 'sniffing' - kw = self._cfg.generate_kwargs("hazard.settings") + kw = self.cfg.generate_kwargs("hazard.settings") gm = open_grid( str(path), @@ -184,15 +197,15 @@ def _set_model_srs(self): _srs = gm.get_srs() if _srs is None: - if "hazard.crs" in self._cfg: + if "hazard.crs" in self.cfg: _srs = osr.SpatialReference() - _srs.SetFromUserInput(self._cfg.get("hazard.crs")) + _srs.SetFromUserInput(self.cfg.get("hazard.crs")) self.srs = _srs # Simple check to see if it's not None check_global_crs( self.srs, - self._cfg.filepath.name, + self.cfg.filepath.name, path.name, ) From dfc25d520ea75ee3819d097e84ba860d1708cc2c Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 27 Sep 2023 16:14:25 +0200 Subject: [PATCH 18/28] Extra exposure column check --- src/delft_fiat/models/geom.py | 52 +++++++++++++++++------------------ 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/delft_fiat/models/geom.py b/src/delft_fiat/models/geom.py index 4b665e94..c56d74d6 100644 --- a/src/delft_fiat/models/geom.py +++ b/src/delft_fiat/models/geom.py @@ -53,7 +53,7 @@ def __del__(self): def _clean_up(self): """_summary_""" - _p = self._cfg.get("output.path.tmp") + _p = self.cfg.get("output.path.tmp") for _f in _p.glob("*"): os.unlink(_f) os.rmdir(_p) @@ -61,7 +61,7 @@ def _clean_up(self): def _read_exposure_data(self): """_summary_""" - path = self._cfg.get("exposure.geom.csv") + path = self.cfg.get("exposure.geom.csv") logger.info(f"Reading exposure data ('{path.name}')") data = open_exp(path, index="Object ID") ##checks @@ -70,13 +70,13 @@ def _read_exposure_data(self): ## Information for output _ex = None - if self._cfg["hazard.risk"]: + if self.cfg["hazard.risk"]: _ex = ["Risk (EAD)"] cols = data.create_all_columns( - self._cfg.get("hazard.band_names"), + self.cfg.get("hazard.band_names"), _ex, ) - self._cfg["output.new_columns"] = cols + self.cfg["output.new_columns"] = cols ## When all is done, add it self.exposure_data = data @@ -85,9 +85,9 @@ def _read_exposure_geoms(self): """_summary_""" _d = {} - _found = [item for item in list(self._cfg) if "exposure.geom.file" in item] + _found = [item for item in list(self.cfg) if "exposure.geom.file" in item] for file in _found: - path = self._cfg.get(file) + path = self.cfg.get(file) logger.info( f"Reading exposure geometry '{file.split('.')[-1]}' ('{path.name}')" ) @@ -113,7 +113,7 @@ def _read_exposure_geoms(self): # check if it falls within the extent of the hazard map check_geom_extent( data.bounds, - self._hazard_grid.bounds, + self.hazard_grid.bounds, ) # Add to the dict @@ -129,23 +129,23 @@ def resolve( # Setup some local referenced datasets and metadata _exp = self.exposure_data _gm = self.exposure_geoms - _risk = self._cfg.get("hazard.risk") - _rp_coef = self._cfg.get("hazard.rp_coefficients") + _risk = self.cfg.get("hazard.risk") + _rp_coef = self.cfg.get("hazard.rp_coefficients") # Reverse the _rp_coef to let them coincide with the acquired # values from the temporary files if _rp_coef: _rp_coef.reverse() - _new_cols = self._cfg["output.new_columns"] + _new_cols = self.cfg["output.new_columns"] _files = {} # Define the outgoing file out_csv = "output.csv" - if "output.csv.name" in self._cfg: - out_csv = self._cfg["output.csv.name"] + if "output.csv.name" in self.cfg: + out_csv = self.cfg["output.csv.name"] # Setup the write and write the header of the file writer = BufferTextHandler( - Path(self._cfg["output.path"], out_csv), + Path(self.cfg["output.path"], out_csv), buffer_size=100000, ) header = b"" @@ -155,7 +155,7 @@ def resolve( writer.write(header) # Get all the temporary data paths - _paths = Path(self._cfg.get("output.path.tmp")).glob("*.dat") + _paths = Path(self.cfg.get("output.path.tmp")).glob("*.dat") # Open the temporary files lazy for p in _paths: @@ -169,12 +169,12 @@ def resolve( # Define outgoing dataset out_geom = f"spatial{_add}.gpkg" - if f"output.geom.name{_add}" in self._cfg: - out_geom = self._cfg[f"output.geom.name{_add}"] + if f"output.geom.name{_add}" in self.cfg: + out_geom = self.cfg[f"output.geom.name{_add}"] # Setup the geometry writer geom_writer = GeomMemFileHandler( - Path(self._cfg["output.path"], out_geom), + Path(self.cfg["output.path"], out_geom), self.srs, gm.layer.GetLayerDefn(), ) @@ -241,11 +241,11 @@ def run( """_summary_""" # Get band names for logging - _nms = self._cfg.get("hazard.band_names") + _nms = self.cfg.get("hazard.band_names") # Setup the mp logger for missing stuff _receiver = setup_mp_log( - self._queue, "missing", log_level=2, dst=self._cfg.get("output.path") + self._queue, "missing", log_level=2, dst=self.cfg.get("output.path") ) logger.info("Starting the calculations") @@ -266,7 +266,7 @@ def run( ) fs = Pool.submit( geom_worker, - self._cfg, + self.cfg, self._queue, self.hazard_grid, idx + 1, @@ -289,7 +289,7 @@ def run( p = Process( target=geom_worker, args=( - self._cfg, + self.cfg, self._queue, self.hazard_grid, 1, @@ -309,20 +309,20 @@ def run( _receiver.close_handlers() if _receiver.count > 0: logger.warning( - f"Some objects had missing data. For more info: 'missing.log' in '{self._cfg.get('output.path')}'" + f"Some objects had missing data. For more info: 'missing.log' in '{self.cfg.get('output.path')}'" ) else: os.unlink( - Path(self._cfg.get("output.path"), "missing.log"), + Path(self.cfg.get("output.path"), "missing.log"), ) logger.info("Producing model output from temporary files") # Patch output from the seperate processes back together self.resolve() - logger.info(f"Output generated in: '{self._cfg['output.path']}'") + logger.info(f"Output generated in: '{self.cfg['output.path']}'") if not self._keep_temp: logger.info("Deleting temporary files...") self._clean_up() - logger.info("All done!") + logger.info("Geom calculation are done!") From 144bbee7bd59400b78a916913efd95ddf5080b6c Mon Sep 17 00:00:00 2001 From: Brendan Date: Wed, 27 Sep 2023 16:14:38 +0200 Subject: [PATCH 19/28] Grid model in progess --- src/delft_fiat/models/grid.py | 67 +++++++++++++--- src/delft_fiat/models/util.py | 147 ++++++++++++++++++++++++++++------ 2 files changed, 178 insertions(+), 36 deletions(-) diff --git a/src/delft_fiat/models/grid.py b/src/delft_fiat/models/grid.py index 5fe8be15..38eed65f 100644 --- a/src/delft_fiat/models/grid.py +++ b/src/delft_fiat/models/grid.py @@ -5,7 +5,8 @@ from delft_fiat.models.util import grid_worker_exact from delft_fiat.util import NEWLINE_CHAR -from concurrent.futures import ProcessPoolExecutor +import time +from concurrent.futures import ProcessPoolExecutor, wait from multiprocessing import Process, get_context logger = spawn_logger("fiat.model.grid") @@ -31,9 +32,16 @@ def _clean_up(self): def _read_exposure_grid(self): """_summary_""" - file = self._cfg.get("exposure.grid.file") + file = self.cfg.get("exposure.grid.file") logger.info(f"Reading exposure grid ('{file.name}')") - kw = self._cfg.generate_kwargs("exposure.grid.settings") + # Set the extra arguments from the settings file + kw = {} + kw.update( + self.cfg.generate_kwargs("global.grid"), + ) + kw.update( + self.cfg.generate_kwargs("exposure.grid.settings"), + ) data = open_grid(file, **kw) ## checks logger.info("Executing exposure data checks...") @@ -42,15 +50,54 @@ def _read_exposure_grid(self): self.exposure_grid = data def resolve(): + """_summary_""" + pass def run(self): """_summary_""" - grid_worker_exact( - self._cfg, - self.hazard_grid, - 1, - self.vulnerability_data, - self.exposure_grid, - ) + _nms = self.cfg.get("hazard.band_names") + + if self.hazard_grid.count > 1: + pcount = min(self.max_threads, self.hazard_grid.count) + futures = [] + with ProcessPoolExecutor(max_workers=pcount) as Pool: + _s = time.time() + for idx in range(self.hazard_grid.count): + logger.info( + f"Submitting a job for the calculations in regards to band: '{_nms[idx]}'" + ) + fs = Pool.submit( + grid_worker_exact, + self.cfg, + self.hazard_grid, + idx + 1, + self.vulnerability_data, + self.exposure_grid, + ) + futures.append(fs) + logger.info("Busy...") + wait(futures) + + else: + logger.info(f"Submitting a job for the calculations in a seperate process") + _s = time.time() + p = Process( + target=grid_worker_exact, + args=( + self.cfg, + self.hazard_grid, + 1, + self.vulnerability_data, + self.exposure_grid, + ), + ) + p.start() + logger.info("Busy...") + p.join() + _e = time.time() - _s + + logger.info(f"Calculations time: {round(_e, 2)} seconds") + logger.info(f"Output generated in: '{self.cfg['output.path']}'") + logger.info("Grid calculation are done!") diff --git a/src/delft_fiat/models/util.py b/src/delft_fiat/models/util.py index 7421eca7..a00503a7 100644 --- a/src/delft_fiat/models/util.py +++ b/src/delft_fiat/models/util.py @@ -13,7 +13,7 @@ def geom_worker( cfg: "ConfigReader", queue: "queue.Queue", - haz: "GridSource", + haz: GridSource, idx: int, vul: "Table", exp: "TableLazy", @@ -107,21 +107,30 @@ def geom_worker( def grid_worker_exact( cfg: "ConfigReader", - haz: "GridSource", + haz: GridSource, idx: int, vul: "Table", - exp: "GridSource", + exp: GridSource, ): """_summary_""" - exp_nd = exp[1].nodata - dmf = exp[1].get_metadata_item("damage_function") + # Set some variables for the calculations + exp_bands = [] + write_bands = [] + exp_nds = [] + dmfs = [] + # Extract the hazard band as an object + haz_band = haz[idx] + # Set the output directory + _out = cfg.get("output.path") + if cfg.get("hazard.risk"): + _out = cfg.get("hazard.path.risk") + # Create the outgoing netcdf containing every exposure damages out_src = GridSource( - "C:\\temp\\output.nc", + Path(_out, "output.nc"), mode="w", ) - out_src.create( exp.shape, exp.count, @@ -131,46 +140,132 @@ def grid_worker_exact( out_src.set_srs(exp.get_srs()) out_src.set_geotransform(exp.get_geotransform()) - for _bandn in range(exp.count): - write_band = out_src[_bandn + 1] - write_band.src.SetNoDataValue(exp_nd) - - for (_, h_ch), (_w, e_ch) in zip(haz[idx], exp[_bandn + 1]): - out_ch = full(e_ch.shape, exp_nd) + # td_out = GridSource( + # Path( + # _out, + # "total_damages.nc", + # ), + # mode="w", + # ) + # td_out.create( + # exp.shape, + # 1, + # exp.dtype, + # options=["FORMAT=NC4", "COMPRESS=DEFLATE"], + # ) + # td_out.set_geotransform(exp.get_geotransform()) + # td_out.set_srs(exp.get_srs()) + # td_band = td_out[1] + # td_noval = -0.5 * 2**128 + # td_band.src.SetNoDataValue(td_noval) + + for idx in range(exp.count): + exp_bands.append(exp[idx + 1]) + write_bands.append(out_src[idx + 1]) + exp_nds.append(exp_bands[idx].nodata) + write_bands[idx].src.SetNoDataValue(exp_nds[idx]) + dmfs.append(exp_bands[idx].get_metadata_item("damage_function")) + + for _w in haz_band.create_windows(): + # td_ch = td_band[_w] + h_ch = haz_band[_w] + + for idx, exp_band in enumerate(exp_bands): + e_ch = exp_band[_w] + + out_ch = full(e_ch.shape, exp_nds[idx]) e_ch = ravel(e_ch) - _coords = where(e_ch != exp_nd)[0] + _coords = where(e_ch != exp_nds[idx])[0] if len(_coords) == 0: - write_band.src.WriteArray(out_ch, *_w[:2]) + write_bands[idx].src.WriteArray(out_ch, *_w[:2]) continue e_ch = e_ch[_coords] - h_ch = ravel(h_ch) - h_ch = h_ch[_coords] - _hcoords = where(h_ch != haz[idx].nodata)[0] + h_1d = ravel(h_ch) + h_1d = h_1d[_coords] + _hcoords = where(h_1d != haz_band.nodata)[0] if len(_hcoords) == 0: - write_band.src.WriteArray(out_ch, *_w[:2]) + write_bands[idx].src.WriteArray(out_ch, *_w[:2]) continue _coords = _coords[_hcoords] e_ch = e_ch[_hcoords] - h_ch = h_ch[_hcoords] - h_ch.clip(min(vul.index), max(vul.index)) + h_1d = h_1d[_hcoords] + h_1d.clip(min(vul.index), max(vul.index)) - dmm = [vul[round(float(n), 2), dmf] for n in h_ch] + dmm = [vul[round(float(n), 2), dmfs[idx]] for n in h_1d] e_ch = e_ch * dmm idx2d = unravel_index(_coords, *[exp._chunk]) out_ch[idx2d] = e_ch - write_band.write_chunk(out_ch, _w[:2]) + write_bands[idx].write_chunk(out_ch, _w[:2]) + + # td_1d = td_ch[idx2d] + # td_1d[where(td_1d == td_noval)] = 0 + # td_1d += e_ch + # td_ch[idx2d] = td_1d + + # td_band.write_chunk(td_ch, _w[:2]) + + for _w in write_bands[:]: + w = _w + write_bands.remove(_w) + w.flush() + w = None + + exp_bands = None + # td_band.flush() + # td_band = None + # td_out.flush() + # td_out = None + + # for _bandn in range(exp.count): + # exp_band = exp[_bandn + 1] + # exp_nd = exp_band.nodata + # dmf = exp_band.get_metadata_item("damage_function") + + # write_band = out_src[_bandn + 1] + # write_band.src.SetNoDataValue(exp_nd) - write_band.flush() - write_band = None + # for (_, h_ch), (_w, e_ch) in zip(haz_band, exp_band): + # out_ch = full(e_ch.shape, exp_nd) + # e_ch = ravel(e_ch) + # _coords = where(e_ch != exp_nd)[0] + # if len(_coords) == 0: + # write_band.src.WriteArray(out_ch, *_w[:2]) + # continue + + # e_ch = e_ch[_coords] + # h_ch = ravel(h_ch) + # h_ch = h_ch[_coords] + # _hcoords = where(h_ch != haz_band.nodata)[0] + + # if len(_hcoords) == 0: + # write_band.src.WriteArray(out_ch, *_w[:2]) + # continue + + # _coords = _coords[_hcoords] + # e_ch = e_ch[_hcoords] + # h_ch = h_ch[_hcoords] + # h_ch.clip(min(vul.index), max(vul.index)) + + # dmm = [vul[round(float(n), 2), dmf] for n in h_ch] + # e_ch = e_ch * dmm + + # idx2d = unravel_index(_coords, *[exp._chunk]) + # out_ch[idx2d] = e_ch + + # write_band.write_chunk(out_ch, _w[:2]) + + # write_band.flush() + # write_band = None + # exp_band = None out_src.flush() out_src = None - pass + haz_band = None def grid_worker_loose(): From 3f7c9461ada6bf9286f95e2e2e4ae3b63e8382a8 Mon Sep 17 00:00:00 2001 From: Brendan Date: Sat, 7 Oct 2023 17:23:05 +0200 Subject: [PATCH 20/28] Removed whitespace as result of merge --- src/fiat/io.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/fiat/io.py b/src/fiat/io.py index abdea39d..0fbd15db 100644 --- a/src/fiat/io.py +++ b/src/fiat/io.py @@ -1,4 +1,3 @@ - from fiat.error import DriverNotFoundError from fiat.util import ( DoNotCall, From 94fe0cea59fd6484ae2e3f5d39ae824496cf5830 Mon Sep 17 00:00:00 2001 From: Brendan Date: Sun, 8 Oct 2023 17:28:41 +0200 Subject: [PATCH 21/28] Restructered settings, added damage_function to grid --- .testdata/create_test_data.py | 72 ++++++++++++++++++++++++++++------- 1 file changed, 58 insertions(+), 14 deletions(-) diff --git a/.testdata/create_test_data.py b/.testdata/create_test_data.py index e429e57a..1969ff65 100644 --- a/.testdata/create_test_data.py +++ b/.testdata/create_test_data.py @@ -1,5 +1,6 @@ from fiat.log import Log +import copy import gc import math import os @@ -212,6 +213,7 @@ def create_exposure_grid(): for x, y in product(oneD, oneD): data[x, y] = 2000 + ((x + y) * 100) band.WriteArray(data) + band.SetMetadataItem("damage_function", "struct_1") band.FlushCache() src.FlushCache() @@ -320,13 +322,11 @@ def create_settings(): "elevation_reference": "DEM", }, "exposure": { - "grid": { - "file": "exposure/spatial.nc", - "crs": "EPSG:4326", + "csv": { + "file": "exposure/spatial.csv", }, "geom": { "file1": "exposure/spatial.gpkg", - "csv": "exposure/spatial.csv", "crs": "EPSG:4326", }, }, @@ -345,20 +345,61 @@ def create_settings(): with open(Path(p, "settings.toml"), "wb") as f: tomli_w.dump(doc, f) - doc["output"]["path"] = "output/event_2g" - doc["output"]["geom"]["name2"] = "spatial2.gpkg" - doc["exposure"]["geom"]["file2"] = "exposure/spatial2.gpkg" + doc2g = copy.deepcopy(doc) + doc2g["output"]["path"] = "output/event_2g" + doc2g["output"]["geom"]["name2"] = "spatial2.gpkg" + doc2g["exposure"]["geom"]["file2"] = "exposure/spatial2.gpkg" with open(Path(p, "settings_2g.toml"), "wb") as f: - tomli_w.dump(doc, f) + tomli_w.dump(doc2g, f) - doc["output"]["path"] = "output/event_missing" - del doc["output"]["geom"]["name2"] - del doc["exposure"]["geom"]["file2"] - doc["output"]["geom"]["name1"] = "spatial_missing.gpkg" - doc["exposure"]["geom"]["file1"] = "exposure/spatial_missing.gpkg" + doc_m = copy.deepcopy(doc) + doc_m["output"]["path"] = "output/event_missing" + doc_m["output"]["geom"]["name1"] = "spatial_missing.gpkg" + doc_m["exposure"]["geom"]["file1"] = "exposure/spatial_missing.gpkg" with open(Path(p, "settings_missing.toml"), "wb") as f: + tomli_w.dump(doc_m, f) + + +def create_settings_grid(): + doc = { + "global": { + "crs": "EPSG:4326", + "keep_temp_files": True, + }, + "output": { + "path": "output/event", + "grid": {"name": "output.nc"}, + }, + "hazard": { + "file": "hazard/event_map.nc", + "crs": "EPSG:4326", + "risk": False, + "elevation_reference": "DEM", + }, + "exposure": { + "csv": { + "file": "exposure/spatial.csv", + }, + "grid": { + "file": "exposure/spatial.nc", + "crs": "EPSG:4326", + }, + }, + "vulnerability": { + "file": "vulnerability/vulnerability_curves.csv", + "step_size": 0.01, + }, + "categorical_bins": { + "low": 0.25, + "medium-low": 0.5, + "medium-high": 0.75, + "high": 1, + }, + } + + with open(Path(p, "settings_grid.toml"), "wb") as f: tomli_w.dump(doc, f) @@ -387,13 +428,15 @@ def create_settings_risk(): }, }, "exposure": { + "csv ": { + "file": "exposure/spatial.csv", + }, "grid": { "file": "exposure/spatial.nc", "crs": "EPSG:4326", }, "geom": { "file1": "exposure/spatial.gpkg", - "csv": "exposure/spatial.csv", "crs": "EPSG:4326", }, }, @@ -449,6 +492,7 @@ def log_base(b, x): create_hazard_map() create_risk_map() create_settings() + create_settings_grid() create_settings_risk() create_vulnerability() gc.collect() From 579eca1c1a206f2b1b7db05b1b42859941aeed99 Mon Sep 17 00:00:00 2001 From: Brendan Date: Sun, 8 Oct 2023 17:30:17 +0200 Subject: [PATCH 22/28] Furthered Grid --- src/fiat/check.py | 18 ++++++++++-------- src/fiat/io.py | 13 +++++++++++++ src/fiat/main.py | 6 +++++- src/fiat/models/geom.py | 7 +++++-- src/fiat/models/util.py | 6 +++--- 5 files changed, 36 insertions(+), 14 deletions(-) diff --git a/src/fiat/check.py b/src/fiat/check.py index 5a29d27e..4d99ef2f 100644 --- a/src/fiat/check.py +++ b/src/fiat/check.py @@ -44,18 +44,20 @@ def check_config_geom( ): """_summary_""" - _req_cols = [ + _req_fields = [ + "exposure.csv.file", "exposure.geom.crs", - "exposure.geom.csv", "exposure.geom.file1", ] - _all_geom = [item for item in cfg if item.startswith("exposure.geom")] + _all_geom = [ + item for item in cfg if item.startswith(("exposure.geom", "exposure.csv")) + ] if len(_all_geom) == 0: return False - _check = [item in _all_geom for item in _req_cols] + _check = [item in _all_geom for item in _req_fields] if not all(_check): - _missing = [item for item, b in zip(_req_cols, _check) if not b] + _missing = [item for item, b in zip(_req_fields, _check) if not b] logger.warning( f"Info for the geometry model was found, but not all. {_missing} was/ were missing" ) @@ -69,7 +71,7 @@ def check_config_grid( ): """_summary_""" - _req_cols = [ + _req_fields = [ "exposure.grid.crs", "exposure.grid.file", ] @@ -77,9 +79,9 @@ def check_config_grid( if len(_all_grid) == 0: return False - _check = [item in _all_grid for item in _req_cols] + _check = [item in _all_grid for item in _req_fields] if not all(_check): - _missing = [item for item, b in zip(_req_cols, _check) if not b] + _missing = [item for item, b in zip(_req_fields, _check) if not b] logger.warning( f"Info for the grid (raster) model was found, but not all. {_missing} was/ were missing" ) diff --git a/src/fiat/io.py b/src/fiat/io.py index 0fbd15db..3affb9da 100644 --- a/src/fiat/io.py +++ b/src/fiat/io.py @@ -819,6 +819,19 @@ def bounds(self): return self.layer.GetExtent() + @property + @_BaseIO._check_state + def fields(self): + """_summary_""" + + if self.layer is not None: + _flds = self.layer.GetLayerDefn() + fh = [ + _flds.GetFieldDefn(_i).GetName() for _i in range(_flds.GetFieldCount()) + ] + _flds = None + return fh + @_BaseIO._check_mode @_BaseIO._check_state def add_feature( diff --git a/src/fiat/main.py b/src/fiat/main.py index ac7563e6..c2d4e35a 100644 --- a/src/fiat/main.py +++ b/src/fiat/main.py @@ -1,9 +1,11 @@ from fiat.cfg import ConfigReader -from fiat.log import Log +from fiat.log import spawn_logger from fiat.models import GeomModel, GridModel from pathlib import Path +logger = spawn_logger("fiat.main") + class FIAT: def __init__(self, cfg: ConfigReader): @@ -36,9 +38,11 @@ def run(self): _models = self.cfg.get_model_type() if _models[0]: + logger.info("Setting up geom model..") model = GeomModel(self.cfg) model.run() if _models[1]: + logger.info("Setting up grid model..") model = GridModel(self.cfg) model.run() diff --git a/src/fiat/models/geom.py b/src/fiat/models/geom.py index 6111fa66..f4ad47f4 100644 --- a/src/fiat/models/geom.py +++ b/src/fiat/models/geom.py @@ -61,9 +61,12 @@ def _clean_up(self): def _read_exposure_data(self): """_summary_""" - path = self.cfg.get("exposure.geom.csv") + path = self.cfg.get("exposure.csv.file") logger.info(f"Reading exposure data ('{path.name}')") - data = open_exp(path, index="Object ID") + index_col = "Object ID" + if self.cfg.get("exposure.csv.index") is not None: + index_col = self.cfg.get("exposure.csv.index") + data = open_exp(path, index=index_col) ##checks logger.info("Executing exposure data checks...") check_exp_columns(data.columns) diff --git a/src/fiat/models/util.py b/src/fiat/models/util.py index a5e28a3a..225f2c6e 100644 --- a/src/fiat/models/util.py +++ b/src/fiat/models/util.py @@ -126,6 +126,7 @@ def grid_worker_exact( _out = cfg.get("output.path") if cfg.get("hazard.risk"): _out = cfg.get("hazard.path.risk") + # Create the outgoing netcdf containing every exposure damages out_src = GridSource( Path(_out, "output.nc"), @@ -139,7 +140,7 @@ def grid_worker_exact( ) out_src.set_srs(exp.get_srs()) out_src.set_geotransform(exp.get_geotransform()) - + # Create the outgoing total damage grid # td_out = GridSource( # Path( # _out, @@ -166,9 +167,8 @@ def grid_worker_exact( write_bands[idx].src.SetNoDataValue(exp_nds[idx]) dmfs.append(exp_bands[idx].get_metadata_item("damage_function")) - for _w in haz_band.create_windows(): + for _w, h_ch in haz_band: # td_ch = td_band[_w] - h_ch = haz_band[_w] for idx, exp_band in enumerate(exp_bands): e_ch = exp_band[_w] From 8b0013ac6cc071bccdae7962ffb26afff4b94523 Mon Sep 17 00:00:00 2001 From: Brendan Date: Sun, 8 Oct 2023 17:41:42 +0200 Subject: [PATCH 23/28] Added total damages to the mix --- src/fiat/models/util.py | 98 ++++++++++++----------------------------- 1 file changed, 28 insertions(+), 70 deletions(-) diff --git a/src/fiat/models/util.py b/src/fiat/models/util.py index 225f2c6e..a4499b16 100644 --- a/src/fiat/models/util.py +++ b/src/fiat/models/util.py @@ -141,24 +141,24 @@ def grid_worker_exact( out_src.set_srs(exp.get_srs()) out_src.set_geotransform(exp.get_geotransform()) # Create the outgoing total damage grid - # td_out = GridSource( - # Path( - # _out, - # "total_damages.nc", - # ), - # mode="w", - # ) - # td_out.create( - # exp.shape, - # 1, - # exp.dtype, - # options=["FORMAT=NC4", "COMPRESS=DEFLATE"], - # ) - # td_out.set_geotransform(exp.get_geotransform()) - # td_out.set_srs(exp.get_srs()) - # td_band = td_out[1] - # td_noval = -0.5 * 2**128 - # td_band.src.SetNoDataValue(td_noval) + td_out = GridSource( + Path( + _out, + "total_damages.nc", + ), + mode="w", + ) + td_out.create( + exp.shape, + 1, + exp.dtype, + options=["FORMAT=NC4", "COMPRESS=DEFLATE"], + ) + td_out.set_geotransform(exp.get_geotransform()) + td_out.set_srs(exp.get_srs()) + td_band = td_out[1] + td_noval = -0.5 * 2**128 + td_band.src.SetNoDataValue(td_noval) for idx in range(exp.count): exp_bands.append(exp[idx + 1]) @@ -168,7 +168,7 @@ def grid_worker_exact( dmfs.append(exp_bands[idx].get_metadata_item("damage_function")) for _w, h_ch in haz_band: - # td_ch = td_band[_w] + td_ch = td_band[_w] for idx, exp_band in enumerate(exp_bands): e_ch = exp_band[_w] @@ -202,12 +202,12 @@ def grid_worker_exact( write_bands[idx].write_chunk(out_ch, _w[:2]) - # td_1d = td_ch[idx2d] - # td_1d[where(td_1d == td_noval)] = 0 - # td_1d += e_ch - # td_ch[idx2d] = td_1d + td_1d = td_ch[idx2d] + td_1d[where(td_1d == td_noval)] = 0 + td_1d += e_ch + td_ch[idx2d] = td_1d - # td_band.write_chunk(td_ch, _w[:2]) + td_band.write_chunk(td_ch, _w[:2]) for _w in write_bands[:]: w = _w @@ -216,52 +216,10 @@ def grid_worker_exact( w = None exp_bands = None - # td_band.flush() - # td_band = None - # td_out.flush() - # td_out = None - - # for _bandn in range(exp.count): - # exp_band = exp[_bandn + 1] - # exp_nd = exp_band.nodata - # dmf = exp_band.get_metadata_item("damage_function") - - # write_band = out_src[_bandn + 1] - # write_band.src.SetNoDataValue(exp_nd) - - # for (_, h_ch), (_w, e_ch) in zip(haz_band, exp_band): - # out_ch = full(e_ch.shape, exp_nd) - # e_ch = ravel(e_ch) - # _coords = where(e_ch != exp_nd)[0] - # if len(_coords) == 0: - # write_band.src.WriteArray(out_ch, *_w[:2]) - # continue - - # e_ch = e_ch[_coords] - # h_ch = ravel(h_ch) - # h_ch = h_ch[_coords] - # _hcoords = where(h_ch != haz_band.nodata)[0] - - # if len(_hcoords) == 0: - # write_band.src.WriteArray(out_ch, *_w[:2]) - # continue - - # _coords = _coords[_hcoords] - # e_ch = e_ch[_hcoords] - # h_ch = h_ch[_hcoords] - # h_ch.clip(min(vul.index), max(vul.index)) - - # dmm = [vul[round(float(n), 2), dmf] for n in h_ch] - # e_ch = e_ch * dmm - - # idx2d = unravel_index(_coords, *[exp._chunk]) - # out_ch[idx2d] = e_ch - - # write_band.write_chunk(out_ch, _w[:2]) - - # write_band.flush() - # write_band = None - # exp_band = None + td_band.flush() + td_band = None + td_out = None + out_src.flush() out_src = None From c5dd5dd899d13178a4c461a2d997bbd987964b82 Mon Sep 17 00:00:00 2001 From: Brendan Date: Mon, 9 Oct 2023 13:13:24 +0200 Subject: [PATCH 24/28] Updated testdata for grid and risk --- .testdata/create_test_data.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/.testdata/create_test_data.py b/.testdata/create_test_data.py index 1969ff65..d3d72f04 100644 --- a/.testdata/create_test_data.py +++ b/.testdata/create_test_data.py @@ -379,9 +379,6 @@ def create_settings_grid(): "elevation_reference": "DEM", }, "exposure": { - "csv": { - "file": "exposure/spatial.csv", - }, "grid": { "file": "exposure/spatial.nc", "crs": "EPSG:4326", @@ -428,13 +425,9 @@ def create_settings_risk(): }, }, "exposure": { - "csv ": { + "csv": { "file": "exposure/spatial.csv", }, - "grid": { - "file": "exposure/spatial.nc", - "crs": "EPSG:4326", - }, "geom": { "file1": "exposure/spatial.gpkg", "crs": "EPSG:4326", From 561e87ff84c08bbecb44831b5a7467e3bbf10468 Mon Sep 17 00:00:00 2001 From: Brendan Date: Mon, 9 Oct 2023 13:52:28 +0200 Subject: [PATCH 25/28] Updated with grid testing --- test/conftest.py | 13 ++++++++++++- test/test_cli.py | 3 ++- test/test_gis.py | 1 - test/test_io.py | 2 -- test/test_log.py | 1 - test/test_logic.py | 2 -- test/test_model.py | 43 ++++++++++++++++++++++++++++++++++++------- 7 files changed, 50 insertions(+), 15 deletions(-) diff --git a/test/conftest.py b/test/conftest.py index 600b3a62..3afe82d1 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -11,6 +11,11 @@ def settings_toml_event(): return Path(Path.cwd(), ".testdata", "settings.toml") +@pytest.fixture +def settings_toml_grid(): + return Path(Path.cwd(), ".testdata", "settings_grid.toml") + + @pytest.fixture def settings_toml_risk(): return Path(Path.cwd(), ".testdata", "settings_risk.toml") @@ -23,11 +28,17 @@ def cfg_event(settings_toml_event): @pytest.fixture -def mis_event(): +def cfg_event_mis(): p = Path(Path.cwd(), ".testdata", "settings_missing.toml") return ConfigReader(p) +@pytest.fixture +def cfg_grid(settings_toml_grid): + c = ConfigReader(settings_toml_grid) + return c + + @pytest.fixture def cfg_risk(settings_toml_risk): c = ConfigReader(settings_toml_risk) diff --git a/test/test_cli.py b/test/test_cli.py index 66310b19..5ea7c390 100644 --- a/test/test_cli.py +++ b/test/test_cli.py @@ -1,4 +1,5 @@ -import pytest +import subprocess +from pathlib import Path def test_cli_run(): diff --git a/test/test_gis.py b/test/test_gis.py index f22704b0..748716f3 100644 --- a/test/test_gis.py +++ b/test/test_gis.py @@ -1,6 +1,5 @@ from fiat.gis import overlay, geom, grid -import pytest from numpy import mean diff --git a/test/test_io.py b/test/test_io.py index 62994feb..eee4c8bb 100644 --- a/test/test_io.py +++ b/test/test_io.py @@ -1,7 +1,5 @@ from fiat.io import open_csv -import pytest - def test_tabel(cfg_event): tb = open_csv(cfg_event.get("vulnerability.file")) diff --git a/test/test_log.py b/test/test_log.py index d04ad7f3..af845c90 100644 --- a/test/test_log.py +++ b/test/test_log.py @@ -1,7 +1,6 @@ from fiat.log import CHandler, Log, MessageFormatter, spawn_logger import io -import pytest from pathlib import Path diff --git a/test/test_logic.py b/test/test_logic.py index 9ed839a3..8f48303a 100644 --- a/test/test_logic.py +++ b/test/test_logic.py @@ -1,7 +1,5 @@ from fiat.models.calc import calc_rp_coef, calc_risk -import pytest - def test_calc_risk(): rps = [1, 2, 5, 25, 50, 100] diff --git a/test/test_model.py b/test/test_model.py index 71b458a7..3068ad9f 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -1,7 +1,7 @@ from fiat.io import open_csv -from fiat.models import GeomModel +from fiat.models import GeomModel, GridModel -import pytest +from osgeo import gdal from pathlib import Path @@ -16,9 +16,9 @@ def test_geom_model(tmpdir, cfg_event): assert float(out[3, "Total Damage"]) == 1038 -def test_geom_missing(tmpdir, mis_event): - mis_event.set_output_dir(str(tmpdir)) - model = GeomModel(mis_event) +def test_geom_missing(tmpdir, cfg_event_mis): + cfg_event_mis.set_output_dir(str(tmpdir)) + model = GeomModel(cfg_event_mis) model.run() assert Path(str(tmpdir), "missing.log").exists() @@ -26,5 +26,34 @@ def test_geom_missing(tmpdir, mis_event): assert sum(1 for _ in missing) == 1 -def test_raster_model(): - pass +def test_geom_risk(tmpdir, cfg_risk): + cfg_risk.set_output_dir(str(tmpdir)) + model = GeomModel(cfg_risk) + model.run() + + out = open_csv(Path(str(tmpdir), "output.csv"), index="Object ID") + assert float(out[2, "Damage: Structure (5Y)"]) == 1804 + assert float(out[4, "Total Damage (10Y)"]) == 3840 + assert float(out[3, "Risk (EAD)"]) == 1022.47 + + +def test_grid_model(tmpdir, cfg_grid): + cfg_grid.set_output_dir(str(tmpdir)) + model = GridModel(cfg_grid) + model.run() + + src = gdal.OpenEx( + str(Path(str(tmpdir), "output.nc")), + ) + arr = src.ReadAsArray() + src = None + assert int(arr[2, 4] * 10) == 14091 + assert int(arr[7, 3] * 10) == 8700 + + src = gdal.OpenEx( + str(Path(str(tmpdir), "total_damages.nc")), + ) + arr = src.ReadAsArray() + src = None + assert int(arr[2, 4] * 10) == 14091 + assert int(arr[7, 3] * 10) == 8700 From 86061386acaf78a20186bfca32022aadb8915a27 Mon Sep 17 00:00:00 2001 From: Brendan Date: Mon, 9 Oct 2023 13:52:41 +0200 Subject: [PATCH 26/28] Added code comments --- src/fiat/models/util.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/fiat/models/util.py b/src/fiat/models/util.py index a4499b16..4d1d23ae 100644 --- a/src/fiat/models/util.py +++ b/src/fiat/models/util.py @@ -154,12 +154,14 @@ def grid_worker_exact( exp.dtype, options=["FORMAT=NC4", "COMPRESS=DEFLATE"], ) + # Set the neccesary attributes td_out.set_geotransform(exp.get_geotransform()) td_out.set_srs(exp.get_srs()) td_band = td_out[1] td_noval = -0.5 * 2**128 td_band.src.SetNoDataValue(td_noval) + # Prepare some stuff for looping for idx in range(exp.count): exp_bands.append(exp[idx + 1]) write_bands.append(out_src[idx + 1]) @@ -167,12 +169,15 @@ def grid_worker_exact( write_bands[idx].src.SetNoDataValue(exp_nds[idx]) dmfs.append(exp_bands[idx].get_metadata_item("damage_function")) + # Going trough the chunks for _w, h_ch in haz_band: td_ch = td_band[_w] + # Per exposure band for idx, exp_band in enumerate(exp_bands): e_ch = exp_band[_w] + # See if there is any exposure data out_ch = full(e_ch.shape, exp_nds[idx]) e_ch = ravel(e_ch) _coords = where(e_ch != exp_nds[idx])[0] @@ -180,6 +185,7 @@ def grid_worker_exact( write_bands[idx].src.WriteArray(out_ch, *_w[:2]) continue + # See if there is overlap with the hazard data e_ch = e_ch[_coords] h_1d = ravel(h_ch) h_1d = h_1d[_coords] @@ -189,6 +195,7 @@ def grid_worker_exact( write_bands[idx].src.WriteArray(out_ch, *_w[:2]) continue + # Do the calculations _coords = _coords[_hcoords] e_ch = e_ch[_hcoords] h_1d = h_1d[_hcoords] @@ -200,21 +207,27 @@ def grid_worker_exact( idx2d = unravel_index(_coords, *[exp._chunk]) out_ch[idx2d] = e_ch + # Write it to the band in the outgoing file write_bands[idx].write_chunk(out_ch, _w[:2]) + # Doing the total damages part + # Checking whether it has values or not td_1d = td_ch[idx2d] td_1d[where(td_1d == td_noval)] = 0 td_1d += e_ch td_ch[idx2d] = td_1d + # Write the total damages chunk td_band.write_chunk(td_ch, _w[:2]) + # Flush the cache and dereference for _w in write_bands[:]: w = _w write_bands.remove(_w) w.flush() w = None + # Flush and close all exp_bands = None td_band.flush() td_band = None From caa3d75adc165d7ae1b5c253ee9e83cdf6b3e71f Mon Sep 17 00:00:00 2001 From: Brendan Date: Mon, 9 Oct 2023 13:53:14 +0200 Subject: [PATCH 27/28] Added extra check for damage function in grid --- src/fiat/check.py | 18 ++++++++++++++++++ src/fiat/models/grid.py | 11 ++++++++++- 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/fiat/check.py b/src/fiat/check.py index 4d99ef2f..065536e1 100644 --- a/src/fiat/check.py +++ b/src/fiat/check.py @@ -290,4 +290,22 @@ def check_exp_columns( ) +def check_exp_grid_dmfs( + exp: object, + dmfs: tuple | list, +): + """_summary_""" + + _ef = [_i.get_metadata_item("damage_function") for _i in exp] + _i = None + + _check = [item in dmfs for item in _ef] + if not all(_check): + _missing = [item for item, b in zip(_ef, _check) if not b] + logger.error( + f"Incorrect damage function identifier found in exposure grid: {_missing}", + ) + sys.exit() + + ## Vulnerability diff --git a/src/fiat/models/grid.py b/src/fiat/models/grid.py index e9e376e7..2648df3e 100644 --- a/src/fiat/models/grid.py +++ b/src/fiat/models/grid.py @@ -1,4 +1,7 @@ -from fiat.check import * +from fiat.check import ( + check_exp_grid_dmfs, + check_grid_exact, +) from fiat.io import BufferTextHandler, open_grid from fiat.log import spawn_logger from fiat.models.base import BaseModel @@ -45,7 +48,13 @@ def _read_exposure_grid(self): data = open_grid(file, **kw) ## checks logger.info("Executing exposure data checks...") + # Check exact overlay of exposure and hazard check_grid_exact(self.hazard_grid, data) + # Check if all damage functions are correct + check_exp_grid_dmfs( + data, + self.vulnerability_data.columns, + ) self.exposure_grid = data From c62b03bdaaa8bb96353010ed4b2280fbb70356ed Mon Sep 17 00:00:00 2001 From: Brendan Date: Mon, 9 Oct 2023 14:16:57 +0200 Subject: [PATCH 28/28] Sort temporary files --- src/fiat/models/geom.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fiat/models/geom.py b/src/fiat/models/geom.py index f4ad47f4..2134400f 100644 --- a/src/fiat/models/geom.py +++ b/src/fiat/models/geom.py @@ -161,7 +161,7 @@ def resolve( _paths = Path(self.cfg.get("output.path.tmp")).glob("*.dat") # Open the temporary files lazy - for p in _paths: + for p in sorted(_paths): _d = open_csv(p, index=_exp.meta["index_name"], large=True) _files[p.stem] = _d _d = None