From 830a181abb258dd25b5a2b71b294ae425d90f6ed Mon Sep 17 00:00:00 2001 From: songololo Date: Tue, 21 Oct 2025 17:41:00 +0100 Subject: [PATCH 1/3] f32 --- demos/athens-demo.py | 1 + tests/test_common.py | 8 +- umep/class_configs.py | 114 ++-- umep/common.py | 51 +- .../Solweig_2025a_calc_forprocessing.py | 555 ++++++++++++++---- umep/functions/SOLWEIGpython/Solweig_run.py | 77 ++- .../SOLWEIGpython/UTCI_calculations.py | 476 ++++++++------- .../functions/SOLWEIGpython/solweig_runner.py | 21 +- .../SOLWEIGpython/solweig_runner_core.py | 28 +- umep/functions/dailyshading.py | 84 +-- umep/shadow_generator_algorithm.py | 18 +- umep/skyviewfactor_algorithm.py | 42 +- umep/solweig_algorithm.py | 48 +- umep/wall_heightaspect_algorithm.py | 6 +- 14 files changed, 982 insertions(+), 547 deletions(-) diff --git a/demos/athens-demo.py b/demos/athens-demo.py index 818e1b4..723f81b 100644 --- a/demos/athens-demo.py +++ b/demos/athens-demo.py @@ -42,6 +42,7 @@ cdsm_rast, cdsm_transf.to_gdal(), CRS.from_epsg(working_crs).to_wkt(), + coerce_f64_to_f32=True, ) # %% # wall info for SOLWEIG diff --git a/tests/test_common.py b/tests/test_common.py index 1d84c37..60852a2 100644 --- a/tests/test_common.py +++ b/tests/test_common.py @@ -17,8 +17,8 @@ def test_save_and_load_raster_roundtrip(tmp_path): crs_wkt = pyproj.CRS.from_epsg(4326).to_wkt() # save and load - save_raster(str(out), data, trf, crs_wkt, no_data_val=-9999) - rast, trf_out, crs_out, nodata = load_raster(str(out)) + save_raster(str(out), data, trf, crs_wkt, no_data_val=-9999, coerce_f64_to_f32=True) + rast, trf_out, crs_out, nodata = load_raster(str(out), coerce_f64_to_f32=True) np.testing.assert_array_equal(rast, data) assert isinstance(trf_out, list) and len(trf_out) == 6 @@ -34,12 +34,12 @@ def test_load_raster_with_bbox(tmp_path): data = np.arange(100, dtype=np.float32).reshape(10, 10) trf = _make_gt(10, 10) crs_wkt = pyproj.CRS.from_epsg(4326).to_wkt() - save_raster(str(out), data, trf, crs_wkt, -9999) + save_raster(str(out), data, trf, crs_wkt, -9999, coerce_f64_to_f32=True) # bbox in spatial coords: [minx, miny, maxx, maxy] # For our geotransform bounds = [0,0,10,10] bbox = [2, 2, 5, 5] - rast_crop, trf_crop, crs_crop, nd = load_raster(str(out), bbox=bbox) + rast_crop, trf_crop, crs_crop, nd = load_raster(str(out), bbox=bbox, coerce_f64_to_f32=True) # Expected slice computed from implementation mapping assert rast_crop.shape == (3, 3) diff --git a/umep/class_configs.py b/umep/class_configs.py index 49267ea..fece9fc 100644 --- a/umep/class_configs.py +++ b/umep/class_configs.py @@ -234,9 +234,9 @@ def __init__( data_len = len(self.YYYY) self.dectime = self.DOY + self.hours / 24 + self.minu / (60 * 24.0) if data_len == 1: - halftimestepdec = 0 + halftimestepdec = 0.0 else: - halftimestepdec = (self.dectime[1] - self.dectime[0]) / 2.0 + halftimestepdec = float((self.dectime[1] - self.dectime[0]) / 2.0) # convert from f32 time = { "sec": 0, "UTC": UTC, @@ -244,15 +244,15 @@ def __init__( sunmaximum = 0.0 # initialize arrays - self.altitude = np.empty(data_len) - self.azimuth = np.empty(data_len) - self.zen = np.empty(data_len) - self.jday = np.empty(data_len) - self.leafon = np.empty(data_len) - self.psi = np.empty(data_len) - self.altmax = np.empty(data_len) - self.Twater = np.empty(data_len) - self.CI = np.empty(data_len) + self.altitude = np.empty(data_len, dtype=np.float32) + self.azimuth = np.empty(data_len, dtype=np.float32) + self.zen = np.empty(data_len, dtype=np.float32) + self.jday = np.empty(data_len, dtype=np.float32) + self.leafon = np.empty(data_len, dtype=np.float32) + self.psi = np.empty(data_len, dtype=np.float32) + self.altmax = np.empty(data_len, dtype=np.float32) + self.Twater = np.empty(data_len, dtype=np.float32) + self.CI = np.empty(data_len, dtype=np.float32) sunmax = dict() @@ -395,22 +395,32 @@ def __init__(self, model_configs: SolweigConfig): with zipfile.ZipFile(svf_path_str, "r") as zip_ref: zip_ref.extractall(in_path_str) # Load SVF rasters - self.svf, _, _, _ = common.load_raster(in_path_str + "/" + "svf.tif") - self.svf_east, _, _, _ = common.load_raster(in_path_str + "/" + "svfE.tif") - self.svf_south, _, _, _ = common.load_raster(in_path_str + "/" + "svfS.tif") - self.svf_west, _, _, _ = common.load_raster(in_path_str + "/" + "svfW.tif") - self.svf_north, _, _, _ = common.load_raster(in_path_str + "/" + "svfN.tif") + self.svf, _, _, _ = common.load_raster(in_path_str + "/" + "svf.tif", coerce_f64_to_f32=True) + self.svf_east, _, _, _ = common.load_raster(in_path_str + "/" + "svfE.tif", coerce_f64_to_f32=True) + self.svf_south, _, _, _ = common.load_raster(in_path_str + "/" + "svfS.tif", coerce_f64_to_f32=True) + self.svf_west, _, _, _ = common.load_raster(in_path_str + "/" + "svfW.tif", coerce_f64_to_f32=True) + self.svf_north, _, _, _ = common.load_raster(in_path_str + "/" + "svfN.tif", coerce_f64_to_f32=True) if model_configs.use_veg_dem: - self.svf_veg, _, _, _ = common.load_raster(in_path_str + "/" + "svfveg.tif") - self.svf_veg_east, _, _, _ = common.load_raster(in_path_str + "/" + "svfEveg.tif") - self.svf_veg_south, _, _, _ = common.load_raster(in_path_str + "/" + "svfSveg.tif") - self.svf_veg_west, _, _, _ = common.load_raster(in_path_str + "/" + "svfWveg.tif") - self.svf_veg_north, _, _, _ = common.load_raster(in_path_str + "/" + "svfNveg.tif") - self.svf_veg_blocks_bldg_sh, _, _, _ = common.load_raster(in_path_str + "/" + "svfaveg.tif") - self.svf_veg_blocks_bldg_sh_east, _, _, _ = common.load_raster(in_path_str + "/" + "svfEaveg.tif") - self.svf_veg_blocks_bldg_sh_south, _, _, _ = common.load_raster(in_path_str + "/" + "svfSaveg.tif") - self.svf_veg_blocks_bldg_sh_west, _, _, _ = common.load_raster(in_path_str + "/" + "svfWaveg.tif") - self.svf_veg_blocks_bldg_sh_north, _, _, _ = common.load_raster(in_path_str + "/" + "svfNaveg.tif") + self.svf_veg, _, _, _ = common.load_raster(in_path_str + "/" + "svfveg.tif", coerce_f64_to_f32=True) + self.svf_veg_east, _, _, _ = common.load_raster(in_path_str + "/" + "svfEveg.tif", coerce_f64_to_f32=True) + self.svf_veg_south, _, _, _ = common.load_raster(in_path_str + "/" + "svfSveg.tif", coerce_f64_to_f32=True) + self.svf_veg_west, _, _, _ = common.load_raster(in_path_str + "/" + "svfWveg.tif", coerce_f64_to_f32=True) + self.svf_veg_north, _, _, _ = common.load_raster(in_path_str + "/" + "svfNveg.tif", coerce_f64_to_f32=True) + self.svf_veg_blocks_bldg_sh, _, _, _ = common.load_raster( + in_path_str + "/" + "svfaveg.tif", coerce_f64_to_f32=True + ) + self.svf_veg_blocks_bldg_sh_east, _, _, _ = common.load_raster( + in_path_str + "/" + "svfEaveg.tif", coerce_f64_to_f32=True + ) + self.svf_veg_blocks_bldg_sh_south, _, _, _ = common.load_raster( + in_path_str + "/" + "svfSaveg.tif", coerce_f64_to_f32=True + ) + self.svf_veg_blocks_bldg_sh_west, _, _, _ = common.load_raster( + in_path_str + "/" + "svfWaveg.tif", coerce_f64_to_f32=True + ) + self.svf_veg_blocks_bldg_sh_north, _, _, _ = common.load_raster( + in_path_str + "/" + "svfNaveg.tif", coerce_f64_to_f32=True + ) logger.info("Vegetation SVF data loaded.") else: self.svf_veg = np.ones_like(self.svf) @@ -538,7 +548,9 @@ def __init__( amax_local_perc: float = 99.9, ): # Load DSM - self.dsm, self.trf_arr, self.crs_wkt, self.nd_val = common.load_raster(model_configs.dsm_path, bbox=None) + self.dsm, self.trf_arr, self.crs_wkt, self.nd_val = common.load_raster( + model_configs.dsm_path, bbox=None, coerce_f64_to_f32=True + ) logger.info("DSM loaded from %s", model_configs.dsm_path) self.scale = 1 / self.trf_arr[1] self.rows = self.dsm.shape[0] @@ -552,7 +564,9 @@ def __init__( # WALLS # heights - self.wallheight, wh_trf, wh_crs, _ = common.load_raster(model_configs.wh_path, bbox=None) + self.wallheight, wh_trf, wh_crs, _ = common.load_raster( + model_configs.wh_path, bbox=None, coerce_f64_to_f32=True + ) if not self.wallheight.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for wall heights and DSM.") if not np.allclose(self.trf_arr, wh_trf): @@ -561,7 +575,9 @@ def __init__( raise ValueError("Mismatching CRS for wall heights and DSM.") logger.info("Wall heights loaded") # aspects - self.wallaspect, wa_trf, wa_crs, _ = common.load_raster(model_configs.wa_path, bbox=None) + self.wallaspect, wa_trf, wa_crs, _ = common.load_raster( + model_configs.wa_path, bbox=None, coerce_f64_to_f32=True + ) if not self.wallaspect.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for wall aspects and DSM.") if not np.allclose(self.trf_arr, wa_trf): @@ -574,7 +590,7 @@ def __init__( # TODO: Is DEM always provided? if model_configs.dem_path: dem_path_str = str(common.check_path(model_configs.dem_path)) - self.dem, dem_trf, dem_crs, dem_nd_val = common.load_raster(dem_path_str, bbox=None) + self.dem, dem_trf, dem_crs, dem_nd_val = common.load_raster(dem_path_str, bbox=None, coerce_f64_to_f32=True) if not self.dem.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for DEM and CDSM.") if dem_crs is not None and dem_crs != self.crs_wkt: @@ -592,7 +608,9 @@ def __init__( # Vegetation if model_configs.use_veg_dem: - self.cdsm, vegdsm_trf, vegdsm_crs, _ = common.load_raster(model_configs.cdsm_path, bbox=None) + self.cdsm, vegdsm_trf, vegdsm_crs, _ = common.load_raster( + model_configs.cdsm_path, bbox=None, coerce_f64_to_f32=True + ) if not self.cdsm.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for DSM and CDSM.") if vegdsm_crs is not None and vegdsm_crs != self.crs_wkt: @@ -602,7 +620,9 @@ def __init__( logger.info("Vegetation DSM loaded from %s", model_configs.cdsm_path) # Tree DSM if model_configs.tdsm_path: - self.tdsm, vegdsm2_trf, vegdsm2_crs, _ = common.load_raster(model_configs.tdsm_path, bbox=None) + self.tdsm, vegdsm2_trf, vegdsm2_crs, _ = common.load_raster( + model_configs.tdsm_path, bbox=None, coerce_f64_to_f32=True + ) if not self.tdsm.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for DSM and CDSM.") if vegdsm2_crs is not None and vegdsm2_crs != self.crs_wkt: @@ -635,7 +655,7 @@ def __init__( ) else: logger.info("Vegetation DEM not used; vegetation arrays set to None.") - self.bush = np.zeros([self.rows, self.cols]) + self.bush = np.zeros([self.rows, self.cols], dtype=np.float32) self.svfbuveg = svf_data.svf common.save_raster( @@ -644,6 +664,7 @@ def __init__( self.trf_arr, self.crs_wkt, self.nd_val, + coerce_f64_to_f32=True, ) common.save_raster( model_configs.output_dir + "/input-cdsm.tif", @@ -651,6 +672,7 @@ def __init__( self.trf_arr, self.crs_wkt, self.nd_val, + coerce_f64_to_f32=True, ) common.save_raster( model_configs.output_dir + "/input-tdsm.tif", @@ -658,6 +680,7 @@ def __init__( self.trf_arr, self.crs_wkt, self.nd_val, + coerce_f64_to_f32=True, ) common.save_raster( model_configs.output_dir + "/input-svfbuveg.tif", @@ -665,6 +688,7 @@ def __init__( self.trf_arr, self.crs_wkt, self.nd_val, + coerce_f64_to_f32=True, ) common.save_raster( model_configs.output_dir + "/input-bush.tif", @@ -672,12 +696,13 @@ def __init__( self.trf_arr, self.crs_wkt, self.nd_val, + coerce_f64_to_f32=True, ) # Land cover if model_configs.use_landcover: lc_path_str = str(common.check_path(model_configs.lc_path)) - self.lcgrid, lc_trf, lc_crs, _ = common.load_raster(lc_path_str, bbox=None) + self.lcgrid, lc_trf, lc_crs, _ = common.load_raster(lc_path_str, bbox=None, coerce_f64_to_f32=True) if not self.lcgrid.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for land cover and DSM.") if lc_crs is not None and lc_crs != self.crs_wkt: @@ -720,6 +745,7 @@ def __init__( self.trf_arr, self.crs_wkt, self.nd_val, + coerce_f64_to_f32=True, ) logger.info("Buildings raster saved to %s/buildings.tif", model_configs.output_dir) @@ -781,7 +807,7 @@ def __init__( self.asvf = np.arccos(np.sqrt(svf_data.svf)) # Empty array for steradians - self.steradians = np.zeros(self.shmat.shape[2]) + self.steradians = np.zeros(self.shmat.shape[2], dtype=np.float32) else: # no anisotropic sky # downstream functions only access these if use_aniso is True @@ -823,13 +849,13 @@ def __init__(self, use_landcover: bool, model_params, raster_data: RasterData): This is a vectorized version that avoids looping over pixels. """ # Initialization of maps - self.Knight = np.zeros((raster_data.rows, raster_data.cols)) - self.Tgmap1 = np.zeros((raster_data.rows, raster_data.cols)) - self.Tgmap1E = np.zeros((raster_data.rows, raster_data.cols)) - self.Tgmap1S = np.zeros((raster_data.rows, raster_data.cols)) - self.Tgmap1W = np.zeros((raster_data.rows, raster_data.cols)) - self.Tgmap1N = np.zeros((raster_data.rows, raster_data.cols)) - self.TgOut1 = np.zeros((raster_data.rows, raster_data.cols)) + self.Knight = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) + self.Tgmap1 = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) + self.Tgmap1E = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) + self.Tgmap1S = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) + self.Tgmap1W = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) + self.Tgmap1N = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) + self.TgOut1 = np.zeros((raster_data.rows, raster_data.cols), dtype=np.float32) # Set up the Tg maps based on whether land cover is used if use_landcover is False: @@ -955,7 +981,7 @@ def __init__( self.voxelMaps = None self.voxelTable = None self.timeStep = 0 - self.walls_scheme = np.ones((raster_data.rows, raster_data.cols)) * 10.0 - self.dirwalls_scheme = np.ones((raster_data.rows, raster_data.cols)) * 10.0 + self.walls_scheme = np.ones((raster_data.rows, raster_data.cols), dtype=np.float32) * 10.0 + self.dirwalls_scheme = np.ones((raster_data.rows, raster_data.cols), dtype=np.float32) * 10.0 self.met_for_xarray = None logger.info("Wall scheme not used; default wall data initialized.") diff --git a/umep/common.py b/umep/common.py index f8a7b95..a756cbd 100644 --- a/umep/common.py +++ b/umep/common.py @@ -61,8 +61,28 @@ def check_path(path_str: str | Path, make_dir: bool = False) -> Path: def save_raster( - out_path_str: str, data_arr: np.ndarray, trf_arr: list[float], crs_wkt: str, no_data_val: float = -9999 + out_path_str: str, + data_arr: np.ndarray, + trf_arr: list[float], + crs_wkt: str, + no_data_val: float = -9999, + coerce_f64_to_f32: bool = True, ): + """ + Save raster to GeoTIFF. + + Args: + out_path_str: Output file path + data_arr: 2D numpy array to save + trf_arr: GDAL-style geotransform [top_left_x, pixel_width, rotation, top_left_y, rotation, pixel_height] + crs_wkt: CRS in WKT format + no_data_val: No-data value to use + coerce_f64_to_f32: If True, convert float64 arrays to float32 before saving (default: True for memory efficiency) + """ + # Only convert float64 to float32, leave ints/bools unchanged + if coerce_f64_to_f32 and data_arr.dtype == np.float64: + data_arr = data_arr.astype(np.float32) + attempts = 2 while attempts > 0: attempts -= 1 @@ -119,8 +139,20 @@ def save_raster( def load_raster( - path_str: str, bbox: list[int] | None = None, band: int = 0 + path_str: str, bbox: list[int] | None = None, band: int = 0, coerce_f64_to_f32: bool = True ) -> tuple[np.ndarray, list[float], str | None, float | None]: + """ + Load raster, optionally crop to bbox. + + Args: + path_str: Path to raster file + bbox: Optional bounding box [minx, miny, maxx, maxy] + band: Band index to read (0-based) + coerce_f64_to_f32: If True, coerce array to float32 (default: True for memory efficiency) + + Returns: + Tuple of (array, transform, crs_wkt, no_data_value) + """ # Load raster, optionally crop to bbox path = check_path(path_str, make_dir=False) if not path.exists(): @@ -150,9 +182,15 @@ def load_raster( if rast.ndim == 3: if band < 0 or band >= rast.shape[0]: raise IndexError(f"Requested band {band} out of range; raster has {rast.shape[0]} band(s)") - rast_arr = rast[band].astype(float) + rast_arr = rast[band] + # Only convert float64 to float32, leave ints/bools unchanged + if coerce_f64_to_f32 and rast_arr.dtype == np.float64: + rast_arr = rast_arr.astype(np.float32) else: - rast_arr = rast.astype(float) + rast_arr = rast + # Only convert float64 to float32, leave ints/bools unchanged + if coerce_f64_to_f32 and rast_arr.dtype == np.float64: + rast_arr = rast_arr.astype(np.float32) else: dataset = gdal.Open(str(path)) if dataset is None: @@ -164,7 +202,10 @@ def load_raster( if rb is None: dataset = None raise IndexError(f"Requested band {band} out of range in GDAL dataset") - rast_arr = rb.ReadAsArray().astype(float) + rast_arr = rb.ReadAsArray() + # Only convert float64 to float32, leave ints/bools unchanged + if coerce_f64_to_f32 and rast_arr.dtype == np.float64: + rast_arr = rast_arr.astype(np.float32) no_data_val = rb.GetNoDataValue() if bbox is not None: min_x, min_y, max_x, max_y = bbox diff --git a/umep/functions/SOLWEIGpython/Solweig_2025a_calc_forprocessing.py b/umep/functions/SOLWEIGpython/Solweig_2025a_calc_forprocessing.py index 1abcf2f..fc63ce2 100644 --- a/umep/functions/SOLWEIGpython/Solweig_2025a_calc_forprocessing.py +++ b/umep/functions/SOLWEIGpython/Solweig_2025a_calc_forprocessing.py @@ -1,48 +1,133 @@ -from __future__ import absolute_import +from copy import deepcopy import numpy as np -from .daylen import daylen + from ...util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import clearnessindex_2013b +from ...util.SEBESOLWEIGCommonFiles.create_patches import create_patches from ...util.SEBESOLWEIGCommonFiles.diffusefraction import diffusefraction +from ...util.SEBESOLWEIGCommonFiles.Perez_v3 import Perez_v3 from ...util.SEBESOLWEIGCommonFiles.shadowingfunction_wallheight_13 import shadowingfunction_wallheight_13 from ...util.SEBESOLWEIGCommonFiles.shadowingfunction_wallheight_23 import shadowingfunction_wallheight_23 -from .gvf_2018a import gvf_2018a +from .anisotropic_sky import anisotropic_sky as ani_sky from .cylindric_wedge import cylindric_wedge -from .TsWaveDelay_2015a import TsWaveDelay_2015a -from .Kup_veg_2015a import Kup_veg_2015a +from .daylen import daylen +from .gvf_2018a import gvf_2018a + # from .Lside_veg_v2015a import Lside_veg_v2015a # from .Kside_veg_v2019a import Kside_veg_v2019a from .Kside_veg_v2022a import Kside_veg_v2022a -from ...util.SEBESOLWEIGCommonFiles.Perez_v3 import Perez_v3 -from ...util.SEBESOLWEIGCommonFiles.create_patches import create_patches +from .Kup_veg_2015a import Kup_veg_2015a # Anisotropic longwave -from .Lcyl_v2022a import Lcyl_v2022a from .Lside_veg_v2022a import Lside_veg_v2022a -from .anisotropic_sky import anisotropic_sky as ani_sky from .patch_radiation import patch_steradians -from copy import deepcopy -import time +from .TsWaveDelay_2015a import TsWaveDelay_2015a # Wall surface temperature scheme from .wall_surface_temperature import wall_surface_temperature -def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, svfveg, svfNveg, svfEveg, svfSveg, - svfWveg, svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg, vegdem, vegdem2, albedo_b, absK, absL, - ewall, Fside, Fup, Fcyl, altitude, azimuth, zen, jday, usevegdem, onlyglobal, buildings, location, psi, - landcover, lc_grid, dectime, altmax, dirwalls, walls, cyl, elvis, Ta, RH, radG, radD, radI, P, - amaxvalue, bush, Twater, TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST, - TmaxLST_wall, first, second, svfalfa, svfbuveg, firstdaytime, timeadd, timestepdec, Tgmap1, - Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, CI, TgOut1, diffsh, shmat, vegshmat, vbshvegshmat, anisotropic_sky, asvf, patch_option, - voxelMaps, voxelTable, ws, wallScheme, timeStep, steradians, walls_scheme, dirwalls_scheme): - -#def Solweig_2021a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, svfveg, svfNveg, svfEveg, svfSveg, -# svfWveg, svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg, vegdem, vegdem2, albedo_b, absK, absL, -# ewall, Fside, Fup, Fcyl, altitude, azimuth, zen, jday, usevegdem, onlyglobal, buildings, location, psi, -# landcover, lc_grid, dectime, altmax, dirwalls, walls, cyl, elvis, Ta, RH, radG, radD, radI, P, -# amaxvalue, bush, Twater, TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST, -# TmaxLST_wall, first, second, svfalfa, svfbuveg, firstdaytime, timeadd, timestepdec, Tgmap1, -# Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, CI, TgOut1, diffsh, ani): + +def Solweig_2025a_calc( + i, + dsm, + scale, + rows, + cols, + svf, + svfN, + svfW, + svfE, + svfS, + svfveg, + svfNveg, + svfEveg, + svfSveg, + svfWveg, + svfaveg, + svfEaveg, + svfSaveg, + svfWaveg, + svfNaveg, + vegdem, + vegdem2, + albedo_b, + absK, + absL, + ewall, + Fside, + Fup, + Fcyl, + altitude, + azimuth, + zen, + jday, + usevegdem, + onlyglobal, + buildings, + location, + psi, + landcover, + lc_grid, + dectime, + altmax, + dirwalls, + walls, + cyl, + elvis, + Ta, + RH, + radG, + radD, + radI, + P, + amaxvalue, + bush, + Twater, + TgK, + Tstart, + alb_grid, + emis_grid, + TgK_wall, + Tstart_wall, + TmaxLST, + TmaxLST_wall, + first, + second, + svfalfa, + svfbuveg, + firstdaytime, + timeadd, + timestepdec, + Tgmap1, + Tgmap1E, + Tgmap1S, + Tgmap1W, + Tgmap1N, + CI, + TgOut1, + diffsh, + shmat, + vegshmat, + vbshvegshmat, + anisotropic_sky, + asvf, + patch_option, + voxelMaps, + voxelTable, + ws, + wallScheme, + timeStep, + steradians, + walls_scheme, + dirwalls_scheme, +): + # def Solweig_2021a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, svfveg, svfNveg, svfEveg, svfSveg, + # svfWveg, svfaveg, svfEaveg, svfSaveg, svfWaveg, svfNaveg, vegdem, vegdem2, albedo_b, absK, absL, + # ewall, Fside, Fup, Fcyl, altitude, azimuth, zen, jday, usevegdem, onlyglobal, buildings, location, psi, + # landcover, lc_grid, dectime, altmax, dirwalls, walls, cyl, elvis, Ta, RH, radG, radD, radI, P, + # amaxvalue, bush, Twater, TgK, Tstart, alb_grid, emis_grid, TgK_wall, Tstart_wall, TmaxLST, + # TmaxLST_wall, first, second, svfalfa, svfbuveg, firstdaytime, timeadd, timestepdec, Tgmap1, + # Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, CI, TgOut1, diffsh, ani): # This is the core function of the SOLWEIG model # 2016-Aug-28 @@ -93,47 +178,47 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s # amaxvalue = max height of buildings # bush = grid representing bushes # Twater = temperature of water (daily) - # TgK, Tstart, TgK_wall, Tstart_wall, TmaxLST,TmaxLST_wall, + # TgK, Tstart, TgK_wall, Tstart_wall, TmaxLST,TmaxLST_wall, # alb_grid, emis_grid = albedo and emmissivity on ground # first, second = conneted to old Ts model (source area based on Smidt et al.) # svfalfa = SVF recalculated to angle - # svfbuveg = complete SVF - # firstdaytime, timeadd, timestepdec, Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, + # svfbuveg = complete SVF + # firstdaytime, timeadd, timestepdec, Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, # CI = Clearness index # TgOut1 = old Ts model # diffsh, ani = Used in anisotrpic models (Wallenberg et al. 2019, 2022) # # # Core program start # # # # Instrument offset in degrees - t = 0. + t = 0.0 # Stefan Bolzmans Constant SBC = 5.67051e-8 # Degrees to radians - deg2rad = np.pi/180 + deg2rad = np.pi / 180 # Find sunrise decimal hour - new from 2014a - _, _, _, SNUP = daylen(jday, location['latitude']) + _, _, _, SNUP = daylen(jday, location["latitude"]) # Vapor pressure - ea = 6.107 * 10 ** ((7.5 * Ta) / (237.3 + Ta)) * (RH / 100.) + ea = 6.107 * 10 ** ((7.5 * Ta) / (237.3 + Ta)) * (RH / 100.0) # Determination of clear - sky emissivity from Prata (1996) msteg = 46.5 * (ea / (Ta + 273.15)) esky = (1 - (1 + msteg) * np.exp(-((1.2 + 3.0 * msteg) ** 0.5))) + elvis # -0.04 old error from Jonsson et al.2006 - if altitude > 0: # # # # # # DAYTIME # # # # # # + if altitude > 0: # # # # # # DAYTIME # # # # # # # Clearness Index on Earth's surface after Crawford and Dunchon (1999) with a correction # factor for low sun elevations after Lindberg et al.(2008) - I0, CI, Kt, I0et, CIuncorr = clearnessindex_2013b(zen, jday, Ta, RH / 100., radG, location, P) - if (CI > 1) or (CI == np.inf): + I0, CI, Kt, I0et, CIuncorr = clearnessindex_2013b(zen, jday, Ta, RH / 100.0, radG, location, P) + if (CI > 1) or (np.inf == CI): CI = 1 # Estimation of radD and radI if not measured after Reindl et al.(1990) if onlyglobal == 1: - I0, CI, Kt, I0et, CIuncorr = clearnessindex_2013b(zen, jday, Ta, RH / 100., radG, location, P) - if (CI > 1) or (CI == np.inf): + I0, CI, Kt, I0et, CIuncorr = clearnessindex_2013b(zen, jday, Ta, RH / 100.0, radG, location, P) + if (CI > 1) or (np.inf == CI): CI = 1 radI, radD = diffusefraction(radG, altitude, Kt, Ta, RH) @@ -142,15 +227,15 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s # Anisotropic Diffuse Radiation after Perez et al. 1993 if anisotropic_sky == 1: patchchoice = 1 - zenDeg = zen*(180/np.pi) + zenDeg = zen * (180 / np.pi) # Relative luminance - lv, pc_, pb_ = Perez_v3(zenDeg, azimuth, radD, radI, jday, patchchoice, patch_option) + lv, pc_, pb_ = Perez_v3(zenDeg, azimuth, radD, radI, jday, patchchoice, patch_option) # Total relative luminance from sky, i.e. from each patch, into each cell - aniLum = np.zeros((rows, cols)) + aniLum = np.zeros((rows, cols), dtype=np.float32) for idx in range(lv.shape[0]): - aniLum += diffsh[:,:,idx] * lv[idx,2] + aniLum += diffsh[:, :, idx] * lv[idx, 2] - dRad = aniLum * radD # Total diffuse radiation from sky into each cell + dRad = aniLum * radD # Total diffuse radiation from sky into each cell else: dRad = radD * svfbuveg patchchoice = 1 @@ -158,29 +243,46 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s # Shadow images if usevegdem == 1: - vegsh, sh, _, wallsh, wallsun, wallshve, _, facesun, wallsh_ = shadowingfunction_wallheight_23(dsm, vegdem, vegdem2, - azimuth, altitude, scale, amaxvalue, bush, walls, dirwalls * np.pi / 180., walls_scheme, dirwalls_scheme * np.pi/180.) + vegsh, sh, _, wallsh, wallsun, wallshve, _, facesun, wallsh_ = shadowingfunction_wallheight_23( + dsm, + vegdem, + vegdem2, + azimuth, + altitude, + scale, + amaxvalue, + bush, + walls, + dirwalls * np.pi / 180.0, + walls_scheme, + dirwalls_scheme * np.pi / 180.0, + ) shadow = sh - (1 - vegsh) * (1 - psi) else: - sh, wallsh, wallsun, facesh, facesun = shadowingfunction_wallheight_13(dsm, azimuth, altitude, scale, - walls, dirwalls * np.pi / 180.) + sh, wallsh, wallsun, facesh, facesun = shadowingfunction_wallheight_13( + dsm, azimuth, altitude, scale, walls, dirwalls * np.pi / 180.0 + ) shadow = sh # # # Surface temperature parameterisation during daytime # # # # # new using max sun alt.instead of dfm # Tgamp = (TgK * altmax - Tstart) + Tstart # Old - Tgamp = TgK * altmax + Tstart # Fixed 2021 + Tgamp = TgK * altmax + Tstart # Fixed 2021 # Tgampwall = (TgK_wall * altmax - (Tstart_wall)) + (Tstart_wall) # Old Tgampwall = TgK_wall * altmax + Tstart_wall - Tg = Tgamp * np.sin((((dectime - np.floor(dectime)) - SNUP / 24) / (TmaxLST / 24 - SNUP / 24)) * np.pi / 2) # 2015 a, based on max sun altitude - Tgwall = Tgampwall * np.sin((((dectime - np.floor(dectime)) - SNUP / 24) / (TmaxLST_wall / 24 - SNUP / 24)) * np.pi / 2) # 2015a, based on max sun altitude + Tg = Tgamp * np.sin( + (((dectime - np.floor(dectime)) - SNUP / 24) / (TmaxLST / 24 - SNUP / 24)) * np.pi / 2 + ) # 2015 a, based on max sun altitude + Tgwall = Tgampwall * np.sin( + (((dectime - np.floor(dectime)) - SNUP / 24) / (TmaxLST_wall / 24 - SNUP / 24)) * np.pi / 2 + ) # 2015a, based on max sun altitude if Tgwall < 0: # temporary for removing low Tg during morning 20130205 # Tg = 0 Tgwall = 0 # New estimation of Tg reduction for non - clear situation based on Reindl et al.1990 - radI0, _ = diffusefraction(I0, altitude, 1., Ta, RH) + radI0, _ = diffusefraction(I0, altitude, 1.0, Ta, RH) corr = 0.1473 * np.log(90 - (zen / np.pi * 180)) + 0.3454 # 20070329 correction of lat, Lindberg et al. 2008 CI_Tg = (radG / radI0) + (1 - corr) if (CI_Tg > 1) or (CI_Tg == np.inf): @@ -190,7 +292,7 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s CI_TgG = (radG / radG0) + (1 - corr) if (CI_TgG > 1) or (CI_TgG == np.inf): CI_TgG = 1 - + # Tg = Tg * CI_Tg # new estimation # Tgwall = Tgwall * CI_Tg Tg = Tg * CI_TgG # new estimation @@ -199,10 +301,47 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s Tg[Tg < 0] = 0 # temporary for removing low Tg during morning 20130205 # # # # Ground View Factors # # # # - gvfLup, gvfalb, gvfalbnosh, gvfLupE, gvfalbE, gvfalbnoshE, gvfLupS, gvfalbS, gvfalbnoshS, gvfLupW, gvfalbW,\ - gvfalbnoshW, gvfLupN, gvfalbN, gvfalbnoshN, gvfSum, gvfNorm = gvf_2018a(wallsun, walls, buildings, scale, shadow, first, - second, dirwalls, Tg, Tgwall, Ta, emis_grid, ewall, alb_grid, SBC, albedo_b, rows, cols, - Twater, lc_grid, landcover) + ( + gvfLup, + gvfalb, + gvfalbnosh, + gvfLupE, + gvfalbE, + gvfalbnoshE, + gvfLupS, + gvfalbS, + gvfalbnoshS, + gvfLupW, + gvfalbW, + gvfalbnoshW, + gvfLupN, + gvfalbN, + gvfalbnoshN, + gvfSum, + gvfNorm, + ) = gvf_2018a( + wallsun, + walls, + buildings, + scale, + shadow, + first, + second, + dirwalls, + Tg, + Tgwall, + Ta, + emis_grid, + ewall, + alb_grid, + SBC, + albedo_b, + rows, + cols, + Twater, + lc_grid, + landcover, + ) # # # # Lup, daytime # # # # # Surface temperature wave delay - new as from 2014a @@ -211,50 +350,102 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s LupS, timeaddnotused, Tgmap1S = TsWaveDelay_2015a(gvfLupS, firstdaytime, timeadd, timestepdec, Tgmap1S) LupW, timeaddnotused, Tgmap1W = TsWaveDelay_2015a(gvfLupW, firstdaytime, timeadd, timestepdec, Tgmap1W) LupN, timeaddnotused, Tgmap1N = TsWaveDelay_2015a(gvfLupN, firstdaytime, timeadd, timestepdec, Tgmap1N) - + # # For Tg output in POIs TgTemp = Tg * shadow + Ta - TgOut, timeadd, TgOut1 = TsWaveDelay_2015a(TgTemp, firstdaytime, timeadd, timestepdec, TgOut1) #timeadd only here v2021a + TgOut, timeadd, TgOut1 = TsWaveDelay_2015a( + TgTemp, firstdaytime, timeadd, timestepdec, TgOut1 + ) # timeadd only here v2021a # Building height angle from svf F_sh = cylindric_wedge(zen, svfalfa, rows, cols) # Fraction shadow on building walls based on sun alt and svf F_sh[np.isnan(F_sh)] = 0.5 # # # # # # # Calculation of shortwave daytime radiative fluxes # # # # # # # - Kdown = radI * shadow * np.sin(altitude * (np.pi / 180)) + dRad + albedo_b * (1 - svfbuveg) * \ - (radG * (1 - F_sh) + radD * F_sh) # *sin(altitude(i) * (pi / 180)) + Kdown = ( + radI * shadow * np.sin(altitude * (np.pi / 180)) + + dRad + + albedo_b * (1 - svfbuveg) * (radG * (1 - F_sh) + radD * F_sh) + ) # *sin(altitude(i) * (pi / 180)) + + Kup, KupE, KupS, KupW, KupN = Kup_veg_2015a( + radI, + radD, + radG, + altitude, + svfbuveg, + albedo_b, + F_sh, + gvfalb, + gvfalbE, + gvfalbS, + gvfalbW, + gvfalbN, + gvfalbnosh, + gvfalbnoshE, + gvfalbnoshS, + gvfalbnoshW, + gvfalbnoshN, + ) + + Keast, Ksouth, Kwest, Knorth, KsideI, KsideD, Kside = Kside_veg_v2022a( + radI, + radD, + radG, + shadow, + svfS, + svfW, + svfN, + svfE, + svfEveg, + svfSveg, + svfWveg, + svfNveg, + azimuth, + altitude, + psi, + t, + albedo_b, + F_sh, + KupE, + KupS, + KupW, + KupN, + cyl, + lv, + anisotropic_sky, + diffsh, + rows, + cols, + asvf, + shmat, + vegshmat, + vbshvegshmat, + ) - Kup, KupE, KupS, KupW, KupN = Kup_veg_2015a(radI, radD, radG, altitude, svfbuveg, albedo_b, F_sh, gvfalb, - gvfalbE, gvfalbS, gvfalbW, gvfalbN, gvfalbnosh, gvfalbnoshE, gvfalbnoshS, gvfalbnoshW, gvfalbnoshN) - - Keast, Ksouth, Kwest, Knorth, KsideI, KsideD, Kside = Kside_veg_v2022a(radI, radD, radG, shadow, svfS, svfW, svfN, svfE, - svfEveg, svfSveg, svfWveg, svfNveg, azimuth, altitude, psi, t, albedo_b, F_sh, KupE, KupS, KupW, - KupN, cyl, lv, anisotropic_sky, diffsh, rows, cols, asvf, shmat, vegshmat, vbshvegshmat) - firstdaytime = 0 else: # # # # # # # NIGHTTIME # # # # # # # # - Tgwall = 0 # CI_Tg = -999 # F_sh = [] # Nocturnal K fluxes set to 0 - Knight = np.zeros((rows, cols)) - Kdown = np.zeros((rows, cols)) - Kwest = np.zeros((rows, cols)) - Kup = np.zeros((rows, cols)) - Keast = np.zeros((rows, cols)) - Ksouth = np.zeros((rows, cols)) - Knorth = np.zeros((rows, cols)) - KsideI = np.zeros((rows, cols)) - KsideD = np.zeros((rows, cols)) - F_sh = np.zeros((rows, cols)) - Tg = np.zeros((rows, cols)) - shadow = np.zeros((rows, cols)) + Knight = np.zeros((rows, cols), dtype=np.float32) + Kdown = np.zeros((rows, cols), dtype=np.float32) + Kwest = np.zeros((rows, cols), dtype=np.float32) + Kup = np.zeros((rows, cols), dtype=np.float32) + Keast = np.zeros((rows, cols), dtype=np.float32) + Ksouth = np.zeros((rows, cols), dtype=np.float32) + Knorth = np.zeros((rows, cols), dtype=np.float32) + KsideI = np.zeros((rows, cols), dtype=np.float32) + KsideD = np.zeros((rows, cols), dtype=np.float32) + F_sh = np.zeros((rows, cols), dtype=np.float32) + Tg = np.zeros((rows, cols), dtype=np.float32) + shadow = np.zeros((rows, cols), dtype=np.float32) CI_Tg = deepcopy(CI) CI_TgG = deepcopy(CI) - dRad = np.zeros((rows,cols)) - Kside = np.zeros((rows,cols)) + dRad = np.zeros((rows, cols), dtype=np.float32) + Kside = np.zeros((rows, cols), dtype=np.float32) # # # # Lup # # # # Lup = SBC * emis_grid * ((Knight + Ta + Tg + 273.15) ** 4) @@ -274,35 +465,70 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s firstdaytime = 1 # # # # Ldown # # # # - Ldown = (svf + svfveg - 1) * esky * SBC * ((Ta + 273.15) ** 4) + (2 - svfveg - svfaveg) * ewall * SBC * \ - ((Ta + 273.15) ** 4) + (svfaveg - svf) * ewall * SBC * ((Ta + 273.15 + Tgwall) ** 4) + \ - (2 - svf - svfveg) * (1 - ewall) * esky * SBC * ((Ta + 273.15) ** 4) # Jonsson et al.(2006) + Ldown = ( + (svf + svfveg - 1) * esky * SBC * ((Ta + 273.15) ** 4) + + (2 - svfveg - svfaveg) * ewall * SBC * ((Ta + 273.15) ** 4) + + (svfaveg - svf) * ewall * SBC * ((Ta + 273.15 + Tgwall) ** 4) + + (2 - svf - svfveg) * (1 - ewall) * esky * SBC * ((Ta + 273.15) ** 4) + ) # Jonsson et al.(2006) # Ldown = Ldown - 25 # Shown by Jonsson et al.(2006) and Duarte et al.(2006) if CI < 0.95: # non - clear conditions c = 1 - CI - Ldown = Ldown * (1 - c) + c * ((svf + svfveg - 1) * SBC * ((Ta + 273.15) ** 4) + (2 - svfveg - svfaveg) * - ewall * SBC * ((Ta + 273.15) ** 4) + (svfaveg - svf) * ewall * SBC * ((Ta + 273.15 + Tgwall) ** 4) + - (2 - svf - svfveg) * (1 - ewall) * SBC * ((Ta + 273.15) ** 4)) # NOT REALLY TESTED!!! BUT MORE CORRECT? + Ldown = Ldown * (1 - c) + c * ( + (svf + svfveg - 1) * SBC * ((Ta + 273.15) ** 4) + + (2 - svfveg - svfaveg) * ewall * SBC * ((Ta + 273.15) ** 4) + + (svfaveg - svf) * ewall * SBC * ((Ta + 273.15 + Tgwall) ** 4) + + (2 - svf - svfveg) * (1 - ewall) * SBC * ((Ta + 273.15) ** 4) + ) # NOT REALLY TESTED!!! BUT MORE CORRECT? # # # # Lside # # # # - Least, Lsouth, Lwest, Lnorth = Lside_veg_v2022a(svfS, svfW, svfN, svfE, svfEveg, svfSveg, svfWveg, svfNveg, - svfEaveg, svfSaveg, svfWaveg, svfNaveg, azimuth, altitude, Ta, Tgwall, SBC, ewall, Ldown, - esky, t, F_sh, CI, LupE, LupS, LupW, LupN, anisotropic_sky) + Least, Lsouth, Lwest, Lnorth = Lside_veg_v2022a( + svfS, + svfW, + svfN, + svfE, + svfEveg, + svfSveg, + svfWveg, + svfNveg, + svfEaveg, + svfSaveg, + svfWaveg, + svfNaveg, + azimuth, + altitude, + Ta, + Tgwall, + SBC, + ewall, + Ldown, + esky, + t, + F_sh, + CI, + LupE, + LupS, + LupW, + LupN, + anisotropic_sky, + ) # New parameterization scheme for wall temperatures if wallScheme == 1: # albedo_g = 0.15 #TODO Change to correct if altitude < 0: wallsh_ = 0 - voxelTable = wall_surface_temperature(voxelTable, wallsh_, altitude, azimuth, timeStep, radI, radD, radG, Ldown, Lup, Ta, esky) + voxelTable = wall_surface_temperature( + voxelTable, wallsh_, altitude, azimuth, timeStep, radI, radD, radG, Ldown, Lup, Ta, esky + ) # Anisotropic sky - if (anisotropic_sky == 1): - if 'lv' not in locals(): + if anisotropic_sky == 1: + if "lv" not in locals(): # Creating skyvault of patches of constant radians (Tregeneza and Sharples, 1993) skyvaultalt, skyvaultazi, _, _, _, _, _ = create_patches(patch_option) - patch_emissivities = np.zeros(skyvaultalt.shape[0]) + patch_emissivities = np.zeros(skyvaultalt.shape[0], dtype=np.float32) x = np.transpose(np.atleast_2d(skyvaultalt)) y = np.transpose(np.atleast_2d(skyvaultazi)) @@ -320,19 +546,72 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s # Create lv from L_patches if nighttime, i.e. lv does not exist if altitude < 0: # CI = deepcopy(CI) - lv = deepcopy(L_patches); KupE = 0; KupS = 0; KupW = 0; KupN = 0 + lv = deepcopy(L_patches) + KupE = 0 + KupS = 0 + KupW = 0 + KupN = 0 # Adjust sky emissivity under semi-cloudy/hazy/cloudy/overcast conditions, i.e. CI lower than 0.95 if CI < 0.95: - esky_c = CI * esky + (1 - CI) * 1. + esky_c = CI * esky + (1 - CI) * 1.0 esky = esky_c - Ldown, Lside, Lside_sky, Lside_veg, Lside_sh, Lside_sun, Lside_ref, Least_, Lwest_, Lnorth_, Lsouth_, \ - Keast, Ksouth, Kwest, Knorth, KsideI, KsideD, Kside, steradians, skyalt = ani_sky(shmat, vegshmat, vbshvegshmat, altitude, azimuth, asvf, cyl, esky, - L_patches, wallScheme, voxelTable, voxelMaps, steradians, Ta, Tgwall, ewall, Lup, radI, radD, radG, lv, - albedo_b, 0, diffsh, shadow, KupE, KupS, KupW, KupN, i) + ( + Ldown, + Lside, + Lside_sky, + Lside_veg, + Lside_sh, + Lside_sun, + Lside_ref, + Least_, + Lwest_, + Lnorth_, + Lsouth_, + Keast, + Ksouth, + Kwest, + Knorth, + KsideI, + KsideD, + Kside, + steradians, + skyalt, + ) = ani_sky( + shmat, + vegshmat, + vbshvegshmat, + altitude, + azimuth, + asvf, + cyl, + esky, + L_patches, + wallScheme, + voxelTable, + voxelMaps, + steradians, + Ta, + Tgwall, + ewall, + Lup, + radI, + radD, + radG, + lv, + albedo_b, + 0, + diffsh, + shadow, + KupE, + KupS, + KupW, + KupN, + i, + ) else: - Lside = np.zeros((rows, cols)) + Lside = np.zeros((rows, cols), dtype=np.float32) L_patches = None # Box and anisotropic longwave @@ -344,20 +623,23 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s # # # # Calculation of radiant flux density and Tmrt # # # # # Human body considered as a cylinder with isotropic all-sky diffuse - if cyl == 1 and anisotropic_sky == 0: - Sstr = absK * (KsideI * Fcyl + (Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside) + absL * \ - ((Ldown + Lup) * Fup + (Lnorth + Least + Lsouth + Lwest) * Fside) - # Human body considered as a cylinder with Perez et al. (1993) (anisotropic sky diffuse) + if cyl == 1 and anisotropic_sky == 0: + Sstr = absK * (KsideI * Fcyl + (Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside) + absL * ( + (Ldown + Lup) * Fup + (Lnorth + Least + Lsouth + Lwest) * Fside + ) + # Human body considered as a cylinder with Perez et al. (1993) (anisotropic sky diffuse) # and Martin and Berdahl (1984) (anisotropic sky longwave) elif cyl == 1 and anisotropic_sky == 1: - Sstr = absK * (Kside * Fcyl + (Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside) + absL * \ - ((Ldown + Lup) * Fup + Lside * Fcyl + (Lnorth + Least + Lsouth + Lwest) * Fside) + Sstr = absK * (Kside * Fcyl + (Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside) + absL * ( + (Ldown + Lup) * Fup + Lside * Fcyl + (Lnorth + Least + Lsouth + Lwest) * Fside + ) # Knorth = nan Ksouth = nan Kwest = nan Keast = nan - else: # Human body considered as a standing cube - Sstr = absK * ((Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside) + absL * \ - ((Ldown + Lup) * Fup + (Lnorth + Least + Lsouth + Lwest) * Fside) + else: # Human body considered as a standing cube + Sstr = absK * ((Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside) + absL * ( + (Ldown + Lup) * Fup + (Lnorth + Least + Lsouth + Lwest) * Fside + ) - Tmrt = np.sqrt(np.sqrt((Sstr / (absL * SBC)))) - 273.2 + Tmrt = np.sqrt(np.sqrt(Sstr / (absL * SBC))) - 273.2 # Add longwave to cardinal directions for output in POI if (cyl == 1) and (anisotropic_sky == 1): @@ -366,7 +648,46 @@ def Solweig_2025a_calc(i, dsm, scale, rows, cols, svf, svfN, svfW, svfE, svfS, s Lnorth += Lnorth_ Lsouth += Lsouth_ - return Tmrt, Kdown, Kup, Ldown, Lup, Tg, ea, esky, I0, CI, shadow, firstdaytime, timestepdec, \ - timeadd, Tgmap1, Tgmap1E, Tgmap1S, Tgmap1W, Tgmap1N, Keast, Ksouth, Kwest, Knorth, Least, \ - Lsouth, Lwest, Lnorth, KsideI, TgOut1, TgOut, radI, radD, \ - Lside, L_patches, CI_Tg, CI_TgG, KsideD, dRad, Kside, steradians, voxelTable + return ( + Tmrt, + Kdown, + Kup, + Ldown, + Lup, + Tg, + ea, + esky, + I0, + CI, + shadow, + firstdaytime, + timestepdec, + timeadd, + Tgmap1, + Tgmap1E, + Tgmap1S, + Tgmap1W, + Tgmap1N, + Keast, + Ksouth, + Kwest, + Knorth, + Least, + Lsouth, + Lwest, + Lnorth, + KsideI, + TgOut1, + TgOut, + radI, + radD, + Lside, + L_patches, + CI_Tg, + CI_TgG, + KsideD, + dRad, + Kside, + steradians, + voxelTable, + ) diff --git a/umep/functions/SOLWEIGpython/Solweig_run.py b/umep/functions/SOLWEIGpython/Solweig_run.py index 4b591aa..a1fdf5a 100644 --- a/umep/functions/SOLWEIGpython/Solweig_run.py +++ b/umep/functions/SOLWEIGpython/Solweig_run.py @@ -62,7 +62,9 @@ def solweig_run(configPath, feedback): standAlone = int(configDict["standalone"]) # Load DSM - dsm_arr, dsm_trf_arr, dsm_crs_wkt, dsm_nd_val = common.load_raster(configDict["filepath_dsm"], bbox=None) + dsm_arr, dsm_trf_arr, dsm_crs_wkt, dsm_nd_val = common.load_raster( + configDict["filepath_dsm"], bbox=None, coerce_f64_to_f32=True + ) # trf is a list: [top left x, w-e pixel size, rotation, top left y, rotation, n-s pixel size] scale = 1 / dsm_trf_arr[1] # pixel resolution in metres left_x = dsm_trf_arr[0] @@ -88,9 +90,9 @@ def solweig_run(configPath, feedback): trunkratio = param["Tree_settings"]["Value"]["Trunk_ratio"] usevegdem = int(configDict["usevegdem"]) if usevegdem == 1: - vegdsm, _, _, _ = common.load_raster(configDict["filepath_cdsm"], bbox=None) + vegdsm, _, _, _ = common.load_raster(configDict["filepath_cdsm"], bbox=None, coerce_f64_to_f32=True) if configDict["filepath_tdsm"] != "": - vegdsm2, _, _, _ = common.load_raster(configDict["filepath_tdsm"], bbox=None) + vegdsm2, _, _, _ = common.load_raster(configDict["filepath_tdsm"], bbox=None, coerce_f64_to_f32=True) else: vegdsm2 = vegdsm * trunkratio else: @@ -100,14 +102,14 @@ def solweig_run(configPath, feedback): # Land cover landcover = int(configDict["landcover"]) if landcover == 1: - lcgrid, _, _, _ = common.load_raster(configDict["filepath_lc"], bbox=None) + lcgrid, _, _, _ = common.load_raster(configDict["filepath_lc"], bbox=None, coerce_f64_to_f32=True) else: lcgrid = 0 # DEM for buildings #TODO: fix nodata in standalone demforbuild = int(configDict["demforbuild"]) if demforbuild == 1: - dem, _, _, dem_nd_val = common.load_raster(configDict["filepath_dem"], bbox=None) + dem, _, _, dem_nd_val = common.load_raster(configDict["filepath_dem"], bbox=None, coerce_f64_to_f32=True) # response to issue and #230 dem[dem == dem_nd_val] = 0.0 if dem.min() < 0: @@ -121,24 +123,44 @@ def solweig_run(configPath, feedback): zip.extractall(configDict["working_dir"]) zip.close() - svf, _, _, _ = common.load_raster(configDict["working_dir"] + "/svf.tif", bbox=None) - svfN, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfN.tif", bbox=None) - svfS, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfS.tif", bbox=None) - svfE, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfE.tif", bbox=None) - svfW, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfW.tif", bbox=None) + svf, _, _, _ = common.load_raster(configDict["working_dir"] + "/svf.tif", bbox=None, coerce_f64_to_f32=True) + svfN, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfN.tif", bbox=None, coerce_f64_to_f32=True) + svfS, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfS.tif", bbox=None, coerce_f64_to_f32=True) + svfE, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfE.tif", bbox=None, coerce_f64_to_f32=True) + svfW, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfW.tif", bbox=None, coerce_f64_to_f32=True) if usevegdem == 1: - svfveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfveg.tif", bbox=None) - svfNveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfNveg.tif", bbox=None) - svfSveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfSveg.tif", bbox=None) - svfEveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfEveg.tif", bbox=None) - svfWveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfWveg.tif", bbox=None) - - svfaveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfaveg.tif", bbox=None) - svfNaveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfNaveg.tif", bbox=None) - svfSaveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfSaveg.tif", bbox=None) - svfEaveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfEaveg.tif", bbox=None) - svfWaveg, _, _, _ = common.load_raster(configDict["working_dir"] + "/svfWaveg.tif", bbox=None) + svfveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfNveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfNveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfSveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfSveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfEveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfEveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfWveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfWveg.tif", bbox=None, coerce_f64_to_f32=True + ) + + svfaveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfaveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfNaveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfNaveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfSaveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfSaveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfEaveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfEaveg.tif", bbox=None, coerce_f64_to_f32=True + ) + svfWaveg, _, _, _ = common.load_raster( + configDict["working_dir"] + "/svfWaveg.tif", bbox=None, coerce_f64_to_f32=True + ) else: svfveg = np.ones((rows, cols)) svfNveg = np.ones((rows, cols)) @@ -156,8 +178,8 @@ def solweig_run(configPath, feedback): # %matlab crazyness around 0 svfalfa = np.arcsin(np.exp(np.log(1.0 - tmp) / 2.0)) - wallheight, _, _, _ = common.load_raster(configDict["filepath_wh"], bbox=None) - wallaspect, _, _, _ = common.load_raster(configDict["filepath_wa"], bbox=None) + wallheight, _, _, _ = common.load_raster(configDict["filepath_wh"], bbox=None, coerce_f64_to_f32=True) + wallaspect, _, _, _ = common.load_raster(configDict["filepath_wa"], bbox=None, coerce_f64_to_f32=True) # Metdata headernum = 1 @@ -318,6 +340,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) # Import shadow matrices (Anisotropic sky) @@ -825,6 +848,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) if configDict["outputkup"] == "1": common.save_raster( @@ -833,6 +857,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) if configDict["outputkdown"] == "1": common.save_raster( @@ -841,6 +866,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) if configDict["outputlup"] == "1": common.save_raster( @@ -849,6 +875,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) if configDict["outputldown"] == "1": common.save_raster( @@ -857,6 +884,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) if configDict["outputsh"] == "1": common.save_raster( @@ -865,6 +893,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) if configDict["outputkdiff"] == "1": common.save_raster( @@ -873,6 +902,7 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) # Sky view image of patches @@ -961,4 +991,5 @@ def solweig_run(configPath, feedback): dsm_trf_arr, dsm_crs_wkt, dsm_nd_val, + coerce_f64_to_f32=True, ) diff --git a/umep/functions/SOLWEIGpython/UTCI_calculations.py b/umep/functions/SOLWEIGpython/UTCI_calculations.py index dbc5959..9a853a5 100644 --- a/umep/functions/SOLWEIGpython/UTCI_calculations.py +++ b/umep/functions/SOLWEIGpython/UTCI_calculations.py @@ -1,222 +1,225 @@ import numpy as np -def utci_polynomial(D_Tmrt, Ta, va, Pa): +def utci_polynomial(D_Tmrt, Ta, va, Pa): # calculate 6th order polynomial as approximation - UTCI_approx = Ta + \ - (6.07562052E-01) + \ - (-2.27712343E-02) * Ta + \ - (8.06470249E-04) * Ta * Ta + \ - (-1.54271372E-04) * Ta * Ta * Ta + \ - (-3.24651735E-06) * Ta * Ta * Ta * Ta + \ - (7.32602852E-08) * Ta * Ta * Ta * Ta * Ta + \ - (1.35959073E-09) * Ta * Ta * Ta * Ta * Ta * Ta + \ - (-2.25836520E+00) * va + \ - (8.80326035E-02) * Ta * va + \ - (2.16844454E-03) * Ta * Ta * va + \ - (-1.53347087E-05) * Ta * Ta * Ta * va + \ - (-5.72983704E-07) * Ta * Ta * Ta * Ta * va + \ - (-2.55090145E-09) * Ta * Ta * Ta * Ta * Ta * va + \ - (-7.51269505E-01) * va * va + \ - (-4.08350271E-03) * Ta * va * va + \ - (-5.21670675E-05) * Ta * Ta * va * va + \ - (1.94544667E-06) * Ta * Ta * Ta * va * va + \ - (1.14099531E-08) * Ta * Ta * Ta * Ta * va * va + \ - (1.58137256E-01) * va * va * va + \ - (-6.57263143E-05) * Ta * va * va * va + \ - (2.22697524E-07) * Ta * Ta * va * va * va + \ - (-4.16117031E-08) * Ta * Ta * Ta * va * va * va + \ - (-1.27762753E-02) * va * va * va * va + \ - (9.66891875E-06) * Ta * va * va * va * va + \ - (2.52785852E-09) * Ta * Ta * va * va * va * va + \ - (4.56306672E-04) * va * va * va * va * va + \ - (-1.74202546E-07) * Ta * va * va * va * va * va + \ - (-5.91491269E-06) * va * va * va * va * va * va + \ - (3.98374029E-01) * D_Tmrt + \ - (1.83945314E-04) * Ta * D_Tmrt + \ - (-1.73754510E-04) * Ta * Ta * D_Tmrt + \ - (-7.60781159E-07) * Ta * Ta * Ta * D_Tmrt + \ - (3.77830287E-08) * Ta * Ta * Ta * Ta * D_Tmrt + \ - (5.43079673E-10) * Ta * Ta * Ta * Ta * Ta * D_Tmrt + \ - (-2.00518269E-02) * va * D_Tmrt + \ - (8.92859837E-04) * Ta * va * D_Tmrt + \ - (3.45433048E-06) * Ta * Ta * va * D_Tmrt + \ - (-3.77925774E-07) * Ta * Ta * Ta * va * D_Tmrt + \ - (-1.69699377E-09) * Ta * Ta * Ta * Ta * va * D_Tmrt + \ - (1.69992415E-04) * va * va * D_Tmrt + \ - (-4.99204314E-05) * Ta * va * va * D_Tmrt + \ - (2.47417178E-07) * Ta * Ta * va * va * D_Tmrt + \ - (1.07596466E-08) * Ta * Ta * Ta * va * va * D_Tmrt + \ - (8.49242932E-05) * va * va * va * D_Tmrt + \ - (1.35191328E-06) * Ta * va * va * va * D_Tmrt + \ - (-6.21531254E-09) * Ta * Ta * va * va * va * D_Tmrt + \ - (-4.99410301E-06) * va * va * va * va * D_Tmrt + \ - (-1.89489258E-08) * Ta * va * va * va * va * D_Tmrt + \ - (8.15300114E-08) * va * va * va * va * va * D_Tmrt + \ - (7.55043090E-04) * D_Tmrt * D_Tmrt + \ - (-5.65095215E-05) * Ta * D_Tmrt * D_Tmrt + \ - (-4.52166564E-07) * Ta * Ta * D_Tmrt * D_Tmrt + \ - (2.46688878E-08) * Ta * Ta * Ta * D_Tmrt * D_Tmrt + \ - (2.42674348E-10) * Ta * Ta * Ta * Ta * D_Tmrt * D_Tmrt + \ - (1.54547250E-04) * va * D_Tmrt * D_Tmrt + \ - (5.24110970E-06) * Ta * va * D_Tmrt * D_Tmrt + \ - (-8.75874982E-08) * Ta * Ta * va * D_Tmrt * D_Tmrt + \ - (-1.50743064E-09) * Ta * Ta * Ta * va * D_Tmrt * D_Tmrt + \ - (-1.56236307E-05) * va * va * D_Tmrt * D_Tmrt + \ - (-1.33895614E-07) * Ta * va * va * D_Tmrt * D_Tmrt + \ - (2.49709824E-09) * Ta * Ta * va * va * D_Tmrt * D_Tmrt + \ - (6.51711721E-07) * va * va * va * D_Tmrt * D_Tmrt + \ - (1.94960053E-09) * Ta * va * va * va * D_Tmrt * D_Tmrt + \ - (-1.00361113E-08) * va * va * va * va * D_Tmrt * D_Tmrt + \ - (-1.21206673E-05) * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-2.18203660E-07) * Ta * D_Tmrt * D_Tmrt * D_Tmrt + \ - (7.51269482E-09) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt + \ - (9.79063848E-11) * Ta * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt + \ - (1.25006734E-06) * va * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-1.81584736E-09) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-3.52197671E-10) * Ta * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-3.36514630E-08) * va * va * D_Tmrt * D_Tmrt * D_Tmrt + \ - (1.35908359E-10) * Ta * va * va * D_Tmrt * D_Tmrt * D_Tmrt + \ - (4.17032620E-10) * va * va * va * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-1.30369025E-09) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (4.13908461E-10) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (9.22652254E-12) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-5.08220384E-09) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-2.24730961E-11) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (1.17139133E-10) * va * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (6.62154879E-10) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (4.03863260E-13) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (1.95087203E-12) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (-4.73602469E-12) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + \ - (5.12733497E+00) * Pa + \ - (-3.12788561E-01) * Ta * Pa + \ - (-1.96701861E-02) * Ta * Ta * Pa + \ - (9.99690870E-04) * Ta * Ta * Ta * Pa + \ - (9.51738512E-06) * Ta * Ta * Ta * Ta * Pa + \ - (-4.66426341E-07) * Ta * Ta * Ta * Ta * Ta * Pa + \ - (5.48050612E-01) * va * Pa + \ - (-3.30552823E-03) * Ta * va * Pa + \ - (-1.64119440E-03) * Ta * Ta * va * Pa + \ - (-5.16670694E-06) * Ta * Ta * Ta * va * Pa + \ - (9.52692432E-07) * Ta * Ta * Ta * Ta * va * Pa + \ - (-4.29223622E-02) * va * va * Pa + \ - (5.00845667E-03) * Ta * va * va * Pa + \ - (1.00601257E-06) * Ta * Ta * va * va * Pa + \ - (-1.81748644E-06) * Ta * Ta * Ta * va * va * Pa + \ - (-1.25813502E-03) * va * va * va * Pa + \ - (-1.79330391E-04) * Ta * va * va * va * Pa + \ - (2.34994441E-06) * Ta * Ta * va * va * va * Pa + \ - (1.29735808E-04) * va * va * va * va * Pa + \ - (1.29064870E-06) * Ta * va * va * va * va * Pa + \ - (-2.28558686E-06) * va * va * va * va * va * Pa + \ - (-3.69476348E-02) * D_Tmrt * Pa + \ - (1.62325322E-03) * Ta * D_Tmrt * Pa + \ - (-3.14279680E-05) * Ta * Ta * D_Tmrt * Pa + \ - (2.59835559E-06) * Ta * Ta * Ta * D_Tmrt * Pa + \ - (-4.77136523E-08) * Ta * Ta * Ta * Ta * D_Tmrt * Pa + \ - (8.64203390E-03) * va * D_Tmrt * Pa + \ - (-6.87405181E-04) * Ta * va * D_Tmrt * Pa + \ - (-9.13863872E-06) * Ta * Ta * va * D_Tmrt * Pa + \ - (5.15916806E-07) * Ta * Ta * Ta * va * D_Tmrt * Pa + \ - (-3.59217476E-05) * va * va * D_Tmrt * Pa + \ - (3.28696511E-05) * Ta * va * va * D_Tmrt * Pa + \ - (-7.10542454E-07) * Ta * Ta * va * va * D_Tmrt * Pa + \ - (-1.24382300E-05) * va * va * va * D_Tmrt * Pa + \ - (-7.38584400E-09) * Ta * va * va * va * D_Tmrt * Pa + \ - (2.20609296E-07) * va * va * va * va * D_Tmrt * Pa + \ - (-7.32469180E-04) * D_Tmrt * D_Tmrt * Pa + \ - (-1.87381964E-05) * Ta * D_Tmrt * D_Tmrt * Pa + \ - (4.80925239E-06) * Ta * Ta * D_Tmrt * D_Tmrt * Pa + \ - (-8.75492040E-08) * Ta * Ta * Ta * D_Tmrt * D_Tmrt * Pa + \ - (2.77862930E-05) * va * D_Tmrt * D_Tmrt * Pa + \ - (-5.06004592E-06) * Ta * va * D_Tmrt * D_Tmrt * Pa + \ - (1.14325367E-07) * Ta * Ta * va * D_Tmrt * D_Tmrt * Pa + \ - (2.53016723E-06) * va * va * D_Tmrt * D_Tmrt * Pa + \ - (-1.72857035E-08) * Ta * va * va * D_Tmrt * D_Tmrt * Pa + \ - (-3.95079398E-08) * va * va * va * D_Tmrt * D_Tmrt * Pa + \ - (-3.59413173E-07) * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (7.04388046E-07) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (-1.89309167E-08) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (-4.79768731E-07) * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (7.96079978E-09) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (1.62897058E-09) * va * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (3.94367674E-08) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (-1.18566247E-09) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (3.34678041E-10) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (-1.15606447E-10) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + \ - (-2.80626406E+00) * Pa * Pa + \ - (5.48712484E-01) * Ta * Pa * Pa + \ - (-3.99428410E-03) * Ta * Ta * Pa * Pa + \ - (-9.54009191E-04) * Ta * Ta * Ta * Pa * Pa + \ - (1.93090978E-05) * Ta * Ta * Ta * Ta * Pa * Pa + \ - (-3.08806365E-01) * va * Pa * Pa + \ - (1.16952364E-02) * Ta * va * Pa * Pa + \ - (4.95271903E-04) * Ta * Ta * va * Pa * Pa + \ - (-1.90710882E-05) * Ta * Ta * Ta * va * Pa * Pa + \ - (2.10787756E-03) * va * va * Pa * Pa + \ - (-6.98445738E-04) * Ta * va * va * Pa * Pa + \ - (2.30109073E-05) * Ta * Ta * va * va * Pa * Pa + \ - (4.17856590E-04) * va * va * va * Pa * Pa + \ - (-1.27043871E-05) * Ta * va * va * va * Pa * Pa + \ - (-3.04620472E-06) * va * va * va * va * Pa * Pa + \ - (5.14507424E-02) * D_Tmrt * Pa * Pa + \ - (-4.32510997E-03) * Ta * D_Tmrt * Pa * Pa + \ - (8.99281156E-05) * Ta * Ta * D_Tmrt * Pa * Pa + \ - (-7.14663943E-07) * Ta * Ta * Ta * D_Tmrt * Pa * Pa + \ - (-2.66016305E-04) * va * D_Tmrt * Pa * Pa + \ - (2.63789586E-04) * Ta * va * D_Tmrt * Pa * Pa + \ - (-7.01199003E-06) * Ta * Ta * va * D_Tmrt * Pa * Pa + \ - (-1.06823306E-04) * va * va * D_Tmrt * Pa * Pa + \ - (3.61341136E-06) * Ta * va * va * D_Tmrt * Pa * Pa + \ - (2.29748967E-07) * va * va * va * D_Tmrt * Pa * Pa + \ - (3.04788893E-04) * D_Tmrt * D_Tmrt * Pa * Pa + \ - (-6.42070836E-05) * Ta * D_Tmrt * D_Tmrt * Pa * Pa + \ - (1.16257971E-06) * Ta * Ta * D_Tmrt * D_Tmrt * Pa * Pa + \ - (7.68023384E-06) * va * D_Tmrt * D_Tmrt * Pa * Pa + \ - (-5.47446896E-07) * Ta * va * D_Tmrt * D_Tmrt * Pa * Pa + \ - (-3.59937910E-08) * va * va * D_Tmrt * D_Tmrt * Pa * Pa + \ - (-4.36497725E-06) * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + \ - (1.68737969E-07) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + \ - (2.67489271E-08) * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + \ - (3.23926897E-09) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + \ - (-3.53874123E-02) * Pa * Pa * Pa + \ - (-2.21201190E-01) * Ta * Pa * Pa * Pa + \ - (1.55126038E-02) * Ta * Ta * Pa * Pa * Pa + \ - (-2.63917279E-04) * Ta * Ta * Ta * Pa * Pa * Pa + \ - (4.53433455E-02) * va * Pa * Pa * Pa + \ - (-4.32943862E-03) * Ta * va * Pa * Pa * Pa + \ - (1.45389826E-04) * Ta * Ta * va * Pa * Pa * Pa + \ - (2.17508610E-04) * va * va * Pa * Pa * Pa + \ - (-6.66724702E-05) * Ta * va * va * Pa * Pa * Pa + \ - (3.33217140E-05) * va * va * va * Pa * Pa * Pa + \ - (-2.26921615E-03) * D_Tmrt * Pa * Pa * Pa + \ - (3.80261982E-04) * Ta * D_Tmrt * Pa * Pa * Pa + \ - (-5.45314314E-09) * Ta * Ta * D_Tmrt * Pa * Pa * Pa + \ - (-7.96355448E-04) * va * D_Tmrt * Pa * Pa * Pa + \ - (2.53458034E-05) * Ta * va * D_Tmrt * Pa * Pa * Pa + \ - (-6.31223658E-06) * va * va * D_Tmrt * Pa * Pa * Pa + \ - (3.02122035E-04) * D_Tmrt * D_Tmrt * Pa * Pa * Pa + \ - (-4.77403547E-06) * Ta * D_Tmrt * D_Tmrt * Pa * Pa * Pa + \ - (1.73825715E-06) * va * D_Tmrt * D_Tmrt * Pa * Pa * Pa + \ - (-4.09087898E-07) * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa * Pa + \ - (6.14155345E-01) * Pa * Pa * Pa * Pa + \ - (-6.16755931E-02) * Ta * Pa * Pa * Pa * Pa + \ - (1.33374846E-03) * Ta * Ta * Pa * Pa * Pa * Pa + \ - (3.55375387E-03) * va * Pa * Pa * Pa * Pa + \ - (-5.13027851E-04) * Ta * va * Pa * Pa * Pa * Pa + \ - (1.02449757E-04) * va * va * Pa * Pa * Pa * Pa + \ - (-1.48526421E-03) * D_Tmrt * Pa * Pa * Pa * Pa + \ - (-4.11469183E-05) * Ta * D_Tmrt * Pa * Pa * Pa * Pa + \ - (-6.80434415E-06) * va * D_Tmrt * Pa * Pa * Pa * Pa + \ - (-9.77675906E-06) * D_Tmrt * D_Tmrt * Pa * Pa * Pa * Pa + \ - (8.82773108E-02) * Pa * Pa * Pa * Pa * Pa + \ - (-3.01859306E-03) * Ta * Pa * Pa * Pa * Pa * Pa + \ - (1.04452989E-03) * va * Pa * Pa * Pa * Pa * Pa + \ - (2.47090539E-04) * D_Tmrt * Pa * Pa * Pa * Pa * Pa + \ - (1.48348065E-03) * Pa * Pa * Pa * Pa * Pa * Pa + UTCI_approx = ( + Ta + + (6.07562052e-01) + + (-2.27712343e-02) * Ta + + (8.06470249e-04) * Ta * Ta + + (-1.54271372e-04) * Ta * Ta * Ta + + (-3.24651735e-06) * Ta * Ta * Ta * Ta + + (7.32602852e-08) * Ta * Ta * Ta * Ta * Ta + + (1.35959073e-09) * Ta * Ta * Ta * Ta * Ta * Ta + + (-2.25836520e00) * va + + (8.80326035e-02) * Ta * va + + (2.16844454e-03) * Ta * Ta * va + + (-1.53347087e-05) * Ta * Ta * Ta * va + + (-5.72983704e-07) * Ta * Ta * Ta * Ta * va + + (-2.55090145e-09) * Ta * Ta * Ta * Ta * Ta * va + + (-7.51269505e-01) * va * va + + (-4.08350271e-03) * Ta * va * va + + (-5.21670675e-05) * Ta * Ta * va * va + + (1.94544667e-06) * Ta * Ta * Ta * va * va + + (1.14099531e-08) * Ta * Ta * Ta * Ta * va * va + + (1.58137256e-01) * va * va * va + + (-6.57263143e-05) * Ta * va * va * va + + (2.22697524e-07) * Ta * Ta * va * va * va + + (-4.16117031e-08) * Ta * Ta * Ta * va * va * va + + (-1.27762753e-02) * va * va * va * va + + (9.66891875e-06) * Ta * va * va * va * va + + (2.52785852e-09) * Ta * Ta * va * va * va * va + + (4.56306672e-04) * va * va * va * va * va + + (-1.74202546e-07) * Ta * va * va * va * va * va + + (-5.91491269e-06) * va * va * va * va * va * va + + (3.98374029e-01) * D_Tmrt + + (1.83945314e-04) * Ta * D_Tmrt + + (-1.73754510e-04) * Ta * Ta * D_Tmrt + + (-7.60781159e-07) * Ta * Ta * Ta * D_Tmrt + + (3.77830287e-08) * Ta * Ta * Ta * Ta * D_Tmrt + + (5.43079673e-10) * Ta * Ta * Ta * Ta * Ta * D_Tmrt + + (-2.00518269e-02) * va * D_Tmrt + + (8.92859837e-04) * Ta * va * D_Tmrt + + (3.45433048e-06) * Ta * Ta * va * D_Tmrt + + (-3.77925774e-07) * Ta * Ta * Ta * va * D_Tmrt + + (-1.69699377e-09) * Ta * Ta * Ta * Ta * va * D_Tmrt + + (1.69992415e-04) * va * va * D_Tmrt + + (-4.99204314e-05) * Ta * va * va * D_Tmrt + + (2.47417178e-07) * Ta * Ta * va * va * D_Tmrt + + (1.07596466e-08) * Ta * Ta * Ta * va * va * D_Tmrt + + (8.49242932e-05) * va * va * va * D_Tmrt + + (1.35191328e-06) * Ta * va * va * va * D_Tmrt + + (-6.21531254e-09) * Ta * Ta * va * va * va * D_Tmrt + + (-4.99410301e-06) * va * va * va * va * D_Tmrt + + (-1.89489258e-08) * Ta * va * va * va * va * D_Tmrt + + (8.15300114e-08) * va * va * va * va * va * D_Tmrt + + (7.55043090e-04) * D_Tmrt * D_Tmrt + + (-5.65095215e-05) * Ta * D_Tmrt * D_Tmrt + + (-4.52166564e-07) * Ta * Ta * D_Tmrt * D_Tmrt + + (2.46688878e-08) * Ta * Ta * Ta * D_Tmrt * D_Tmrt + + (2.42674348e-10) * Ta * Ta * Ta * Ta * D_Tmrt * D_Tmrt + + (1.54547250e-04) * va * D_Tmrt * D_Tmrt + + (5.24110970e-06) * Ta * va * D_Tmrt * D_Tmrt + + (-8.75874982e-08) * Ta * Ta * va * D_Tmrt * D_Tmrt + + (-1.50743064e-09) * Ta * Ta * Ta * va * D_Tmrt * D_Tmrt + + (-1.56236307e-05) * va * va * D_Tmrt * D_Tmrt + + (-1.33895614e-07) * Ta * va * va * D_Tmrt * D_Tmrt + + (2.49709824e-09) * Ta * Ta * va * va * D_Tmrt * D_Tmrt + + (6.51711721e-07) * va * va * va * D_Tmrt * D_Tmrt + + (1.94960053e-09) * Ta * va * va * va * D_Tmrt * D_Tmrt + + (-1.00361113e-08) * va * va * va * va * D_Tmrt * D_Tmrt + + (-1.21206673e-05) * D_Tmrt * D_Tmrt * D_Tmrt + + (-2.18203660e-07) * Ta * D_Tmrt * D_Tmrt * D_Tmrt + + (7.51269482e-09) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt + + (9.79063848e-11) * Ta * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt + + (1.25006734e-06) * va * D_Tmrt * D_Tmrt * D_Tmrt + + (-1.81584736e-09) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt + + (-3.52197671e-10) * Ta * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt + + (-3.36514630e-08) * va * va * D_Tmrt * D_Tmrt * D_Tmrt + + (1.35908359e-10) * Ta * va * va * D_Tmrt * D_Tmrt * D_Tmrt + + (4.17032620e-10) * va * va * va * D_Tmrt * D_Tmrt * D_Tmrt + + (-1.30369025e-09) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (4.13908461e-10) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (9.22652254e-12) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (-5.08220384e-09) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (-2.24730961e-11) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (1.17139133e-10) * va * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (6.62154879e-10) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (4.03863260e-13) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (1.95087203e-12) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (-4.73602469e-12) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt + + (5.12733497e00) * Pa + + (-3.12788561e-01) * Ta * Pa + + (-1.96701861e-02) * Ta * Ta * Pa + + (9.99690870e-04) * Ta * Ta * Ta * Pa + + (9.51738512e-06) * Ta * Ta * Ta * Ta * Pa + + (-4.66426341e-07) * Ta * Ta * Ta * Ta * Ta * Pa + + (5.48050612e-01) * va * Pa + + (-3.30552823e-03) * Ta * va * Pa + + (-1.64119440e-03) * Ta * Ta * va * Pa + + (-5.16670694e-06) * Ta * Ta * Ta * va * Pa + + (9.52692432e-07) * Ta * Ta * Ta * Ta * va * Pa + + (-4.29223622e-02) * va * va * Pa + + (5.00845667e-03) * Ta * va * va * Pa + + (1.00601257e-06) * Ta * Ta * va * va * Pa + + (-1.81748644e-06) * Ta * Ta * Ta * va * va * Pa + + (-1.25813502e-03) * va * va * va * Pa + + (-1.79330391e-04) * Ta * va * va * va * Pa + + (2.34994441e-06) * Ta * Ta * va * va * va * Pa + + (1.29735808e-04) * va * va * va * va * Pa + + (1.29064870e-06) * Ta * va * va * va * va * Pa + + (-2.28558686e-06) * va * va * va * va * va * Pa + + (-3.69476348e-02) * D_Tmrt * Pa + + (1.62325322e-03) * Ta * D_Tmrt * Pa + + (-3.14279680e-05) * Ta * Ta * D_Tmrt * Pa + + (2.59835559e-06) * Ta * Ta * Ta * D_Tmrt * Pa + + (-4.77136523e-08) * Ta * Ta * Ta * Ta * D_Tmrt * Pa + + (8.64203390e-03) * va * D_Tmrt * Pa + + (-6.87405181e-04) * Ta * va * D_Tmrt * Pa + + (-9.13863872e-06) * Ta * Ta * va * D_Tmrt * Pa + + (5.15916806e-07) * Ta * Ta * Ta * va * D_Tmrt * Pa + + (-3.59217476e-05) * va * va * D_Tmrt * Pa + + (3.28696511e-05) * Ta * va * va * D_Tmrt * Pa + + (-7.10542454e-07) * Ta * Ta * va * va * D_Tmrt * Pa + + (-1.24382300e-05) * va * va * va * D_Tmrt * Pa + + (-7.38584400e-09) * Ta * va * va * va * D_Tmrt * Pa + + (2.20609296e-07) * va * va * va * va * D_Tmrt * Pa + + (-7.32469180e-04) * D_Tmrt * D_Tmrt * Pa + + (-1.87381964e-05) * Ta * D_Tmrt * D_Tmrt * Pa + + (4.80925239e-06) * Ta * Ta * D_Tmrt * D_Tmrt * Pa + + (-8.75492040e-08) * Ta * Ta * Ta * D_Tmrt * D_Tmrt * Pa + + (2.77862930e-05) * va * D_Tmrt * D_Tmrt * Pa + + (-5.06004592e-06) * Ta * va * D_Tmrt * D_Tmrt * Pa + + (1.14325367e-07) * Ta * Ta * va * D_Tmrt * D_Tmrt * Pa + + (2.53016723e-06) * va * va * D_Tmrt * D_Tmrt * Pa + + (-1.72857035e-08) * Ta * va * va * D_Tmrt * D_Tmrt * Pa + + (-3.95079398e-08) * va * va * va * D_Tmrt * D_Tmrt * Pa + + (-3.59413173e-07) * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (7.04388046e-07) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (-1.89309167e-08) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (-4.79768731e-07) * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (7.96079978e-09) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (1.62897058e-09) * va * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (3.94367674e-08) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (-1.18566247e-09) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (3.34678041e-10) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (-1.15606447e-10) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa + + (-2.80626406e00) * Pa * Pa + + (5.48712484e-01) * Ta * Pa * Pa + + (-3.99428410e-03) * Ta * Ta * Pa * Pa + + (-9.54009191e-04) * Ta * Ta * Ta * Pa * Pa + + (1.93090978e-05) * Ta * Ta * Ta * Ta * Pa * Pa + + (-3.08806365e-01) * va * Pa * Pa + + (1.16952364e-02) * Ta * va * Pa * Pa + + (4.95271903e-04) * Ta * Ta * va * Pa * Pa + + (-1.90710882e-05) * Ta * Ta * Ta * va * Pa * Pa + + (2.10787756e-03) * va * va * Pa * Pa + + (-6.98445738e-04) * Ta * va * va * Pa * Pa + + (2.30109073e-05) * Ta * Ta * va * va * Pa * Pa + + (4.17856590e-04) * va * va * va * Pa * Pa + + (-1.27043871e-05) * Ta * va * va * va * Pa * Pa + + (-3.04620472e-06) * va * va * va * va * Pa * Pa + + (5.14507424e-02) * D_Tmrt * Pa * Pa + + (-4.32510997e-03) * Ta * D_Tmrt * Pa * Pa + + (8.99281156e-05) * Ta * Ta * D_Tmrt * Pa * Pa + + (-7.14663943e-07) * Ta * Ta * Ta * D_Tmrt * Pa * Pa + + (-2.66016305e-04) * va * D_Tmrt * Pa * Pa + + (2.63789586e-04) * Ta * va * D_Tmrt * Pa * Pa + + (-7.01199003e-06) * Ta * Ta * va * D_Tmrt * Pa * Pa + + (-1.06823306e-04) * va * va * D_Tmrt * Pa * Pa + + (3.61341136e-06) * Ta * va * va * D_Tmrt * Pa * Pa + + (2.29748967e-07) * va * va * va * D_Tmrt * Pa * Pa + + (3.04788893e-04) * D_Tmrt * D_Tmrt * Pa * Pa + + (-6.42070836e-05) * Ta * D_Tmrt * D_Tmrt * Pa * Pa + + (1.16257971e-06) * Ta * Ta * D_Tmrt * D_Tmrt * Pa * Pa + + (7.68023384e-06) * va * D_Tmrt * D_Tmrt * Pa * Pa + + (-5.47446896e-07) * Ta * va * D_Tmrt * D_Tmrt * Pa * Pa + + (-3.59937910e-08) * va * va * D_Tmrt * D_Tmrt * Pa * Pa + + (-4.36497725e-06) * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + + (1.68737969e-07) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + + (2.67489271e-08) * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + + (3.23926897e-09) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa + + (-3.53874123e-02) * Pa * Pa * Pa + + (-2.21201190e-01) * Ta * Pa * Pa * Pa + + (1.55126038e-02) * Ta * Ta * Pa * Pa * Pa + + (-2.63917279e-04) * Ta * Ta * Ta * Pa * Pa * Pa + + (4.53433455e-02) * va * Pa * Pa * Pa + + (-4.32943862e-03) * Ta * va * Pa * Pa * Pa + + (1.45389826e-04) * Ta * Ta * va * Pa * Pa * Pa + + (2.17508610e-04) * va * va * Pa * Pa * Pa + + (-6.66724702e-05) * Ta * va * va * Pa * Pa * Pa + + (3.33217140e-05) * va * va * va * Pa * Pa * Pa + + (-2.26921615e-03) * D_Tmrt * Pa * Pa * Pa + + (3.80261982e-04) * Ta * D_Tmrt * Pa * Pa * Pa + + (-5.45314314e-09) * Ta * Ta * D_Tmrt * Pa * Pa * Pa + + (-7.96355448e-04) * va * D_Tmrt * Pa * Pa * Pa + + (2.53458034e-05) * Ta * va * D_Tmrt * Pa * Pa * Pa + + (-6.31223658e-06) * va * va * D_Tmrt * Pa * Pa * Pa + + (3.02122035e-04) * D_Tmrt * D_Tmrt * Pa * Pa * Pa + + (-4.77403547e-06) * Ta * D_Tmrt * D_Tmrt * Pa * Pa * Pa + + (1.73825715e-06) * va * D_Tmrt * D_Tmrt * Pa * Pa * Pa + + (-4.09087898e-07) * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa * Pa + + (6.14155345e-01) * Pa * Pa * Pa * Pa + + (-6.16755931e-02) * Ta * Pa * Pa * Pa * Pa + + (1.33374846e-03) * Ta * Ta * Pa * Pa * Pa * Pa + + (3.55375387e-03) * va * Pa * Pa * Pa * Pa + + (-5.13027851e-04) * Ta * va * Pa * Pa * Pa * Pa + + (1.02449757e-04) * va * va * Pa * Pa * Pa * Pa + + (-1.48526421e-03) * D_Tmrt * Pa * Pa * Pa * Pa + + (-4.11469183e-05) * Ta * D_Tmrt * Pa * Pa * Pa * Pa + + (-6.80434415e-06) * va * D_Tmrt * Pa * Pa * Pa * Pa + + (-9.77675906e-06) * D_Tmrt * D_Tmrt * Pa * Pa * Pa * Pa + + (8.82773108e-02) * Pa * Pa * Pa * Pa * Pa + + (-3.01859306e-03) * Ta * Pa * Pa * Pa * Pa * Pa + + (1.04452989e-03) * va * Pa * Pa * Pa * Pa * Pa + + (2.47090539e-04) * D_Tmrt * Pa * Pa * Pa * Pa * Pa + + (1.48348065e-03) * Pa * Pa * Pa * Pa * Pa * Pa + ) return UTCI_approx + def utci_calculator(Ta, RH, Tmrt, va10m): # Program for calculating UTCI Temperature (UTCI) # released for public use after termination of COST Action 730 @@ -229,17 +232,27 @@ def utci_calculator(Ta, RH, Tmrt, va10m): UTCI_approx = -999 else: # saturation vapour pressure (es) - g = np.array([-2.8365744E3, - 6.028076559E3, 1.954263612E1, - 2.737830188E-2, - 1.6261698E-5, 7.0229056E-10, - 1.8680009E-13, 2.7150305]) + g = np.array( + [ + -2.8365744e3, + -6.028076559e3, + 1.954263612e1, + -2.737830188e-2, + 1.6261698e-5, + 7.0229056e-10, + -1.8680009e-13, + 2.7150305, + ] + ) tk = Ta + 273.15 # ! air temp in K es = g[7] * np.log(tk) for i in range(0, 7): - es = es + g[i] * tk ** (i + 1 - 3.) + es = es + g[i] * tk ** (i + 1 - 3.0) es = np.exp(es) * 0.01 - ehPa = es * RH / 100. + ehPa = es * RH / 100.0 D_Tmrt = Tmrt - Ta Pa = ehPa / 10.0 # use vapour pressure in kPa @@ -250,6 +263,7 @@ def utci_calculator(Ta, RH, Tmrt, va10m): return UTCI_approx + def utci_calculator_grid(Ta, RH, Tmrt, va10m, feedback): # Program for calculating UTCI Temperature (UTCI) # released for public use after termination of COST Action 730 @@ -263,28 +277,39 @@ def utci_calculator_grid(Ta, RH, Tmrt, va10m, feedback): # If Tair or relative humidity is NaN, set all UTCI to -999 if Ta <= -999 or RH <= -999: - UTCI_approx = np.ones((rows, cols)) * -999 + UTCI_approx = np.ones((rows, cols), dtype=np.float32) * -999 # Else, create an empty raster for UTCI and start estimating UTCI else: - UTCI_approx = np.zeros((rows, cols)) + UTCI_approx = np.zeros((rows, cols), dtype=np.float32) # saturation vapour pressure (es) - g = np.array([-2.8365744E3, - 6.028076559E3, 1.954263612E1, - 2.737830188E-2, - 1.6261698E-5, 7.0229056E-10, - 1.8680009E-13, 2.7150305]) + g = np.array( + [ + -2.8365744e3, + -6.028076559e3, + 1.954263612e1, + -2.737830188e-2, + 1.6261698e-5, + 7.0229056e-10, + -1.8680009e-13, + 2.7150305, + ], + dtype=np.float32, + ) tk = Ta + 273.15 # ! air temp in K es = g[7] * np.log(tk) for i in range(0, 7): - es = es + g[i] * tk ** (i + 1 - 3.) + es = es + g[i] * tk ** (i + 1 - 3.0) es = np.exp(es) * 0.01 - ehPa = es * RH / 100. + ehPa = es * RH / 100.0 Pa = ehPa / 10.0 # use vapour pressure in kPa # For progress counter index = 0 - total = 100. / (rows * cols) + total = 100.0 / (rows * cols) # Start UTCI approximation for iy in np.arange(rows): @@ -292,10 +317,9 @@ def utci_calculator_grid(Ta, RH, Tmrt, va10m, feedback): feedback.setProgressText("Calculation cancelled") break for ix in np.arange(cols): - D_Tmrt = Tmrt[iy, ix] - Ta # va = va10m[iy, ix] - + # Set UTCI to NaN (-999) if Tmrt or wind speed is NaN if Tmrt[iy, ix] <= -999 or va10m[iy, ix] <= -999: UTCI_approx[iy, ix] = -9999 @@ -307,4 +331,4 @@ def utci_calculator_grid(Ta, RH, Tmrt, va10m, feedback): index = index + 1 feedback.setProgress(int(index * total)) - return UTCI_approx \ No newline at end of file + return UTCI_approx diff --git a/umep/functions/SOLWEIGpython/solweig_runner.py b/umep/functions/SOLWEIGpython/solweig_runner.py index ab16cb5..af341e6 100644 --- a/umep/functions/SOLWEIGpython/solweig_runner.py +++ b/umep/functions/SOLWEIGpython/solweig_runner.py @@ -159,7 +159,7 @@ def load_epw_weather(self) -> EnvironData: def load_met_weather(self, header_rows: int = 1, delim: str = " ") -> EnvironData: """Load weather data from a MET file.""" met_path_str = str(common.check_path(self.config.met_path)) - met_data = np.loadtxt(met_path_str, skiprows=header_rows, delimiter=delim) + met_data = np.loadtxt(met_path_str, skiprows=header_rows, delimiter=delim, dtype=np.float32) return EnvironData( self.config, self.params, @@ -201,7 +201,7 @@ def hemispheric_image(self): """ n_patches = self.shadow_mats.shmat.shape[2] n_pois = self.poi_pixel_xys.shape[0] - patch_characteristics = np.zeros((n_patches, n_pois)) + patch_characteristics = np.zeros((n_patches, n_pois), dtype=np.float32) # Get POI indices as integer arrays poi_y = self.poi_pixel_xys[:, 2].astype(int) @@ -378,12 +378,12 @@ def run(self) -> None: if np.unique(self.environ_data.DOY).shape[0] > 1: unique_days = np.unique(self.environ_data.DOY) first_unique_day = self.environ_data.DOY[unique_days[0] == self.environ_data.DOY] - I0_array = np.zeros_like(first_unique_day) + I0_array = np.zeros_like(first_unique_day, dtype=np.float32) else: first_unique_day = self.environ_data.DOY.copy() - I0_array = np.zeros_like(self.environ_data.DOY) + I0_array = np.zeros_like(self.environ_data.DOY, dtype=np.float32) # For Tmrt plot - tmrt_agg = np.zeros((self.raster_data.rows, self.raster_data.cols)) + tmrt_agg = np.zeros((self.raster_data.rows, self.raster_data.cols), dtype=np.float32) # Number of iterations num = len(self.environ_data.Ta) # Prepare progress tracking @@ -633,6 +633,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) if self.config.output_kup: common.save_raster( @@ -641,6 +642,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) if self.config.output_kdown: common.save_raster( @@ -649,6 +651,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) if self.config.output_lup: common.save_raster( @@ -657,6 +660,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) if self.config.output_ldown: common.save_raster( @@ -665,6 +669,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) if self.config.output_sh: common.save_raster( @@ -673,6 +678,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) if self.config.output_kdiff: common.save_raster( @@ -681,6 +687,7 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) # Sky view image of patches @@ -771,7 +778,8 @@ def run(self) -> None: self.location["altitude"], self.shadow_mats.patch_option, ] - ] + ], + dtype=np.float32, ) np.savetxt( self.config.output_dir + "/treeplantersettings.txt", @@ -790,4 +798,5 @@ def run(self) -> None: self.raster_data.trf_arr, self.raster_data.crs_wkt, self.raster_data.nd_val, + coerce_f64_to_f32=True, ) diff --git a/umep/functions/SOLWEIGpython/solweig_runner_core.py b/umep/functions/SOLWEIGpython/solweig_runner_core.py index 7fdfaae..b99ec50 100644 --- a/umep/functions/SOLWEIGpython/solweig_runner_core.py +++ b/umep/functions/SOLWEIGpython/solweig_runner_core.py @@ -42,7 +42,7 @@ def load_poi_data(self) -> Tuple[Any, Any]: poi_path_str = str(common.check_path(self.config.poi_path)) pois_gdf = gpd.read_file(poi_path_str) trf = Affine.from_gdal(*self.raster_data.trf_arr) - self.poi_pixel_xys = np.zeros((len(pois_gdf), 3)) - 999 + self.poi_pixel_xys = np.zeros((len(pois_gdf), 3), dtype=np.float32) - 999 self.poi_names = [] for n, (idx, row) in enumerate(pois_gdf.iterrows()): self.poi_names.append(idx) @@ -80,7 +80,7 @@ def load_woi_data(self) -> Tuple[Any, Any]: """Load walls of interest (WOIs) from a file.""" woi_gdf = gpd.read_file(self.config.woi_file) trf = Affine.from_gdal(*self.raster_data.trf_arr) - self.woi_pixel_xys = np.zeros((len(woi_gdf), 3)) - 999 + self.woi_pixel_xys = np.zeros((len(woi_gdf), 3), dtype=np.float32) - 999 self.woi_names = [] for n, (idx, row) in enumerate(woi_gdf.iterrows()): self.woi_names.append(idx) @@ -155,7 +155,7 @@ def load_epw_weather(self) -> EnvironData: "Wind": filtered_df["wind_speed"], "RH": filtered_df["relative_humidity"], "Tair": filtered_df["temp_air"], - "pres": filtered_df["atmospheric_pressure"].astype(float), # Pascal, ensure float + "pres": filtered_df["atmospheric_pressure"].astype(np.float32), # Pascal, ensure float32 "rain": -999, "Kdown": filtered_df["ghi"], "snow": filtered_df["snow_depth"], @@ -182,17 +182,17 @@ def load_epw_weather(self) -> EnvironData: return EnvironData( self.config, self.params, - YYYY=umep_df["iy"].to_numpy(), - DOY=umep_df["id"].to_numpy(), - hours=umep_df["it"].to_numpy(), - minu=umep_df["imin"].to_numpy(), - Ta=umep_df["Tair"].to_numpy(), - RH=umep_df["RH"].to_numpy(), - radG=umep_df["Kdown"].to_numpy(), - radD=umep_df["ldown"].to_numpy(), - radI=umep_df["Kdiff"].to_numpy(), - P=umep_df["pres"].to_numpy() / 100.0, # convert from Pa to hPa, - Ws=umep_df["Wind"].to_numpy(), + YYYY=umep_df["iy"].to_numpy(dtype=np.float32), + DOY=umep_df["id"].to_numpy(dtype=np.float32), + hours=umep_df["it"].to_numpy(dtype=np.float32), + minu=umep_df["imin"].to_numpy(dtype=np.float32), + Ta=umep_df["Tair"].to_numpy(dtype=np.float32), + RH=umep_df["RH"].to_numpy(dtype=np.float32), + radG=umep_df["Kdown"].to_numpy(dtype=np.float32), + radD=umep_df["ldown"].to_numpy(dtype=np.float32), + radI=umep_df["Kdiff"].to_numpy(dtype=np.float32), + P=umep_df["pres"].to_numpy(dtype=np.float32) / 100.0, # convert from Pa to hPa, + Ws=umep_df["Wind"].to_numpy(dtype=np.float32), location=self.location, UTC=self.config.utc, ) diff --git a/umep/functions/dailyshading.py b/umep/functions/dailyshading.py index 0863d89..7da9530 100644 --- a/umep/functions/dailyshading.py +++ b/umep/functions/dailyshading.py @@ -60,7 +60,7 @@ def dailyshading( vegdem2[vegdem2 == dsm] = 0 # Bush separation - bush = np.logical_not((vegdem2 * vegdem)) * vegdem + bush = np.logical_not(vegdem2 * vegdem) * vegdem shtot = np.zeros((dsm_height, dsm_width)) @@ -71,7 +71,7 @@ def dailyshading( alt = np.zeros(itera) azi = np.zeros(itera) - hour = int(0) + hour = 0 index = 0 time = dict() time["UTC"] = UTC @@ -95,13 +95,7 @@ def dailyshading( doy = day_of_year(year, month, day) - ut_time = ( - doy - - 1.0 - + ((hour - dst) / 24.0) - + (minu / (60.0 * 24.0)) - + (0.0 / (60.0 * 60.0 * 24.0)) - ) + ut_time = doy - 1.0 + ((hour - dst) / 24.0) + (minu / (60.0 * 24.0)) + (0.0 / (60.0 * 60.0 * 24.0)) if ut_time < 0: year = year - 1 @@ -131,39 +125,28 @@ def dailyshading( if time["hour"] == 24: time["hour"] = 0 - time_vector = dt.datetime( - year, month, day, time["hour"], time["min"], time["sec"] - ) + time_vector = dt.datetime(year, month, day, time["hour"], time["min"], time["sec"]) timestr = time_vector.strftime("%Y%m%d_%H%M") if alt[i] > 0: if wallshadow == 1: # Include wall shadows (Issue #121) if usevegdem == 1: - vegsh, sh, _, wallsh, _, wallshve, _, _ = ( - shadowingfunction_wallheight_23( - dsm, - vegdem, - vegdem2, - azi[i], - alt[i], - scale, - amaxvalue, - bush, - walls, - dirwalls * np.pi / 180.0, - ) + vegsh, sh, _, wallsh, _, wallshve, _, _ = shadowingfunction_wallheight_23( + dsm, + vegdem, + vegdem2, + azi[i], + alt[i], + scale, + amaxvalue, + bush, + walls, + dirwalls * np.pi / 180.0, ) # create output folders sh = sh - (1 - vegsh) * (1 - psi) if onetime == 0: - filenamewallshve = ( - folder - + "/facade_shdw_veg/facade_shdw_veg_" - + timestr - + "_LST.tif" - ) - common.save_raster( - filenamewallshve, wallshve, dsm_transf, dsm_crs - ) + filenamewallshve = folder + "/facade_shdw_veg/facade_shdw_veg_" + timestr + "_LST.tif" + common.save_raster(filenamewallshve, wallshve, dsm_transf, dsm_crs, coerce_f64_to_f32=True) else: sh, wallsh, _, _, _ = shadowingfunction_wallheight_13( dsm, azi[i], alt[i], scale, walls, dirwalls * np.pi / 180.0 @@ -171,23 +154,14 @@ def dailyshading( # shtot = shtot + sh if onetime == 0: - filename = ( - folder + "/shadow_ground/shadow_ground_" + timestr + "_LST.tif" - ) - common.save_raster(filename, sh, dsm_transf, dsm_crs) - filenamewallsh = ( - folder - + "/facade_shdw_bldgs/facade_shdw_bldgs_" - + timestr - + "_LST.tif" - ) - common.save_raster(filenamewallsh, wallsh, dsm_transf, dsm_crs) + filename = folder + "/shadow_ground/shadow_ground_" + timestr + "_LST.tif" + common.save_raster(filename, sh, dsm_transf, dsm_crs, coerce_f64_to_f32=True) + filenamewallsh = folder + "/facade_shdw_bldgs/facade_shdw_bldgs_" + timestr + "_LST.tif" + common.save_raster(filenamewallsh, wallsh, dsm_transf, dsm_crs, coerce_f64_to_f32=True) else: if usevegdem == 0: - sh = shadow.shadowingfunctionglobalradiation( - dsm, azi[i], alt[i], scale, 0 - ) + sh = shadow.shadowingfunctionglobalradiation(dsm, azi[i], alt[i], scale, 0) # shtot = shtot + sh else: shadowresult = shadow.shadowingfunction_20( @@ -208,7 +182,7 @@ def dailyshading( if onetime == 0: filename = folder + "/Shadow_" + timestr + "_LST.tif" - common.save_raster(filename, sh, dsm_transf, dsm_crs) + common.save_raster(filename, sh, dsm_transf, dsm_crs, coerce_f64_to_f32=True) shtot = shtot + sh index += 1 @@ -217,15 +191,11 @@ def dailyshading( if wallshadow == 1: if onetime == 1: - filenamewallsh = ( - folder + "/facade_shdw_bldgs/facade_shdw_bldgs_" + timestr + "_LST.tif" - ) - common.save_raster(filenamewallsh, wallsh, dsm_transf, dsm_crs) + filenamewallsh = folder + "/facade_shdw_bldgs/facade_shdw_bldgs_" + timestr + "_LST.tif" + common.save_raster(filenamewallsh, wallsh, dsm_transf, dsm_crs, coerce_f64_to_f32=True) if usevegdem == 1: - filenamewallshve = ( - folder + "/facade_shdw_veg/facade_shdw_veg_" + timestr + "_LST.tif" - ) - common.save_raster(filenamewallshve, wallshve, dsm_transf, dsm_crs) + filenamewallshve = folder + "/facade_shdw_veg/facade_shdw_veg_" + timestr + "_LST.tif" + common.save_raster(filenamewallshve, wallshve, dsm_transf, dsm_crs, coerce_f64_to_f32=True) shadowresult = {"shfinal": shfinal, "time_vector": time_vector} diff --git a/umep/shadow_generator_algorithm.py b/umep/shadow_generator_algorithm.py index 244c78c..55e8214 100644 --- a/umep/shadow_generator_algorithm.py +++ b/umep/shadow_generator_algorithm.py @@ -65,7 +65,7 @@ def generate_shadows( trans_veg: float = 3, trunk_zone_ht_perc: float = 0.25, ): - dsm, dsm_transf, dsm_crs, _dsm_nd = common.load_raster(dsm_path, bbox) + dsm, dsm_transf, dsm_crs, _dsm_nd = common.load_raster(dsm_path, bbox, coerce_f64_to_f32=True) dsm_height, dsm_width = dsm.shape # y rows by x cols dsm_scale = 1 / dsm_transf[1] # y is flipped - so return max for lower row @@ -85,7 +85,9 @@ def generate_shadows( if veg_dsm_path is not None: usevegdem = 1 - veg_dsm, veg_dsm_transf, veg_dsm_crs, _veg_dsm_nd = common.load_raster(veg_dsm_path, bbox) + veg_dsm, veg_dsm_transf, veg_dsm_crs, _veg_dsm_nd = common.load_raster( + veg_dsm_path, bbox, coerce_f64_to_f32=True + ) veg_dsm_height, veg_dsm_width = veg_dsm.shape if not (veg_dsm_width == dsm_width) & (veg_dsm_height == dsm_height): raise ValueError("Error in Vegetation Canopy DSM: All rasters must be of same extent and resolution") @@ -102,11 +104,11 @@ def generate_shadows( if wall_aspect_path and wall_ht_path: print("Facade shadow scheme activated") wallsh = 1 - wh_rast, wh_transf, wh_crs, _wh_nd = common.load_raster(wall_ht_path, bbox) + wh_rast, wh_transf, wh_crs, _wh_nd = common.load_raster(wall_ht_path, bbox, coerce_f64_to_f32=True) wh_height, wh_width = wh_rast.shape if not (wh_width == dsm_width) & (wh_height == dsm_height): raise ValueError("Error in Wall height raster: All rasters must be of same extent and resolution") - wa_rast, wa_transf, wa_crs, _wa_nd = common.load_raster(wall_aspect_path, bbox) + wa_rast, wa_transf, wa_crs, _wa_nd = common.load_raster(wall_aspect_path, bbox, coerce_f64_to_f32=True) wa_height, wa_width = wa_rast.shape if not (wa_width == dsm_width) & (wa_height == dsm_height): raise ValueError("Error in Wall aspect raster: All rasters must be of same extent and resolution") @@ -168,4 +170,10 @@ def generate_shadows( ) shfinal = shadowresult["shfinal"] - common.save_raster(out_path_str + "/shadow_composite.tif", shfinal, dsm_transf, dsm_crs) + common.save_raster( + out_path_str + "/shadow_composite.tif", + shfinal, + dsm_transf, + dsm_crs, + coerce_f64_to_f32=True, + ) diff --git a/umep/skyviewfactor_algorithm.py b/umep/skyviewfactor_algorithm.py index c5e12f1..ab499b7 100644 --- a/umep/skyviewfactor_algorithm.py +++ b/umep/skyviewfactor_algorithm.py @@ -74,13 +74,13 @@ def generate_svf( out_path_str = str(out_path) # Open the DSM file - dsm, dsm_trf, dsm_crs, dsm_nd = common.load_raster(dsm_path, bbox) + dsm, dsm_trf, dsm_crs, dsm_nd = common.load_raster(dsm_path, bbox, coerce_f64_to_f32=True) dsm_pix_size = dsm_trf[1] dsm_scale = 1 / dsm_pix_size dem = None if dem_path is not None: - dem, dem_trf, dem_crs, _dem_nd = common.load_raster(dem_path, bbox) + dem, dem_trf, dem_crs, _dem_nd = common.load_raster(dem_path, bbox, coerce_f64_to_f32=True) assert dem.shape == dsm.shape, "Mismatching raster shapes for DSM and DEM." assert np.allclose(dsm_trf, dem_trf), "Mismatching spatial transform for DSM and DEM." assert dem_crs == dsm_crs, "Mismatching CRS for DSM and DEM." @@ -89,7 +89,7 @@ def generate_svf( cdsm = None if cdsm_path is not None: use_cdsm = True - cdsm, cdsm_trf, cdsm_crs, _cdsm_nd = common.load_raster(cdsm_path, bbox) + cdsm, cdsm_trf, cdsm_crs, _cdsm_nd = common.load_raster(cdsm_path, bbox, coerce_f64_to_f32=True) assert cdsm.shape == dsm.shape, "Mismatching raster shapes for DSM and CDSM." assert np.allclose(dsm_trf, cdsm_trf), "Mismatching spatial transform for DSM and CDSM." assert cdsm_crs == dsm_crs, "Mismatching CRS for DSM and CDSM." @@ -118,6 +118,7 @@ def generate_svf( dsm_trf, dsm_crs, dsm_nd, + coerce_f64_to_f32=True, ) if dem is not None: common.save_raster( @@ -126,6 +127,7 @@ def generate_svf( dsm_trf, dsm_crs, dsm_nd, + coerce_f64_to_f32=True, ) if use_cdsm: common.save_raster( @@ -134,6 +136,7 @@ def generate_svf( dsm_trf, dsm_crs, dsm_nd, + coerce_f64_to_f32=True, ) common.save_raster( out_path_str + "/input-tdsm.tif", @@ -141,6 +144,7 @@ def generate_svf( dsm_trf, dsm_crs, dsm_nd, + coerce_f64_to_f32=True, ) # compute @@ -153,11 +157,11 @@ def generate_svf( svfbuN = ret["svfN"] # Save the rasters using rasterio - common.save_raster(out_path_str + "/" + "svf.tif", svfbu, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfE.tif", svfbuE, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfS.tif", svfbuS, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfW.tif", svfbuW, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfN.tif", svfbuN, dsm_trf, dsm_crs) + common.save_raster(out_path_str + "/" + "svf.tif", svfbu, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfE.tif", svfbuE, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfS.tif", svfbuS, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfW.tif", svfbuW, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfN.tif", svfbuN, dsm_trf, dsm_crs, coerce_f64_to_f32=True) # Create or update the ZIP file zip_filepath = out_path_str + "/" + "svfs.zip" @@ -194,16 +198,16 @@ def generate_svf( svfNaveg = ret["svfNaveg"] # Save vegetation rasters - common.save_raster(out_path_str + "/" + "svfveg.tif", svfveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfEveg.tif", svfEveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfSveg.tif", svfSveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfWveg.tif", svfWveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfNveg.tif", svfNveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfaveg.tif", svfaveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfEaveg.tif", svfEaveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfSaveg.tif", svfSaveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfWaveg.tif", svfWaveg, dsm_trf, dsm_crs) - common.save_raster(out_path_str + "/" + "svfNaveg.tif", svfNaveg, dsm_trf, dsm_crs) + common.save_raster(out_path_str + "/" + "svfveg.tif", svfveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfEveg.tif", svfEveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfSveg.tif", svfSveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfWveg.tif", svfWveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfNveg.tif", svfNveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfaveg.tif", svfaveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfEaveg.tif", svfEaveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfSaveg.tif", svfSaveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfWaveg.tif", svfWaveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) + common.save_raster(out_path_str + "/" + "svfNaveg.tif", svfNaveg, dsm_trf, dsm_crs, coerce_f64_to_f32=True) # Add vegetation rasters to the ZIP file with zipfile.ZipFile(zip_filepath, "a") as zippo: @@ -234,7 +238,7 @@ def generate_svf( svftotal = svfbu - (1 - svfveg) * (1 - trans_veg) # Save the final svftotal raster - common.save_raster(out_path_str + "/" + "svf_total.tif", svftotal, dsm_trf, dsm_crs) + common.save_raster(out_path_str + "/" + "svf_total.tif", svftotal, dsm_trf, dsm_crs, coerce_f64_to_f32=True) # Save shadow matrices as compressed npz shmat = ret["shmat"] diff --git a/umep/solweig_algorithm.py b/umep/solweig_algorithm.py index baed3d4..5e761f8 100644 --- a/umep/solweig_algorithm.py +++ b/umep/solweig_algorithm.py @@ -129,7 +129,7 @@ def generate_solweig( Path.mkdir(out_path / "shadows", parents=True, exist_ok=True) Path.mkdir(out_path / "Tmrt", parents=True, exist_ok=True) - dsm, dsm_transf, dsm_crs, _dsm_nd = common.load_raster(dsm_path, bbox) + dsm, dsm_transf, dsm_crs, _dsm_nd = common.load_raster(dsm_path, bbox, coerce_f64_to_f32=True) dsm_trf_affine = Affine.from_gdal(*dsm_transf) dsm_scale = 1 / dsm_transf[1] dsm_height, dsm_width = dsm.shape # y rows by x cols @@ -152,7 +152,9 @@ def generate_solweig( if veg_dsm_path is not None: usevegdem = 1 - veg_dsm, veg_dsm_transf, veg_dsm_crs, _veg_dsm_nd = common.load_raster(veg_dsm_path, bbox) + veg_dsm, veg_dsm_transf, veg_dsm_crs, _veg_dsm_nd = common.load_raster( + veg_dsm_path, bbox, coerce_f64_to_f32=True + ) veg_dsm_height, veg_dsm_width = veg_dsm.shape if not (veg_dsm_width == dsm_width) & (veg_dsm_height == dsm_height): raise ValueError("Error in Vegetation Canopy DSM: All rasters must be of same extent and resolution") @@ -182,23 +184,23 @@ def generate_solweig( zip.extractall(temp_dir) zip.close() - svf, _, _, _ = common.load_raster(temp_dir + "/svf.tif", bbox) - svfN, _, _, _ = common.load_raster(temp_dir + "/svfN.tif", bbox) - svfS, _, _, _ = common.load_raster(temp_dir + "/svfS.tif", bbox) - svfE, _, _, _ = common.load_raster(temp_dir + "/svfE.tif", bbox) - svfW, _, _, _ = common.load_raster(temp_dir + "/svfW.tif", bbox) + svf, _, _, _ = common.load_raster(temp_dir + "/svf.tif", bbox, coerce_f64_to_f32=True) + svfN, _, _, _ = common.load_raster(temp_dir + "/svfN.tif", bbox, coerce_f64_to_f32=True) + svfS, _, _, _ = common.load_raster(temp_dir + "/svfS.tif", bbox, coerce_f64_to_f32=True) + svfE, _, _, _ = common.load_raster(temp_dir + "/svfE.tif", bbox, coerce_f64_to_f32=True) + svfW, _, _, _ = common.load_raster(temp_dir + "/svfW.tif", bbox, coerce_f64_to_f32=True) if usevegdem == 1: - svfveg, _, _, _ = common.load_raster(temp_dir + "/svfveg.tif", bbox) - svfNveg, _, _, _ = common.load_raster(temp_dir + "/svfNveg.tif", bbox) - svfSveg, _, _, _ = common.load_raster(temp_dir + "/svfSveg.tif", bbox) - svfEveg, _, _, _ = common.load_raster(temp_dir + "/svfEveg.tif", bbox) - svfWveg, _, _, _ = common.load_raster(temp_dir + "/svfWveg.tif", bbox) - svfaveg, _, _, _ = common.load_raster(temp_dir + "/svfaveg.tif", bbox) - svfNaveg, _, _, _ = common.load_raster(temp_dir + "/svfNaveg.tif", bbox) - svfSaveg, _, _, _ = common.load_raster(temp_dir + "/svfSaveg.tif", bbox) - svfEaveg, _, _, _ = common.load_raster(temp_dir + "/svfEaveg.tif", bbox) - svfWaveg, _, _, _ = common.load_raster(temp_dir + "/svfWaveg.tif", bbox) + svfveg, _, _, _ = common.load_raster(temp_dir + "/svfveg.tif", bbox, coerce_f64_to_f32=True) + svfNveg, _, _, _ = common.load_raster(temp_dir + "/svfNveg.tif", bbox, coerce_f64_to_f32=True) + svfSveg, _, _, _ = common.load_raster(temp_dir + "/svfSveg.tif", bbox, coerce_f64_to_f32=True) + svfEveg, _, _, _ = common.load_raster(temp_dir + "/svfEveg.tif", bbox, coerce_f64_to_f32=True) + svfWveg, _, _, _ = common.load_raster(temp_dir + "/svfWveg.tif", bbox, coerce_f64_to_f32=True) + svfaveg, _, _, _ = common.load_raster(temp_dir + "/svfaveg.tif", bbox, coerce_f64_to_f32=True) + svfNaveg, _, _, _ = common.load_raster(temp_dir + "/svfNaveg.tif", bbox, coerce_f64_to_f32=True) + svfSaveg, _, _, _ = common.load_raster(temp_dir + "/svfSaveg.tif", bbox, coerce_f64_to_f32=True) + svfEaveg, _, _, _ = common.load_raster(temp_dir + "/svfEaveg.tif", bbox, coerce_f64_to_f32=True) + svfWaveg, _, _, _ = common.load_raster(temp_dir + "/svfWaveg.tif", bbox, coerce_f64_to_f32=True) else: svfveg = np.ones((dsm_height, dsm_width)) svfNveg = np.ones((dsm_height, dsm_width)) @@ -218,11 +220,11 @@ def generate_solweig( tmp[tmp < 0.0] = 0.0 svfalfa = np.arcsin(np.exp(np.log(1.0 - tmp) / 2.0)) - wh_rast, wh_transf, wh_crs, _wh_nd = common.load_raster(wall_ht_path, bbox) + wh_rast, wh_transf, wh_crs, _wh_nd = common.load_raster(wall_ht_path, bbox, coerce_f64_to_f32=True) wh_height, wh_width = wh_rast.shape if not (wh_width == dsm_width) & (wh_height == dsm_height): raise ValueError("Error in Wall height raster: All rasters must be of same extent and resolution") - wa_rast, wa_transf, wa_crs, _wa_nd = common.load_raster(wall_aspect_path, bbox) + wa_rast, wa_transf, wa_crs, _wa_nd = common.load_raster(wall_aspect_path, bbox, coerce_f64_to_f32=True) wa_height, wa_width = wa_rast.shape if not (wa_width == dsm_width) & (wa_height == dsm_height): raise ValueError("Error in Wall aspect raster: All rasters must be of same extent and resolution") @@ -678,10 +680,7 @@ def generate_solweig( poi_results.append(result_row) common.save_raster( - out_path_str + "/Tmrt/Tmrt_" + time_code + ".tif", - Tmrt, - dsm_transf, - dsm_crs, + out_path_str + "/Tmrt/Tmrt_" + time_code + ".tif", Tmrt, dsm_transf, dsm_crs, coerce_f64_to_f32=True ) common.save_raster( out_path_str @@ -699,6 +698,7 @@ def generate_solweig( shadow, dsm_transf, dsm_crs, + coerce_f64_to_f32=True, ) # After the main loop, write all POI results to a single GPKG file with multi-index @@ -735,6 +735,6 @@ def generate_solweig( umep_df.to_csv(out_path_str + "/metforcing.csv") tmrtplot = tmrtplot / Ta.__len__() # fix average Tmrt instead of sum, 20191022 - common.save_raster(out_path_str + "/Tmrt_average.tif", tmrtplot, dsm_transf, dsm_crs) + common.save_raster(out_path_str + "/Tmrt_average.tif", tmrtplot, dsm_transf, dsm_crs, coerce_f64_to_f32=True) rmtree(temp_dir, ignore_errors=True) diff --git a/umep/wall_heightaspect_algorithm.py b/umep/wall_heightaspect_algorithm.py index b9b36aa..7d4b6cb 100644 --- a/umep/wall_heightaspect_algorithm.py +++ b/umep/wall_heightaspect_algorithm.py @@ -58,7 +58,7 @@ def generate_wall_hts( wall_limit: float = 1, ): """ """ - dsm_rast, dsm_transf, dsm_crs, _dsm_nd = common.load_raster(dsm_path, bbox) + dsm_rast, dsm_transf, dsm_crs, _dsm_nd = common.load_raster(dsm_path, bbox, coerce_f64_to_f32=True) dsm_scale = 1 / dsm_transf[1] out_path = Path(out_dir) @@ -66,7 +66,7 @@ def generate_wall_hts( out_path_str = str(out_path) walls = wa.findwalls(dsm_rast, wall_limit) - common.save_raster(out_path_str + "/" + "wall_hts.tif", walls, dsm_transf, dsm_crs) + common.save_raster(out_path_str + "/" + "wall_hts.tif", walls, dsm_transf, dsm_crs, coerce_f64_to_f32=True) dirwalls = wa.filter1Goodwin_as_aspect_v3(walls, dsm_scale, dsm_rast) - common.save_raster(out_path_str + "/" + "wall_aspects.tif", dirwalls, dsm_transf, dsm_crs) + common.save_raster(out_path_str + "/" + "wall_aspects.tif", dirwalls, dsm_transf, dsm_crs, coerce_f64_to_f32=True) From 7676d6ebbfbec244b8970b823bb91007ad461f4c Mon Sep 17 00:00:00 2001 From: Gareth Simons Date: Wed, 19 Nov 2025 11:02:02 +0000 Subject: [PATCH 2/3] tag for tests 0.0.1a5 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index f1299b5..48b8304 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "umep" -version = "0.0.1b47" +version = "0.0.1a5" description = "urban multi-scale environmental predictor" readme = "README.md" requires-python = ">=3.9, <3.14" From a2198d5a5abec5a72940e06dd0fa68de09911ef4 Mon Sep 17 00:00:00 2001 From: Gareth Simons Date: Wed, 19 Nov 2025 21:52:21 +0000 Subject: [PATCH 3/3] enforcing f32 --- pyproject.toml | 2 +- tests/test_solweig_rasters.py | 63 ++++++++++++++ umep/class_configs.py | 150 +++++++++++++++++++++++----------- 3 files changed, 167 insertions(+), 48 deletions(-) create mode 100644 tests/test_solweig_rasters.py diff --git a/pyproject.toml b/pyproject.toml index 48b8304..af6f029 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "umep" -version = "0.0.1a5" +version = "0.0.1a6" description = "urban multi-scale environmental predictor" readme = "README.md" requires-python = ">=3.9, <3.14" diff --git a/tests/test_solweig_rasters.py b/tests/test_solweig_rasters.py new file mode 100644 index 0000000..edb7817 --- /dev/null +++ b/tests/test_solweig_rasters.py @@ -0,0 +1,63 @@ +from pathlib import Path +from typing import Optional + +import numpy as np + +from umep.class_configs import SolweigConfig +from umep.functions.SOLWEIGpython.solweig_runner_core import SolweigRunCore + +PROJECT_ROOT = Path(__file__).resolve().parents[1] +CONFIG_PATH = PROJECT_ROOT / "demos/data/athens/configsolweig.ini" +PARAMS_PATH = PROJECT_ROOT / "demos/data/athens/parametersforsolweig.json" + + +def _prepare_temp_config(tmp_path: Path) -> Path: + """Return a config file that writes outputs into a temp dir.""" + config = SolweigConfig() + config.from_file(str(CONFIG_PATH)) + + output_dir = tmp_path / "output" + working_dir = tmp_path / "working" + + config.output_dir = str(output_dir) + config.working_dir = str(working_dir) + + path_overrides = { + "dsm_path": PROJECT_ROOT / "demos/data/athens/DSM.tif", + "cdsm_path": PROJECT_ROOT / "temp/athens/CDSM.tif", + "dem_path": PROJECT_ROOT / "demos/data/athens/DEM.tif", + "wh_path": PROJECT_ROOT / "temp/athens/walls/wall_hts.tif", + "wa_path": PROJECT_ROOT / "temp/athens/walls/wall_aspects.tif", + "svf_path": PROJECT_ROOT / "temp/athens/svf/svfs.zip", + "aniso_path": PROJECT_ROOT / "temp/athens/svf/shadowmats.npz", + "epw_path": PROJECT_ROOT / "demos/data/athens/athens_2023.epw", + } + for attr, path in path_overrides.items(): + setattr(config, attr, str(path)) + + tmp_config = tmp_path / "configsolweig.ini" + config.to_file(str(tmp_config)) + return tmp_config + + +def _assert_float32(array: Optional[np.ndarray], name: str): + if array is not None: + assert array.dtype == np.float32, f"{name} dtype was {array.dtype}, expected float32" + + +def test_athens_runner_rasters_are_float32(tmp_path): + tmp_config = _prepare_temp_config(tmp_path) + runner = SolweigRunCore(str(tmp_config), str(PARAMS_PATH)) + raster_data = runner.raster_data + + _assert_float32(raster_data.dsm, "dsm") + _assert_float32(raster_data.wallheight, "wallheight") + _assert_float32(raster_data.wallaspect, "wallaspect") + _assert_float32(raster_data.dem, "dem") + _assert_float32(raster_data.cdsm, "cdsm") + _assert_float32(raster_data.tdsm, "tdsm") + _assert_float32(raster_data.bush, "bush") + _assert_float32(raster_data.svfbuveg, "svfbuveg") + _assert_float32(raster_data.buildings, "buildings") + + runner.test_hook() diff --git a/umep/class_configs.py b/umep/class_configs.py index fece9fc..e0d66c5 100644 --- a/umep/class_configs.py +++ b/umep/class_configs.py @@ -2,7 +2,7 @@ import logging import zipfile from dataclasses import dataclass -from typing import Any, Optional, Tuple, Union +from typing import Any, Optional, Union import numpy as np import scipy.ndimage as ndi @@ -259,6 +259,16 @@ def __init__( # These variables lag across days until updated Twater = None CI = 1 + + def _scalar(value: Any) -> float: + """Convert sun-position outputs to native float scalars.""" + if isinstance(value, (int, float, np.floating)): + return float(value) + arr = np.asarray(value) + if arr.ndim == 0: + return float(arr) + return float(arr.reshape(-1)[0]) + # Iterate over time steps and set vars for i in range(data_len): YMD = datetime.datetime(int(self.YYYY[i]), 1, 1) + datetime.timedelta(int(self.DOY[i]) - 1) @@ -267,8 +277,8 @@ def __init__( fifteen = 0.0 sunmaximum = -90.0 sunmax["zenith"] = 90.0 - while sunmaximum <= 90.0 - sunmax["zenith"]: - sunmaximum = 90.0 - sunmax["zenith"] + while sunmaximum <= 90.0 - _scalar(sunmax["zenith"]): + sunmaximum = 90.0 - _scalar(sunmax["zenith"]) fifteen = fifteen + 15.0 / 1440.0 HM = datetime.timedelta(days=(60 * 10) / 1440.0 + fifteen) YMDHM = YMD + HM @@ -278,7 +288,7 @@ def __init__( time["hour"] = YMDHM.hour time["min"] = YMDHM.minute sunmax = sp.sun_position(time, location) - self.altmax[i] = sunmaximum + self.altmax[i] = float(sunmaximum) # Calculate sun position half = datetime.timedelta(days=halftimestepdec) H = datetime.timedelta(hours=int(self.hours[i])) @@ -290,13 +300,15 @@ def __init__( time["hour"] = YMDHM.hour time["min"] = YMDHM.minute sun = sp.sun_position(time, location) - if (sun["zenith"] > 89.0) & ( - sun["zenith"] <= 90.0 + sun_zenith = _scalar(sun["zenith"]) + if (sun_zenith > 89.0) & ( + sun_zenith <= 90.0 ): # Hopefully fixes weird values in Perez et al. when altitude < 1.0, i.e. close to sunrise/sunset - sun["zenith"] = 89.0 - self.altitude[i] = 90.0 - sun["zenith"] - self.zen[i] = sun["zenith"] * (np.pi / 180.0) - self.azimuth[i] = sun["azimuth"] + sun_zenith = 89.0 + sun_azimuth = _scalar(sun["azimuth"]) + self.altitude[i] = 90.0 - sun_zenith + self.zen[i] = sun_zenith * (np.pi / 180.0) + self.azimuth[i] = sun_azimuth # day of year and check for leap year # if calendar.isleap(time["year"]): # dayspermonth = np.atleast_2d([31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) @@ -436,7 +448,9 @@ def __init__(self, model_configs: SolweigConfig): # Calculate SVF alpha tmp = self.svf + self.svf_veg - 1.0 tmp[tmp < 0.0] = 0.0 - self.svfalfa = np.arcsin(np.exp(np.log(1.0 - tmp) / 2.0)) + eps = np.finfo(np.float32).tiny + safe_term = np.clip(1.0 - tmp, eps, 1.0) + self.svfalfa = np.arcsin(np.exp(np.log(safe_term) / 2.0)) logger.info("SVF data loaded and processed.") @@ -450,6 +464,9 @@ def raster_preprocessing( amax_local_window_m: int = 100, amax_local_perc: float = 99.9, ): + nan32 = np.float32(np.nan) + zero32 = np.float32(0.0) + threshold = np.float32(0.1) # amax if dem is None: amaxvalue = float(np.nanmax(dsm) - np.nanmin(dsm)) @@ -473,6 +490,7 @@ def raster_preprocessing( # CDSM is relative to flat surface without DEM if cdsm is None: cdsm = np.zeros_like(dsm) + tdsm = np.zeros_like(dsm) else: if np.nanmax(cdsm) > 50: logger.warning( @@ -501,20 +519,25 @@ def raster_preprocessing( # Set vegetated pixels to DEM + CDSM otherwise DSM + CDSM if dem is not None: - cdsm = np.where(~np.isnan(dem), dem + cdsm, np.nan) - cdsm = np.where(cdsm - dem < 0.1, 0, cdsm) - tdsm = np.where(~np.isnan(dem), dem + tdsm, np.nan) - tdsm = np.where(tdsm - dem < 0.1, 0, tdsm) + cdsm = np.where(~np.isnan(dem), dem + cdsm, nan32) + cdsm = np.where(cdsm - dem < threshold, zero32, cdsm) + tdsm = np.where(~np.isnan(dem), dem + tdsm, nan32) + tdsm = np.where(tdsm - dem < threshold, zero32, tdsm) else: - cdsm = np.where(~np.isnan(dsm), dsm + cdsm, np.nan) - cdsm = np.where(cdsm - dsm < 0.1, 0, cdsm) - tdsm = np.where(~np.isnan(dsm), dsm + tdsm, np.nan) - tdsm = np.where(tdsm - dsm < 0.1, 0, tdsm) + cdsm = np.where(~np.isnan(dsm), dsm + cdsm, nan32) + cdsm = np.where(cdsm - dsm < threshold, zero32, cdsm) + tdsm = np.where(~np.isnan(dsm), dsm + tdsm, nan32) + tdsm = np.where(tdsm - dsm < threshold, zero32, tdsm) logger.info("Calculated max height for shadows: %.2fm", amaxvalue) if amaxvalue > 100: logger.warning("Max shadow height exceeds 100m, double-check the input rasters for anomalies.") + dsm = np.asarray(dsm, dtype=np.float32) + dem = None if dem is None else np.asarray(dem, dtype=np.float32) + cdsm = np.asarray(cdsm, dtype=np.float32) + tdsm = np.asarray(tdsm, dtype=np.float32) + return dsm, dem, cdsm, tdsm, amaxvalue @@ -523,9 +546,9 @@ class RasterData: amaxvalue: float dsm: np.ndarray - crs_wkt: str - trf_arr: np.ndarray - nd_val: float + crs_wkt: Optional[str] + trf_arr: list[float] + nd_val: Optional[float] scale: float rows: int cols: int @@ -547,14 +570,27 @@ def __init__( amax_local_window_m: int = 100, amax_local_perc: float = 99.9, ): + if model_configs.dsm_path is None: + raise ValueError("DSM path must be provided before initializing raster data.") + if model_configs.wh_path is None or model_configs.wa_path is None: + raise ValueError("Wall height and aspect rasters must be provided.") + if model_configs.output_dir is None: + raise ValueError("Output directory must be configured before initializing raster data.") + output_dir = model_configs.output_dir + dsm_path = model_configs.dsm_path + wh_path = model_configs.wh_path + wa_path = model_configs.wa_path # Load DSM self.dsm, self.trf_arr, self.crs_wkt, self.nd_val = common.load_raster( - model_configs.dsm_path, bbox=None, coerce_f64_to_f32=True + dsm_path, bbox=None, coerce_f64_to_f32=True ) - logger.info("DSM loaded from %s", model_configs.dsm_path) + logger.info("DSM loaded from %s", dsm_path) self.scale = 1 / self.trf_arr[1] self.rows = self.dsm.shape[0] self.cols = self.dsm.shape[1] + one32 = np.float32(1.0) + zero32 = np.float32(0.0) + self.dsm = np.ascontiguousarray(self.dsm, dtype=np.float32) # TODO: is this needed? # if self.dsm.min() < 0: # dsmraise = np.abs(self.dsm.min()) @@ -564,9 +600,7 @@ def __init__( # WALLS # heights - self.wallheight, wh_trf, wh_crs, _ = common.load_raster( - model_configs.wh_path, bbox=None, coerce_f64_to_f32=True - ) + self.wallheight, wh_trf, wh_crs, _ = common.load_raster(wh_path, bbox=None, coerce_f64_to_f32=True) if not self.wallheight.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for wall heights and DSM.") if not np.allclose(self.trf_arr, wh_trf): @@ -574,10 +608,9 @@ def __init__( if not self.crs_wkt == wh_crs: raise ValueError("Mismatching CRS for wall heights and DSM.") logger.info("Wall heights loaded") + self.wallheight = np.ascontiguousarray(self.wallheight, dtype=np.float32) # aspects - self.wallaspect, wa_trf, wa_crs, _ = common.load_raster( - model_configs.wa_path, bbox=None, coerce_f64_to_f32=True - ) + self.wallaspect, wa_trf, wa_crs, _ = common.load_raster(wa_path, bbox=None, coerce_f64_to_f32=True) if not self.wallaspect.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for wall aspects and DSM.") if not np.allclose(self.trf_arr, wa_trf): @@ -585,6 +618,7 @@ def __init__( if not self.crs_wkt == wa_crs: raise ValueError("Mismatching CRS for wall aspects and DSM.") logger.info("Wall aspects loaded") + self.wallaspect = np.ascontiguousarray(self.wallaspect, dtype=np.float32) # DEM # TODO: Is DEM always provided? @@ -598,6 +632,7 @@ def __init__( if not np.allclose(self.trf_arr, dem_trf): raise ValueError("Mismatching spatial transform for DEM and CDSM.") logger.info("DEM loaded from %s", model_configs.dem_path) + self.dem = np.ascontiguousarray(self.dem, dtype=np.float32) # dem[dem == dem_nd_val] = 0.0 # TODO: Check if this is needed re DSM ramifications # if dem.min() < 0: @@ -618,6 +653,7 @@ def __init__( if not np.allclose(self.trf_arr, vegdsm_trf): raise ValueError("Mismatching spatial transform for DSM and CDSM.") logger.info("Vegetation DSM loaded from %s", model_configs.cdsm_path) + self.cdsm = np.ascontiguousarray(self.cdsm, dtype=np.float32) # Tree DSM if model_configs.tdsm_path: self.tdsm, vegdsm2_trf, vegdsm2_crs, _ = common.load_raster( @@ -630,6 +666,7 @@ def __init__( if not np.allclose(self.trf_arr, vegdsm2_trf): raise ValueError("Mismatching spatial transform for DSM and CDSM.") logger.info("Tree DSM loaded from %s", model_configs.tdsm_path) + self.tdsm = np.ascontiguousarray(self.tdsm, dtype=np.float32) else: self.tdsm = None else: @@ -646,20 +683,28 @@ def __init__( amax_local_window_m, amax_local_perc, ) + self.dsm = np.ascontiguousarray(self.dsm, dtype=np.float32) + if self.dem is not None: + self.dem = np.ascontiguousarray(self.dem, dtype=np.float32) + if self.cdsm is not None: + self.cdsm = np.ascontiguousarray(self.cdsm, dtype=np.float32) + if self.tdsm is not None: + self.tdsm = np.ascontiguousarray(self.tdsm, dtype=np.float32) # bushes etc if model_configs.use_veg_dem: + transmissivity = np.float32(model_params.Tree_settings.Value.Transmissivity) self.bush = np.logical_not(self.tdsm * self.cdsm) * self.cdsm - self.svfbuveg = svf_data.svf - (1.0 - svf_data.svf_veg) * ( - 1.0 - model_params.Tree_settings.Value.Transmissivity - ) + self.bush = np.ascontiguousarray(self.bush, dtype=np.float32) + self.svfbuveg = svf_data.svf - (one32 - svf_data.svf_veg) * (one32 - transmissivity) + self.svfbuveg = np.ascontiguousarray(self.svfbuveg, dtype=np.float32) else: logger.info("Vegetation DEM not used; vegetation arrays set to None.") self.bush = np.zeros([self.rows, self.cols], dtype=np.float32) - self.svfbuveg = svf_data.svf + self.svfbuveg = np.ascontiguousarray(svf_data.svf, dtype=np.float32) common.save_raster( - model_configs.output_dir + "/input-dsm.tif", + output_dir + "/input-dsm.tif", self.dsm, self.trf_arr, self.crs_wkt, @@ -667,7 +712,7 @@ def __init__( coerce_f64_to_f32=True, ) common.save_raster( - model_configs.output_dir + "/input-cdsm.tif", + output_dir + "/input-cdsm.tif", self.cdsm, self.trf_arr, self.crs_wkt, @@ -675,7 +720,7 @@ def __init__( coerce_f64_to_f32=True, ) common.save_raster( - model_configs.output_dir + "/input-tdsm.tif", + output_dir + "/input-tdsm.tif", self.tdsm, self.trf_arr, self.crs_wkt, @@ -683,7 +728,7 @@ def __init__( coerce_f64_to_f32=True, ) common.save_raster( - model_configs.output_dir + "/input-svfbuveg.tif", + output_dir + "/input-svfbuveg.tif", self.svfbuveg, self.trf_arr, self.crs_wkt, @@ -691,7 +736,7 @@ def __init__( coerce_f64_to_f32=True, ) common.save_raster( - model_configs.output_dir + "/input-bush.tif", + output_dir + "/input-bush.tif", self.bush, self.trf_arr, self.crs_wkt, @@ -701,7 +746,10 @@ def __init__( # Land cover if model_configs.use_landcover: - lc_path_str = str(common.check_path(model_configs.lc_path)) + lc_path = model_configs.lc_path + if lc_path is None: + raise ValueError("Land cover path must be provided when use_landcover is True.") + lc_path_str = str(common.check_path(lc_path)) self.lcgrid, lc_trf, lc_crs, _ = common.load_raster(lc_path_str, bbox=None, coerce_f64_to_f32=True) if not self.lcgrid.shape == self.dsm.shape: raise ValueError("Mismatching raster shapes for land cover and DSM.") @@ -709,7 +757,8 @@ def __init__( raise ValueError("Mismatching CRS for land cover and DSM.") if not np.allclose(self.trf_arr, lc_trf): raise ValueError("Mismatching spatial transform for land cover and DSM.") - logger.info("Land cover loaded from %s", model_configs.lc_path) + logger.info("Land cover loaded from %s", lc_path) + self.lcgrid = np.ascontiguousarray(self.lcgrid, dtype=np.float32) else: self.lcgrid = None logger.info("Land cover not used; lcgrid set to None.") @@ -725,14 +774,21 @@ def __init__( buildings[buildings == 4] = 1 buildings[buildings == 3] = 1 buildings[buildings == 2] = 0 - self.buildings = buildings + self.buildings = np.ascontiguousarray(buildings, dtype=np.float32) logger.info("Buildings raster created from land cover data.") elif model_configs.use_dem_for_buildings: - buildings = np.where(~np.isnan(self.dem) & ~np.isnan(self.dsm), self.dsm - self.dem, 0.0) + if self.dem is None: + raise ValueError("DEM raster must be available when use_dem_for_buildings is True.") + height_diff = self.dsm - self.dem + buildings = np.where( + ~np.isnan(self.dem) & ~np.isnan(self.dsm), + height_diff, + zero32, + ) # TODO: Check intended logic here - 1 vs 0 buildings[buildings < 2.0] = 1.0 buildings[buildings >= 2.0] = 0.0 - self.buildings = buildings + self.buildings = np.ascontiguousarray(buildings, dtype=np.float32) logger.info("Buildings raster created from DSM and DEM data.") else: self.buildings = None @@ -740,14 +796,14 @@ def __init__( # Save buildings raster if requested if self.buildings is not None and model_configs.save_buildings: common.save_raster( - model_configs.output_dir + "/buildings.tif", + output_dir + "/buildings.tif", self.buildings, self.trf_arr, self.crs_wkt, self.nd_val, coerce_f64_to_f32=True, ) - logger.info("Buildings raster saved to %s/buildings.tif", model_configs.output_dir) + logger.info("Buildings raster saved to %s/buildings.tif", output_dir) class ShadowMatrices: @@ -915,7 +971,7 @@ class WallsData: timeStep: int walls_scheme: np.ndarray dirwalls_scheme: np.ndarray - met_for_xarray: Optional[Tuple[Any]] + met_for_xarray: Optional[Any] def __init__( self,