diff --git a/functions/SOLWEIGpython/Lside_veg.py b/functions/SOLWEIGpython/Lside_veg.py new file mode 100644 index 0000000..da67dc7 --- /dev/null +++ b/functions/SOLWEIGpython/Lside_veg.py @@ -0,0 +1,437 @@ +from __future__ import absolute_import +import numpy as np +from .Lvikt_veg import Lvikt_veg + + +def Lside_veg_v2022a( + svfS, + svfW, + svfN, + svfE, + svfEveg, + svfSveg, + svfWveg, + svfNveg, + svfEaveg, + svfSaveg, + svfWaveg, + svfNaveg, + azimuth, + altitude, + Ta, + Tw, + SBC, + ewall, + Ldown, + esky, + t, + F_sh, + CI, + LupE, + LupS, + LupW, + LupN, + anisotropic_longwave, +): + + # This m-file is the current one that estimates L from the four cardinal points 20100414 + + # Building height angle from svf + svfalfaE = np.arcsin(np.exp((np.log(1 - svfE)) / 2)) + svfalfaS = np.arcsin(np.exp((np.log(1 - svfS)) / 2)) + svfalfaW = np.arcsin(np.exp((np.log(1 - svfW)) / 2)) + svfalfaN = np.arcsin(np.exp((np.log(1 - svfN)) / 2)) + + vikttot = 4.4897 + aziW = azimuth + t + aziN = azimuth - 90 + t + aziE = azimuth - 180 + t + aziS = azimuth - 270 + t + + F_sh = 2 * F_sh - 1 # (cylindric_wedge scaled 0-1) + + c = 1 - CI + Lsky_allsky = esky * SBC * ((Ta + 273.15) ** 4) * (1 - c) + c * SBC * ( + (Ta + 273.15) ** 4 + ) + + ## Least + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfE, svfEveg, svfEaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaE) + betaB = np.arctan(np.tan((svfalfaE) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaE)*(1+F_sh)) #TODO This should be considered in future versions + if (azimuth > (180 - t)) and (azimuth <= (360 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziE * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + # Longwave from ground (see Lcyl_v2022a for remaining fluxes) + if anisotropic_longwave == 1: + Lground = LupE * 0.5 + Least = Lground + else: + Lsky = ((svfE + svfEveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lground = LupE * 0.5 + Lrefl = (Ldown + LupE) * (viktrefl) * (1 - ewall) * 0.5 + Least = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + ## Lsouth + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfS, svfSveg, svfSaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaS) + betaB = np.arctan(np.tan((svfalfaS) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaS)*(1+F_sh)) + if (azimuth <= (90 - t)) or (azimuth > (270 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziS * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + # Longwave from ground (see Lcyl_v2022a for remaining fluxes) + if anisotropic_longwave == 1: + Lground = LupS * 0.5 + Lsouth = Lground + else: + Lsky = ((svfS + svfSveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lground = LupS * 0.5 + Lrefl = (Ldown + LupS) * (viktrefl) * (1 - ewall) * 0.5 + Lsouth = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + ## Lwest + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfW, svfWveg, svfWaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaW) + betaB = np.arctan(np.tan((svfalfaW) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaW)*(1+F_sh)) + if (azimuth > (360 - t)) or (azimuth <= (180 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziW * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + # Longwave from ground (see Lcyl_v2022a for remaining fluxes) + if anisotropic_longwave == 1: + Lground = LupW * 0.5 + Lwest = Lground + else: + Lsky = ((svfW + svfWveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lground = LupW * 0.5 + Lrefl = (Ldown + LupW) * (viktrefl) * (1 - ewall) * 0.5 + Lwest = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + ## Lnorth + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfN, svfNveg, svfNaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaN) + betaB = np.arctan(np.tan((svfalfaN) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaN)*(1+F_sh)) + if (azimuth > (90 - t)) and (azimuth <= (270 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziN * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + # Longwave from ground (see Lcyl_v2022a for remaining fluxes) + if anisotropic_longwave == 1: + Lground = LupN * 0.5 + Lnorth = Lground + else: + Lsky = ((svfN + svfNveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lground = LupN * 0.5 + Lrefl = (Ldown + LupN) * (viktrefl) * (1 - ewall) * 0.5 + Lnorth = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + return Least, Lsouth, Lwest, Lnorth + + +def Lside_veg_v2026( + svfS, + svfW, + svfN, + svfE, + svfEveg, + svfSveg, + svfWveg, + svfNveg, + svfEaveg, + svfSaveg, + svfWaveg, + svfNaveg, + azimuth, + altitude, + Ta, + Tw, + SBC, + ewall, + Ldown, + esky, + t, + F_sh, + CI, + anisotropic_longwave, +): + + # This m-file is the current one that estimates L from the four cardinal points 20100414 + + # Building height angle from svf + svfalfaE = np.arcsin(np.exp((np.log(1 - svfE)) / 2)) + svfalfaS = np.arcsin(np.exp((np.log(1 - svfS)) / 2)) + svfalfaW = np.arcsin(np.exp((np.log(1 - svfW)) / 2)) + svfalfaN = np.arcsin(np.exp((np.log(1 - svfN)) / 2)) + + vikttot = 4.4897 + aziW = azimuth + t + aziN = azimuth - 90 + t + aziE = azimuth - 180 + t + aziS = azimuth - 270 + t + + F_sh = 2 * F_sh - 1 # (cylindric_wedge scaled 0-1) + + c = 1 - CI + Lsky_allsky = esky * SBC * ((Ta + 273.15) ** 4) * (1 - c) + c * SBC * ( + (Ta + 273.15) ** 4 + ) + + ## Least + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfE, svfEveg, svfEaveg, vikttot + ) + + if anisotropic_longwave == 1: + Least = np.zeros_like(Ldown) + Lnorth = np.zeros_like(Ldown) + Lwest = np.zeros_like(Ldown) + Lsouth = np.zeros_like(Ldown) + + return Least, Lsouth, Lwest, Lnorth + else: + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaE) + betaB = np.arctan(np.tan((svfalfaE) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaE)*(1+F_sh)) #TODO This should be considered in future versions + if (azimuth > (180 - t)) and (azimuth <= (360 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziE * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + Lsky = ((svfE + svfEveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lrefl = Ldown * (viktrefl) * (1 - ewall) * 0.5 + Least = Lsky + Lwallsun + Lwallsh + Lveg + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + ## Lsouth + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfS, svfSveg, svfSaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaS) + betaB = np.arctan(np.tan((svfalfaS) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaS)*(1+F_sh)) + if (azimuth <= (90 - t)) or (azimuth > (270 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziS * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + Lsky = ((svfS + svfSveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lrefl = Ldown * (viktrefl) * (1 - ewall) * 0.5 + Lsouth = Lsky + Lwallsun + Lwallsh + Lveg + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + ## Lwest + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfW, svfWveg, svfWaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaW) + betaB = np.arctan(np.tan((svfalfaW) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaW)*(1+F_sh)) + if (azimuth > (360 - t)) or (azimuth <= (180 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziW * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + Lsky = ((svfW + svfWveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lrefl = Ldown * (viktrefl) * (1 - ewall) * 0.5 + Lwest = Lsky + Lwallsun + Lwallsh + Lveg + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + ## Lnorth + [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( + svfN, svfNveg, svfNaveg, vikttot + ) + + if altitude > 0: # daytime + alfaB = np.arctan(svfalfaN) + betaB = np.arctan(np.tan((svfalfaN) * F_sh)) + betasun = ((alfaB - betaB) / 2) + betaB + # betasun = np.arctan(0.5*np.tan(svfalfaN)*(1+F_sh)) + if (azimuth > (90 - t)) and (azimuth <= (270 - t)): + Lwallsun = ( + SBC + * ewall + * ((Ta + 273.15 + Tw * np.sin(aziN * (np.pi / 180))) ** 4) + * viktwall + * (1 - F_sh) + * np.cos(betasun) + * 0.5 + ) + Lwallsh = ( + SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 + ) + else: + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + else: # nighttime + Lwallsun = 0 + Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 + + Lsky = ((svfN + svfNveg - 1) * Lsky_allsky) * viktsky * 0.5 + Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 + Lrefl = Ldown * (viktrefl) * (1 - ewall) * 0.5 + Lnorth = Lsky + Lwallsun + Lwallsh + Lveg + Lrefl + + # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl viktveg viktwall viktsky + + return Least, Lsouth, Lwest, Lnorth diff --git a/functions/SOLWEIGpython/Lside_veg_v2015a.py b/functions/SOLWEIGpython/Lside_veg_v2015a.py deleted file mode 100644 index 6700501..0000000 --- a/functions/SOLWEIGpython/Lside_veg_v2015a.py +++ /dev/null @@ -1,215 +0,0 @@ -from __future__ import absolute_import -import numpy as np -from .Lvikt_veg import Lvikt_veg - - -def Lside_veg_v2015a( - svfS, - svfW, - svfN, - svfE, - svfEveg, - svfSveg, - svfWveg, - svfNveg, - svfEaveg, - svfSaveg, - svfWaveg, - svfNaveg, - azimuth, - altitude, - Ta, - Tw, - SBC, - ewall, - Ldown, - esky, - t, - F_sh, - CI, - LupE, - LupS, - LupW, - LupN, -): - - # This m-file is the current one that estimates L from the four cardinal - # points 20100414 - - # Building height angle from svf - svfalfaE = np.arcsin(np.exp((np.log(1 - svfE)) / 2)) - svfalfaS = np.arcsin(np.exp((np.log(1 - svfS)) / 2)) - svfalfaW = np.arcsin(np.exp((np.log(1 - svfW)) / 2)) - svfalfaN = np.arcsin(np.exp((np.log(1 - svfN)) / 2)) - - vikttot = 4.4897 - aziW = azimuth + t - aziN = azimuth - 90 + t - aziE = azimuth - 180 + t - aziS = azimuth - 270 + t - - F_sh = 2 * F_sh - 1 # (cylindric_wedge scaled 0-1) - - c = 1 - CI - Lsky_allsky = esky * SBC * ((Ta + 273.15) ** 4) * (1 - c) + c * SBC * ( - (Ta + 273.15) ** 4 - ) - - # Least - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfE, svfEveg, svfEaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaE) - betaB = np.arctan(np.tan((svfalfaE) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaE)*(1+F_sh)) #TODO This should - # be considered in future versions - if (azimuth > (180 - t)) and (azimuth <= (360 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziE * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - Lsky = ((svfE + svfEveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupE * 0.5 - Lrefl = (Ldown + LupE) * (viktrefl) * (1 - ewall) * 0.5 - Least = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - # Lsouth - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfS, svfSveg, svfSaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaS) - betaB = np.arctan(np.tan((svfalfaS) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaS)*(1+F_sh)) - if (azimuth <= (90 - t)) or (azimuth > (270 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziS * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - Lsky = ((svfS + svfSveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupS * 0.5 - Lrefl = (Ldown + LupS) * (viktrefl) * (1 - ewall) * 0.5 - Lsouth = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - # Lwest - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfW, svfWveg, svfWaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaW) - betaB = np.arctan(np.tan((svfalfaW) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaW)*(1+F_sh)) - if (azimuth > (360 - t)) or (azimuth <= (180 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziW * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - Lsky = ((svfW + svfWveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupW * 0.5 - Lrefl = (Ldown + LupW) * (viktrefl) * (1 - ewall) * 0.5 - Lwest = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - # Lnorth - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfN, svfNveg, svfNaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaN) - betaB = np.arctan(np.tan((svfalfaN) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaN)*(1+F_sh)) - if (azimuth > (90 - t)) and (azimuth <= (270 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziN * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - Lsky = ((svfN + svfNveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupN * 0.5 - Lrefl = (Ldown + LupN) * (viktrefl) * (1 - ewall) * 0.5 - Lnorth = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - return Least, Lsouth, Lwest, Lnorth diff --git a/functions/SOLWEIGpython/Lside_veg_v2022a.py b/functions/SOLWEIGpython/Lside_veg_v2022a.py deleted file mode 100644 index b124539..0000000 --- a/functions/SOLWEIGpython/Lside_veg_v2022a.py +++ /dev/null @@ -1,236 +0,0 @@ -from __future__ import absolute_import -import numpy as np -from .Lvikt_veg import Lvikt_veg - - -def Lside_veg_v2022a( - svfS, - svfW, - svfN, - svfE, - svfEveg, - svfSveg, - svfWveg, - svfNveg, - svfEaveg, - svfSaveg, - svfWaveg, - svfNaveg, - azimuth, - altitude, - Ta, - Tw, - SBC, - ewall, - Ldown, - esky, - t, - F_sh, - CI, - LupE, - LupS, - LupW, - LupN, - anisotropic_longwave, -): - - # This m-file is the current one that estimates L from the four cardinal - # points 20100414 - - # Building height angle from svf - svfalfaE = np.arcsin(np.exp((np.log(1 - svfE)) / 2)) - svfalfaS = np.arcsin(np.exp((np.log(1 - svfS)) / 2)) - svfalfaW = np.arcsin(np.exp((np.log(1 - svfW)) / 2)) - svfalfaN = np.arcsin(np.exp((np.log(1 - svfN)) / 2)) - - vikttot = 4.4897 - aziW = azimuth + t - aziN = azimuth - 90 + t - aziE = azimuth - 180 + t - aziS = azimuth - 270 + t - - F_sh = 2 * F_sh - 1 # (cylindric_wedge scaled 0-1) - - c = 1 - CI - Lsky_allsky = esky * SBC * ((Ta + 273.15) ** 4) * (1 - c) + c * SBC * ( - (Ta + 273.15) ** 4 - ) - - # Least - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfE, svfEveg, svfEaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaE) - betaB = np.arctan(np.tan((svfalfaE) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaE)*(1+F_sh)) #TODO This should - # be considered in future versions - if (azimuth > (180 - t)) and (azimuth <= (360 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziE * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - # Longwave from ground (see Lcyl_v2022a for remaining fluxes) - if anisotropic_longwave == 1: - Lground = LupE * 0.5 - Least = Lground - else: - Lsky = ((svfE + svfEveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupE * 0.5 - Lrefl = (Ldown + LupE) * (viktrefl) * (1 - ewall) * 0.5 - Least = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - # Lsouth - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfS, svfSveg, svfSaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaS) - betaB = np.arctan(np.tan((svfalfaS) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaS)*(1+F_sh)) - if (azimuth <= (90 - t)) or (azimuth > (270 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziS * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - # Longwave from ground (see Lcyl_v2022a for remaining fluxes) - if anisotropic_longwave == 1: - Lground = LupS * 0.5 - Lsouth = Lground - else: - Lsky = ((svfS + svfSveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupS * 0.5 - Lrefl = (Ldown + LupS) * (viktrefl) * (1 - ewall) * 0.5 - Lsouth = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - # Lwest - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfW, svfWveg, svfWaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaW) - betaB = np.arctan(np.tan((svfalfaW) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaW)*(1+F_sh)) - if (azimuth > (360 - t)) or (azimuth <= (180 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziW * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - # Longwave from ground (see Lcyl_v2022a for remaining fluxes) - if anisotropic_longwave == 1: - Lground = LupW * 0.5 - Lwest = Lground - else: - Lsky = ((svfW + svfWveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupW * 0.5 - Lrefl = (Ldown + LupW) * (viktrefl) * (1 - ewall) * 0.5 - Lwest = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - # Lnorth - [viktveg, viktwall, viktsky, viktrefl] = Lvikt_veg( - svfN, svfNveg, svfNaveg, vikttot - ) - - if altitude > 0: # daytime - alfaB = np.arctan(svfalfaN) - betaB = np.arctan(np.tan((svfalfaN) * F_sh)) - betasun = ((alfaB - betaB) / 2) + betaB - # betasun = np.arctan(0.5*np.tan(svfalfaN)*(1+F_sh)) - if (azimuth > (90 - t)) and (azimuth <= (270 - t)): - Lwallsun = ( - SBC - * ewall - * ((Ta + 273.15 + Tw * np.sin(aziN * (np.pi / 180))) ** 4) - * viktwall - * (1 - F_sh) - * np.cos(betasun) - * 0.5 - ) - Lwallsh = ( - SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * F_sh * 0.5 - ) - else: - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - else: # nighttime - Lwallsun = 0 - Lwallsh = SBC * ewall * ((Ta + 273.15) ** 4) * viktwall * 0.5 - - # Longwave from ground (see Lcyl_v2022a for remaining fluxes) - if anisotropic_longwave == 1: - Lground = LupN * 0.5 - Lnorth = Lground - else: - Lsky = ((svfN + svfNveg - 1) * Lsky_allsky) * viktsky * 0.5 - Lveg = SBC * ewall * ((Ta + 273.15) ** 4) * viktveg * 0.5 - Lground = LupN * 0.5 - Lrefl = (Ldown + LupN) * (viktrefl) * (1 - ewall) * 0.5 - Lnorth = Lsky + Lwallsun + Lwallsh + Lveg + Lground + Lrefl - - # clear alfaB betaB betasun Lsky Lwallsh Lwallsun Lveg Lground Lrefl - # viktveg viktwall viktsky - - return Least, Lsouth, Lwest, Lnorth diff --git a/functions/SOLWEIGpython/Solweig_2026a_calc_forprocessing.py b/functions/SOLWEIGpython/Solweig_2026a_calc_forprocessing.py new file mode 100644 index 0000000..f30d7e3 --- /dev/null +++ b/functions/SOLWEIGpython/Solweig_2026a_calc_forprocessing.py @@ -0,0 +1,962 @@ +""" +@author Fredrik Lindberg, University of Gothenburg +""" + +from __future__ import absolute_import + +import numpy as np +import matplotlib.pyplot as plt +from .daylen import daylen +from ...util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import ( + clearnessindex_2013b, +) +from ...util.SEBESOLWEIGCommonFiles.diffusefraction import diffusefraction +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 .cylindric_wedge import cylindric_wedge +from .TsWaveDelay_2015a import TsWaveDelay_2015a +from .Kup_veg_2015a import Kup_veg_2015a + +# 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 + +# Anisotropic longwave +from .Lcyl_v2022a import Lcyl_v2022a +from .Lside_veg import Lside_veg_v2022a, Lside_veg_v2026 +from .anisotropic_sky import anisotropic_sky as ani_sky +from .patch_radiation import patch_steradians +from copy import deepcopy +import time + +# Wall surface temperature scheme +from .wall_surface_temperature import wall_surface_temperature + +# Ground surface temperature +from .ground_surface import ( + surfaceTemperature_calc, + outgoingLongwave_calc, +) + + +def Solweig_2026a_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, + diffsh, + shmat, + vegshmat, + vbshvegshmat, + anisotropic_sky, + asvf, + patch_option, + voxelMaps, + voxelTable, + ws, + wallScheme, + timeStep, + steradians, + walls_scheme, + dirwalls_scheme, + groundScheme, + outgoingLW, + Tg, + Rn, + Rn_past, + G, + Tm, + cap_grid, + diff_grid, + a1_grid, + a2_grid, + a3_grid, + shadow_past, +): + """ + This is the core function of the SOLWEIG model + 2016-Aug-28 + Fredrik Lindberg, fredrikl@gvc.gu.se + Goteborg Urban Climate Group + Gothenburg University + + :Input variables: + dsm = digital surface model + scale = height to pixel size (2m pixel gives scale = 0.5) + svf,svfN,svfW,svfE,svfS = SVFs for building and ground + svfveg,svfNveg,svfEveg,svfSveg,svfWveg = Veg SVFs blocking sky + svfaveg,svfEaveg,svfSaveg,svfWaveg,svfNaveg = Veg SVFs blocking buildings + vegdem = Vegetation canopy DSM + vegdem2 = Vegetation trunk zone DSM + albedo_b = building wall albedo + absK = human absorption coefficient for shortwave radiation + absL = human absorption coefficient for longwave radiation + ewall = Emissivity of building walls + Fside = The angular factors between a person and the surrounding surfaces + Fup = The angular factors between a person and the surrounding surfaces + Fcyl = The angular factors between a culidric person and the surrounding surfaces + altitude = Sun altitude (degree) + azimuth = Sun azimuth (degree) + zen = Sun zenith angle (radians) + jday = day of year + usevegdem = use vegetation scheme + onlyglobal = calculate dir and diff from global shortwave (Reindl et al. 1990) + buildings = Boolena grid to identify building pixels + location = geographic location + height = height of measurements point (center of gravity of human) + psi = 1 - Transmissivity of shortwave through vegetation + landcover = use landcover scheme !!!NEW IN 2015a!!! + lc_grid = grid with landcoverclasses + lc_class = table with landcover properties + dectime = decimal time + altmax = maximum sun altitude + dirwalls = aspect of walls + walls = one pixel row outside building footprint. height of building walls + cyl = consider man as cylinder instead of cude + elvis = dummy + Ta = air temp + RH + radG = global radiation + radD = diffuse + radI = direct + P = pressure + amaxvalue = max height of buildings + bush = grid representing bushes + Twater = temperature of water (daily) + 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, + 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.0 + + # Stefan Bolzmans Constant + SBC = 5.67051e-8 + + # Degrees to radians + deg2rad = np.pi / 180 + + # Find sunrise decimal hour - new from 2014a + _, _, _, SNUP = daylen(jday, location["latitude"]) + + # Vapor pressure in hPa + 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 # # # # # # + # 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.0, radG, location, P + ) + if (CI > 1) or (CI == np.inf): + 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.0, radG, location, P + ) + if (CI > 1) or (CI == np.inf): + CI = 1 + + radI, radD = diffusefraction(radG, altitude, Kt, Ta, RH) + + # Diffuse Radiation + # Anisotropic Diffuse Radiation after Perez et al. 1993 + if anisotropic_sky == 1: + patchchoice = 1 + zenDeg = zen * (180 / np.pi) + # Relative luminance + 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)) + for idx in range(lv.shape[0]): + aniLum += diffsh[:, :, idx] * lv[idx, 2] + + dRad = ( + aniLum * radD + ) # Total diffuse radiation from sky into each cell + else: + dRad = radD * svfbuveg + patchchoice = 1 + lv = None + + # 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.0, + walls_scheme, + dirwalls_scheme * np.pi / 180.0, + ) + ) + shadow = sh - (1 - vegsh) * (1 - psi) + else: + sh, wallsh, wallsun, facesh, facesun, wallsh_ = ( + shadowingfunction_wallheight_13( + dsm, + azimuth, + altitude, + scale, + walls, + dirwalls * np.pi / 180.0, + walls_scheme, + dirwalls_scheme * np.pi / 180.0, + ) + ) + shadow = sh + + # 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 + + # New estimation of Tgwall with reduction for non-clear situation based on Reindl et al. 1990 + radI0, _ = diffusefraction(I0, altitude, 1.0, Ta, RH) + radG0 = radI0 * (np.sin(altitude * deg2rad)) + _ + corr = ( + 0.1473 * np.log(90 - (zen / np.pi * 180)) + 0.3454 + ) # 20070329 correction of lat, Lindberg et al. 2008 + CI_TgG = (radG / radG0) + (1 - corr) + if (CI_TgG > 1) or (CI_TgG == np.inf): + CI_TgG = 1 + + Tgampwall = TgK_wall * altmax + Tstart_wall + 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 + Tgwall = 0 + Tgwall = Tgwall * CI_TgG + + # # # # Kdown # # # # + 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)) + + # # # # 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 = 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? + + # Surface temperature parameterization during daytime + if groundScheme == 1: + # calculate the ground surface temperature, and relevant heat fluxes + Tg, Rn, Rn_past, G = surfaceTemperature_calc( + Kdown, + Ldown, + Rn, + Rn_past, + G, + Tg, + Tm, + alb_grid, + emis_grid, + cap_grid, + diff_grid, + lc_grid, + a1_grid, + a2_grid, + a3_grid, + timeStep, + RH, + shadow, + shadow_past, + ) + + else: + # using max sun alt instead of dfm + Tgamp = TgK * altmax + Tstart # Fixed 2021 + Tgdiff = Tgamp * np.sin( + ( + ((dectime - np.floor(dectime)) - SNUP / 24) + / (TmaxLST / 24 - SNUP / 24) + ) + * np.pi + / 2 + ) # 2015 a, based on max sun altitude + + Tgdiff = Tgdiff * CI_TgG # new estimation + + # For Tg output in POIs + TgTemp = Tgdiff * shadow + Ta + _, timeadd, Tg = TsWaveDelay_2015a( + TgTemp, firstdaytime, timeadd, timestepdec, Tg + ) # timeadd only here v2021a + + if landcover == 1: + Tg[Tg < 0] = ( + 0 # temporary for removing low Tg during morning 20130205 + ) + + # Calculate the outgoing longwave radiation + if outgoingLW == 1: + # According to the solid angle parameterization + # # # # Lup, daytime # # # # + ( + Lup, + gvfalb, + gvfalbnosh, + LupE, + gvfalbE, + gvfalbnoshE, + LupS, + gvfalbS, + gvfalbnoshS, + LupW, + gvfalbW, + gvfalbnoshW, + LupN, + gvfalbN, + gvfalbnoshN, + gvfLsideW, + gvfLsideS, + gvfLsideE, + gvfLsideN, + ) = outgoingLongwave_calc( + Tg, + Tgwall, + Ta, + Ldown, + emis_grid, + alb_grid, + buildings, + shadow, + wallsun, + walls, + rows, + cols, + 1 / scale, + ) + + else: + ### 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, + ) + + # # # # Lup, daytime # # # # + # Surface temperature wave delay - new as from 2014a + Lup, timeaddnotused, Tgmap1 = TsWaveDelay_2015a( + gvfLup, firstdaytime, timeadd, timestepdec, Tgmap1 + ) + LupE, timeaddnotused, Tgmap1E = TsWaveDelay_2015a( + gvfLupE, firstdaytime, timeadd, timestepdec, Tgmap1E + ) + 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 + ) + + # # # # Kup # # # # + 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 # # # # # # # # + # 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)) + shadow = np.zeros((rows, cols)) + CI_TgG = deepcopy(CI) + dRad = np.zeros((rows, cols)) + Kside = np.zeros((rows, cols)) + wallsun = np.zeros((rows, cols)) + + Tgwall = 0 + + # # # # 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 = 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? + + # Surface temperature parameterization + if groundScheme == 1: + # calculate the ground surface temperature, and relevant heat fluxes + Tg, Rn, Rn_past, G = surfaceTemperature_calc( + Kdown, + Ldown, + Rn, + Rn_past, + G, + Tg, + Tm, + alb_grid, + emis_grid, + cap_grid, + diff_grid, + lc_grid, + a1_grid, + a2_grid, + a3_grid, + timeStep, + RH, + shadow, + shadow_past, + ) + + else: + # In the old scheme the ground surface temperature is equal to the air temperature during nighttime + Tg = np.ones((rows, cols)) * Ta + + # Calculate the outgoing longwave radiation + if outgoingLW == 1: + # According to the solid angle parameterization + # # # # Lup, daytime # # # # + ( + Lup, + gvfalb, + gvfalbnosh, + LupE, + gvfalbE, + gvfalbnoshE, + LupS, + gvfalbS, + gvfalbnoshS, + LupW, + gvfalbW, + gvfalbnoshW, + LupN, + gvfalbN, + gvfalbnoshN, + gvfLsideW, + gvfLsideS, + gvfLsideE, + gvfLsideN, + ) = outgoingLongwave_calc( + Tg, + Tgwall, + Ta, + Ldown, + emis_grid, + alb_grid, + buildings, + shadow, + wallsun, + walls, + rows, + cols, + 1 / scale, + ) + + else: + # # # # Lup, nighttime # # # # + Lup = SBC * emis_grid * ((Knight + Tg + 273.15) ** 4) + LupE = Lup + LupS = Lup + LupW = Lup + LupN = Lup + + I0 = 0 + timeadd = 0 + firstdaytime = 1 + + # # # # Lside # # # # + if groundScheme == 1: + Least = np.copy(gvfLsideE) + Lsouth = np.copy(gvfLsideS) + Lwest = np.copy(gvfLsideW) + Lnorth = np.copy(gvfLsideN) + Least_, Lsouth_, Lwest_, Lnorth_ = Lside_veg_v2026( + svfS, + svfW, + svfN, + svfE, + svfEveg, + svfSveg, + svfWveg, + svfNveg, + svfEaveg, + svfSaveg, + svfWaveg, + svfNaveg, + azimuth, + altitude, + Ta, + Tgwall, + SBC, + ewall, + Ldown, + esky, + t, + F_sh, + CI, + anisotropic_sky, + ) + else: + Least = np.zeros_like(Ldown) + Lnorth = np.zeros_like(Ldown) + Lwest = np.zeros_like(Ldown) + Lsouth = np.zeros_like(Ldown) + 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 += Least_ + Lsouth += Lsouth_ + Lwest += Lwest_ + Lnorth += Lnorth_ + Lside = (Lsouth + Lnorth + Least + Lwest) / 4 + + ### Anisotropic sky + 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]) + + x = np.transpose(np.atleast_2d(skyvaultalt)) + y = np.transpose(np.atleast_2d(skyvaultazi)) + z = np.transpose(np.atleast_2d(patch_emissivities)) + + L_patches = np.append(np.append(x, y, axis=1), z, axis=1) + + else: + L_patches = deepcopy(lv) + + # Calculate steradians for patches if it is the first model iteration + if i == 0: + steradians, skyalt, patch_altitude = patch_steradians(L_patches) + + # 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 + + # 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.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, + ) + Lside += Lside_ + else: + Lside_ = np.zeros((rows, cols)) + L_patches = None + + # Box and anisotropic longwave + if cyl == 0 and anisotropic_sky == 1: + Least += Least_ + Lwest += Lwest_ + Lnorth += Lnorth_ + Lsouth += Lsouth_ + + # Calculation of radiant flux density + # 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) + # 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 + ) + # Human body considered as a standing cube + else: + Sstr = absK * ( + (Kdown + Kup) * Fup + (Knorth + Keast + Ksouth + Kwest) * Fside + ) + absL * ( + (Ldown + Lup) * Fup + (Lnorth + Least + Lsouth + Lwest) * Fside + ) + + # # # # Tmrt # # # # + 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): + Least += Least_ + Lwest += Lwest_ + 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, + radI, + radD, + Lside, + L_patches, + CI_TgG, + KsideD, + dRad, + Kside, + steradians, + voxelTable, + Rn, + Rn_past, + Tm, + G, + ) diff --git a/functions/SOLWEIGpython/Solweig_run.py b/functions/SOLWEIGpython/Solweig_run.py index 01b4996..454ba03 100644 --- a/functions/SOLWEIGpython/Solweig_run.py +++ b/functions/SOLWEIGpython/Solweig_run.py @@ -1,1408 +1,1473 @@ -# This is the main function of the SOLWEIG model -# 2025-Mar-21 -# Fredrik Lindberg, fredrikl@gvc.gu.se -# Goteborg Urban Climate Group -# Gothenburg University - -# imports -from __future__ import absolute_import -from ...util.umep_solweig_export_component import read_solweig_config -from ...util.SEBESOLWEIGCommonFiles.Solweig_v2015_metdata_noload import ( - Solweig_2015a_metdata_noload, -) -from ...util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import ( - clearnessindex_2013b, -) - -from ...functions.SOLWEIGpython import Solweig_2025a_calc_forprocessing as so - -# from ...functions.SOLWEIGpython import WriteMetadataSOLWEIG # Not needed -# anymore? -from ...functions.SOLWEIGpython import PET_calculations as p -from ...functions.SOLWEIGpython import UTCI_calculations as utci -from ...functions.SOLWEIGpython.CirclePlotBar import PolarBarPlot -from ...functions.SOLWEIGpython.wall_surface_temperature import load_walls -from ...functions.SOLWEIGpython.wallOfInterest import pointOfInterest -from ...functions.SOLWEIGpython.patch_characteristics import hemispheric_image -from ...functions.SOLWEIGpython.wallsAsNetCDF import walls_as_netcdf -from ...functions.SOLWEIGpython.Tgmaps_v1 import Tgmaps_v1 -from ...functions import wallalgorithms as wa - -import numpy as np -import json -import zipfile -import pandas as pd -import matplotlib.pylab as plt -from shutil import copyfile - -# imports from osgeo/qgis dependency -try: - from osgeo import gdal - from ...util.misc import saveraster, xy2latlon_fromraster - from qgis.core import QgsRasterLayer - import configparser -except BaseException: - pass - -# imports for standalone -try: - from umep import common - from rasterio.transform import xy, rowcol - import pyproj - from tqdm import tqdm - import geopandas as gpd -except BaseException: - pass - - -def config_to_dict(config: configparser.ConfigParser): - def auto_cast(value: str): - """Try to interpret strings as bool, int, or float.""" - v = value.strip().lower() - if v in ("true", "yes", "on"): - return True - if v in ("false", "no", "off"): - return False - if v.isdigit(): - return int(v) - try: - return float(v) - except ValueError: - return value # fallback: leave as string - - return {k: auto_cast(v) for k, v in config.items()} - - -def solweig_run(configPath, feedback): - """ - Input: - configPath : config file including geodata paths and settings. - feedback : To communicate with qgis gui. Set to None if standalone - """ - - # Load config file - configDict = read_solweig_config(configPath) - - # Load parameters settings for SOLWEIG - with open(configDict["para_json_path"], "r") as jsn: - param = json.load(jsn) - - configDict = config_to_dict(configDict) - - standAlone = int(configDict["standalone"]) - - # reading variables from config and parameters that is not yet presented - cyl = int(configDict["cyl"]) - albedo_b = param["Albedo"]["Effective"]["Value"]["Walls"] - ewall = param["Emissivity"]["Value"]["Walls"] - onlyglobal = int(configDict["onlyglobal"]) - elvis = 0.0 - absK = param["Tmrt_params"]["Value"]["absK"] - absL = param["Tmrt_params"]["Value"]["absL"] - - # Load DSM - if standAlone == 1: - dsm, dsm_transf, dsm_crs = common.load_raster( - configDict["filepath_dsm"], bbox=None - ) - scale = 1 / dsm_transf.a - # dsm_height, dsm_width = dsm.shape # y rows by x cols - # y is flipped - so return max for lower row - minx, miny = xy(dsm_transf, dsm.shape[0], 0) - # Define the source and target CRS - source_crs = pyproj.CRS(dsm_crs) - target_crs = pyproj.CRS(4326) # WGS 84 - # Create a transformer object - transformer = pyproj.Transformer.from_crs( - source_crs, target_crs, always_xy=True - ) - # Perform the transformation - lon, lat = transformer.transform(minx, miny) - nd = -9999 # TODO: extract nodatavalue from rasterio - else: - dsm_wkt = QgsRasterLayer(configDict["filepath_dsm"]).crs().toWkt() - gdal_dsm = gdal.Open(configDict["filepath_dsm"]) - lat, lon, scale, minx, miny = xy2latlon_fromraster(dsm_wkt, gdal_dsm) - dsm = gdal_dsm.ReadAsArray().astype(float) - nd = gdal_dsm.GetRasterBand(1).GetNoDataValue() - - rows = dsm.shape[0] - cols = dsm.shape[1] - - # response to issue #85 - dsm[dsm == nd] = 0.0 - if dsm.min() < 0: - dsmraise = np.abs(dsm.min()) - dsm = dsm + dsmraise - else: - dsmraise = 0 - - alt = np.median(dsm) - if alt < 0: - alt = 3 - - # Vegetation - transVeg = param["Tree_settings"]["Value"]["Transmissivity"] - trunkratio = param["Tree_settings"]["Value"]["Trunk_ratio"] - usevegdem = int(configDict["usevegdem"]) - if usevegdem == 1: - if standAlone == 0: - vegdsm = ( - gdal.Open(configDict["filepath_cdsm"]) - .ReadAsArray() - .astype(float) - ) - else: - vegdsm, _, _ = common.load_raster( - configDict["filepath_cdsm"], bbox=None - ) - if configDict["filepath_tdsm"] is not "": - if standAlone == 0: - vegdsm2 = ( - gdal.Open(configDict["filepath_tdsm"]) - .ReadAsArray() - .astype(float) - ) - else: - vegdsm2, _, _ = common.load_raster( - configDict["filepath_tdsm"], bbox=None - ) - else: - vegdsm2 = vegdsm * trunkratio - else: - vegdsm = 0 - vegdsm2 = 0 - - # Land cover - landcover = int(configDict["landcover"]) - if landcover == 1: - if standAlone == 0: - lcgrid = ( - gdal.Open(configDict["filepath_lc"]) - .ReadAsArray() - .astype(float) - ) - else: - lcgrid, _, _ = common.load_raster( - configDict["filepath_lc"], bbox=None - ) - else: - lcgrid = 0 - - # DEM for buildings #TODO: fix nodata in standalone - demforbuild = int(configDict["demforbuild"]) - if demforbuild == 1: - if standAlone == 0: - # .ReadAsArray().astype(float) - gdal_dem = gdal.Open(configDict["filepath_dem"]) - dem = gdal_dem.ReadAsArray().astype(float) - nd = gdal_dem.GetRasterBand(1).GetNoDataValue() - else: - dem, _, _ = common.load_raster( - configDict["filepath_dem"], bbox=None - ) - nd = -9999 # TODO: standAlone nd exposure - - # response to issue and #230 - dem[dem == nd] = 0.0 - if dem.min() < 0: - demraise = np.abs(dem.min()) - dem = dem + demraise - else: - demraise = 0 - - # SVF - zip = zipfile.ZipFile(configDict["input_svf"], "r") - zip.extractall(configDict["working_dir"]) - zip.close() - - if standAlone == 0: - svf = ( - gdal.Open(configDict["working_dir"] + "/svf.tif") - .ReadAsArray() - .astype(float) - ) - svfN = ( - gdal.Open(configDict["working_dir"] + "/svfN.tif") - .ReadAsArray() - .astype(float) - ) - svfS = ( - gdal.Open(configDict["working_dir"] + "/svfS.tif") - .ReadAsArray() - .astype(float) - ) - svfE = ( - gdal.Open(configDict["working_dir"] + "/svfE.tif") - .ReadAsArray() - .astype(float) - ) - svfW = ( - gdal.Open(configDict["working_dir"] + "/svfW.tif") - .ReadAsArray() - .astype(float) - ) - else: - 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 - ) - - if usevegdem == 1: - if standAlone == 0: - svfveg = ( - gdal.Open(configDict["working_dir"] + "/svfveg.tif") - .ReadAsArray() - .astype(float) - ) - svfNveg = ( - gdal.Open(configDict["working_dir"] + "/svfNveg.tif") - .ReadAsArray() - .astype(float) - ) - svfSveg = ( - gdal.Open(configDict["working_dir"] + "/svfSveg.tif") - .ReadAsArray() - .astype(float) - ) - svfEveg = ( - gdal.Open(configDict["working_dir"] + "/svfEveg.tif") - .ReadAsArray() - .astype(float) - ) - svfWveg = ( - gdal.Open(configDict["working_dir"] + "/svfWveg.tif") - .ReadAsArray() - .astype(float) - ) - - svfaveg = ( - gdal.Open(configDict["working_dir"] + "/svfaveg.tif") - .ReadAsArray() - .astype(float) - ) - svfNaveg = ( - gdal.Open(configDict["working_dir"] + "/svfNaveg.tif") - .ReadAsArray() - .astype(float) - ) - svfSaveg = ( - gdal.Open(configDict["working_dir"] + "/svfSaveg.tif") - .ReadAsArray() - .astype(float) - ) - svfEaveg = ( - gdal.Open(configDict["working_dir"] + "/svfEaveg.tif") - .ReadAsArray() - .astype(float) - ) - svfWaveg = ( - gdal.Open(configDict["working_dir"] + "/svfWaveg.tif") - .ReadAsArray() - .astype(float) - ) - else: - 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 - ) - else: - svfveg = np.ones((rows, cols)) - svfNveg = np.ones((rows, cols)) - svfSveg = np.ones((rows, cols)) - svfEveg = np.ones((rows, cols)) - svfWveg = np.ones((rows, cols)) - svfaveg = np.ones((rows, cols)) - svfNaveg = np.ones((rows, cols)) - svfSaveg = np.ones((rows, cols)) - svfEaveg = np.ones((rows, cols)) - svfWaveg = np.ones((rows, cols)) - - tmp = svf + svfveg - 1.0 - tmp[tmp < 0.0] = 0.0 - # %matlab crazyness around 0 - svfalfa = np.arcsin(np.exp((np.log((1.0 - tmp)) / 2.0))) - - if standAlone == 0: - wallheight = ( - gdal.Open(configDict["filepath_wh"]).ReadAsArray().astype(float) - ) - wallaspect = ( - gdal.Open(configDict["filepath_wa"]).ReadAsArray().astype(float) - ) - else: - wallheight, _, _ = common.load_raster( - configDict["filepath_wh"], bbox=None - ) - wallaspect, _, _ = common.load_raster( - configDict["filepath_wa"], bbox=None - ) - - # Metdata - headernum = 1 - delim = " " - Twater = [] - - metdata = np.loadtxt( - configDict["input_met"], skiprows=headernum, delimiter=delim - ) - - location = {"longitude": lon, "latitude": lat, "altitude": alt} - YYYY, altitude, azimuth, zen, jday, leafon, dectime, altmax = ( - Solweig_2015a_metdata_noload(metdata, location, int(configDict["utc"])) - ) - - DOY = metdata[:, 1] - hours = metdata[:, 2] - minu = metdata[:, 3] - Ta = metdata[:, 11] - RH = metdata[:, 10] - radG = metdata[:, 14] - radD = metdata[:, 21] - radI = metdata[:, 22] - P = metdata[:, 12] - Ws = metdata[:, 9] - - # POIs check - if configDict["poi_file"] is not "": # usePOI: - header = ( - "yyyy id it imin dectime altitude azimuth kdir kdiff kglobal kdown kup keast ksouth " - "kwest knorth ldown lup least lsouth lwest lnorth Ta Tg RH Esky Tmrt " - "I0 CI Shadow SVF_b SVF_bv KsideI PET UTCI CI_Tg CI_TgG KsideD Lside diffDown Kside" - ) - # poiname = [] - # self.parameterAsString(parameters, self.POI_FIELD, context) - poi_field = configDict["poi_field"] - if standAlone == 0: - # self.parameterAsStrings(parameters, self.WOI_FIELD, context) - poi_field = configDict["woi_field"] - poisxy, poiname = pointOfInterest( - configDict["poi_file"], poi_field, scale, gdal_dsm - ) - - else: - pois_gdf = gpd.read_file(configDict["poi_file"]) - numfeat = pois_gdf.shape[0] - poisxy = np.zeros((numfeat, 3)) - 999 - for idx, row in pois_gdf.iterrows(): - # TODO: This produce different result since no standalone round - # coordinates - y, x = rowcol( - dsm_transf, - row["geometry"].centroid.x, - row["geometry"].centroid.y, - ) - poiname.append(row[configDict["poi_field"]]) - poisxy[idx, 0] = idx - poisxy[idx, 1] = x - poisxy[idx, 2] = y - - for k in range(0, poisxy.shape[0]): - poi_save = [] # np.zeros((1, 33)) - data_out = ( - configDict["output_dir"] + "/POI_" + str(poiname[k]) + ".txt" - ) - np.savetxt( - data_out, poi_save, delimiter=" ", header=header, comments="" - ) - - # Num format for POI output - numformat = "%d %d %d %d %.5f " + "%.2f " * 36 - - # Other PET variables - sensorheight = param["Wind_Height"]["Value"]["magl"] - age = param["PET_settings"]["Value"]["Age"] - mbody = param["PET_settings"]["Value"]["Weight"] - ht = param["PET_settings"]["Value"]["Height"] - clo = param["PET_settings"]["Value"]["clo"] - activity = param["PET_settings"]["Value"]["Activity"] - sex = param["PET_settings"]["Value"]["Sex"] - else: - poisxy = None - - # Posture settings - if param["Tmrt_params"]["Value"]["posture"] == "Standing": - Fside = param["Posture"]["Standing"]["Value"]["Fside"] - Fup = param["Posture"]["Standing"]["Value"]["Fup"] - height = param["Posture"]["Standing"]["Value"]["height"] - Fcyl = param["Posture"]["Standing"]["Value"]["Fcyl"] - pos = 1 - else: - Fside = param["Posture"]["Sitting"]["Value"]["Fside"] - Fup = param["Posture"]["Sitting"]["Value"]["Fup"] - height = param["Posture"]["Sitting"]["Value"]["height"] - Fcyl = param["Posture"]["Sitting"]["Value"]["Fcyl"] - pos = 0 - - # Radiative surface influence, Rule of thumb by Schmid et al. (1990). - first = np.round(height) - if first == 0.0: - first = 1.0 - second = np.round((height * 20.0)) - - if usevegdem == 1: - # Conifer or deciduous - if configDict["conifer_bool"]: - leafon = np.ones((1, DOY.shape[0])) - else: - leafon = np.zeros((1, DOY.shape[0])) - if ( - param["Tree_settings"]["Value"]["First_day_leaf"] - > param["Tree_settings"]["Value"]["Last_day_leaf"] - ): - leaf_bool = ( - DOY > param["Tree_settings"]["Value"]["First_day_leaf"] - ) | (DOY < param["Tree_settings"]["Value"]["Last_day_leaf"]) - else: - leaf_bool = ( - DOY > param["Tree_settings"]["Value"]["First_day_leaf"] - ) & (DOY < param["Tree_settings"]["Value"]["Last_day_leaf"]) - leafon[0, leaf_bool] = 1 - - # Vegetation transmittivity of shortwave radiation - psi = leafon * transVeg - psi[leafon == 0] = 0.5 - # amaxvalue - vegmax = vegdsm.max() - amaxvalue = dsm.max() - dsm.min() - amaxvalue = np.maximum(amaxvalue, vegmax) - - # Elevation vegdsms if buildingDEM includes ground heights - vegdsm = vegdsm + dsm - vegdsm[vegdsm == dsm] = 0 - vegdsm2 = vegdsm2 + dsm - vegdsm2[vegdsm2 == dsm] = 0 - - # % Bush separation - bush = np.logical_not((vegdsm2 * vegdsm)) * vegdsm - - # % major bug fixed 20141203 - svfbuveg = svf - (1.0 - svfveg) * (1.0 - transVeg) - else: - psi = leafon * 0.0 + 1.0 - svfbuveg = svf - bush = np.zeros([rows, cols]) - amaxvalue = 0 - - # Initialization of maps - Knight = np.zeros((rows, cols)) - Tgmap1 = np.zeros((rows, cols)) - Tgmap1E = np.zeros((rows, cols)) - Tgmap1S = np.zeros((rows, cols)) - Tgmap1W = np.zeros((rows, cols)) - Tgmap1N = np.zeros((rows, cols)) - - # Create building boolean raster from either land cover or height rasters - if demforbuild == 0: - buildings = np.copy(lcgrid) - buildings[buildings == 7] = 1 - buildings[buildings == 6] = 1 - buildings[buildings == 5] = 1 - buildings[buildings == 4] = 1 - buildings[buildings == 3] = 1 - buildings[buildings == 2] = 0 - else: - buildings = dsm - dem - buildings[buildings < 2.0] = 1.0 - buildings[buildings >= 2.0] = 0.0 - - if int(configDict["savebuild"]) == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/buildings.tif", - buildings, - ) - else: - common.save_raster( - configDict["output_dir"] + "/buildings.tif", - buildings, - dsm_transf, - dsm_crs, - ) - - # Import shadow matrices (Anisotropic sky) - anisotropic_sky = int(configDict["aniso"]) - if anisotropic_sky == 1: # UseAniso - data = np.load(configDict["input_aniso"]) - shmat = data["shadowmat"] - vegshmat = data["vegshadowmat"] - vbshvegshmat = data["vbshmat"] - if usevegdem == 1: - diffsh = np.zeros((rows, cols, shmat.shape[2])) - for i in range(0, shmat.shape[2]): - # changes in psi not implemented yet - diffsh[:, :, i] = shmat[:, :, i] - (1 - vegshmat[:, :, i]) * ( - 1 - transVeg - ) - else: - diffsh = shmat - - # Estimate number of patches based on shadow matrices - if shmat.shape[2] == 145: - patch_option = 1 # patch_option = 1 # 145 patches - elif shmat.shape[2] == 153: - patch_option = 2 # patch_option = 2 # 153 patches - elif shmat.shape[2] == 306: - patch_option = 3 # patch_option = 3 # 306 patches - elif shmat.shape[2] == 612: - patch_option = 4 # patch_option = 4 # 612 patches - - # asvf to calculate sunlit and shaded patches - asvf = np.arccos(np.sqrt(svf)) - - # Empty array for steradians - steradians = np.zeros((shmat.shape[2])) - else: - # anisotropic_sky = 0 - diffsh = None - shmat = None - vegshmat = None - vbshvegshmat = None - asvf = None - patch_option = 0 - steradians = 0 - - # % Ts parameterisation maps - if landcover == 1.0: - # Get land cover properties for Tg wave (land cover scheme based on - # Bogren et al. 2000, explained in Lindberg et al., 2008 and Lindberg, - # Onomura & Grimmond, 2016) - [ - TgK, - Tstart, - alb_grid, - emis_grid, - TgK_wall, - Tstart_wall, - TmaxLST, - TmaxLST_wall, - ] = Tgmaps_v1(lcgrid.copy(), param) - else: - TgK = Knight + param["Ts_deg"]["Value"]["Cobble_stone_2014a"] - Tstart = Knight - param["Tstart"]["Value"]["Cobble_stone_2014a"] - TmaxLST = param["TmaxLST"]["Value"]["Cobble_stone_2014a"] - alb_grid = ( - Knight - + param["Albedo"]["Effective"]["Value"]["Cobble_stone_2014a"] - ) - emis_grid = Knight + param["Emissivity"]["Value"]["Cobble_stone_2014a"] - TgK_wall = param["Ts_deg"]["Value"]["Walls"] - Tstart_wall = param["Tstart"]["Value"]["Walls"] - TmaxLST_wall = param["TmaxLST"]["Value"]["Walls"] - - # Import data for wall temperature parameterization TODO: fix for - # standalone - wallScheme = int(configDict["wallscheme"]) - if wallScheme == 1: - wallData = np.load(configDict["input_wall"]) - voxelMaps = wallData["voxelId"] - voxelTable = wallData["voxelTable"] - # Get wall type from standalone - if standAlone == 1: - wall_type_standalone = { - "Brick_wall": "100", - "Concrete_wall": "101", - "Wood_wall": "102", - } - wall_type = wall_type_standalone[configDict["walltype"]] - else: - # Get wall type set in GUI - wall_type = str(configDict["walltype"]) - # Calculate wall height for wall scheme, i.e. include corners (thicker - # walls) - walls_scheme = wa.findwalls_sp( - dsm, 2, np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]]) - ) - # Calculate wall aspect for wall scheme, i.e. include corners (thicker - # walls) - dirwalls_scheme = wa.filter1Goodwin_as_aspect_v3( - walls_scheme.copy(), scale, dsm, feedback, 100.0 / 180.0 - ) - - # Used in wall temperature parameterization scheme - first_timestep = ( - pd.to_datetime(YYYY[0][0], format="%Y") - + pd.to_timedelta(DOY[0] - 1, unit="d") - + pd.to_timedelta(hours[0], unit="h") - + pd.to_timedelta(minu[0], unit="m") - ) - second_timestep = ( - pd.to_datetime(YYYY[0][1], format="%Y") - + pd.to_timedelta(DOY[1] - 1, unit="d") - + pd.to_timedelta(hours[1], unit="h") - + pd.to_timedelta(minu[1], unit="m") - ) - - timeStep = (second_timestep - first_timestep).seconds - - # Load voxelTable as Pandas DataFrame - voxelTable, dirwalls_scheme = load_walls( - voxelTable, - param, - wall_type, - dirwalls_scheme, - Ta[0], - timeStep, - albedo_b, - ewall, - alb_grid, - landcover, - lcgrid, - dsm, - ) - - # Use wall of interest - woi_file = configDict["woi_file"] - - if woi_file: - # (dsm_minx, dsm_x_size, dsm_x_rotation, dsm_miny, dsm_y_rotation, dsm_y_size) = gdal_dsm.GetGeoTransform() #TODO: fix for standalone - if standAlone == 0: - # self.parameterAsStrings(parameters, self.WOI_FIELD, context) - woi_field = configDict["woi_field"] - woisxy, woiname = pointOfInterest( - configDict["woi_file"], woi_field, scale, gdal_dsm - ) - else: - pois_gdf = gpd.read_file(configDict["poi_file"]) - numfeat = pois_gdf.shape[0] - poisxy = np.zeros((numfeat, 3)) - 999 - for idx, row in pois_gdf.iterrows(): - # TODO: This produce different result since no standalone - # round coordinates - y, x = rowcol( - dsm_transf, - row["geometry"].centroid.x, - row["geometry"].centroid.y, - ) - poiname.append(row[configDict["poi_field"]]) - poisxy[idx, 0] = idx - poisxy[idx, 1] = x - poisxy[idx, 2] = y - else: - woisxy = None - # Create pandas datetime object to be used when createing an xarray - # DataSet where wall temperatures/radiation is stored and eventually - # saved as a NetCDf - if configDict["wallnetcdf"] == 1: - met_for_xarray = ( - pd.to_datetime(YYYY[0][:], format="%Y") - + pd.to_timedelta(DOY - 1, unit="d") - + pd.to_timedelta(hours, unit="h") - + pd.to_timedelta(minu, unit="m") - ) - else: - wallScheme = 0 - voxelMaps = 0 - voxelTable = 0 - timeStep = 0 - # thermal_effusivity = 0 - walls_scheme = np.ones((rows, cols)) * 10.0 - dirwalls_scheme = np.ones((rows, cols)) * 10.0 - - # Initialisation of time related variables - if Ta.__len__() == 1: - timestepdec = 0 - else: - timestepdec = dectime[1] - dectime[0] - timeadd = 0.0 - # timeaddE = 0. - # timeaddS = 0. - # timeaddW = 0. - # timeaddN = 0. - firstdaytime = 1.0 - - # Save hemispheric image - if anisotropic_sky == 1: - if not poisxy is None: - patch_characteristics = hemispheric_image( - poisxy, shmat, vegshmat, vbshvegshmat, voxelMaps, wallScheme - ) - - # If metfile starts at night - CI = 1.0 - - # Main loop - tmrtplot = np.zeros((rows, cols)) - TgOut1 = np.zeros((rows, cols)) - - # Initiate array for I0 values - if np.unique(DOY).shape[0] > 1: - unique_days = np.unique(DOY) - first_unique_day = DOY[DOY == unique_days[0]] - I0_array = np.zeros((first_unique_day.shape[0])) - else: - first_unique_day = DOY.copy() - I0_array = np.zeros((DOY.shape[0])) - - if standAlone == 1: - progress = tqdm(total=Ta.__len__()) - - for i in np.arange(0, Ta.__len__()): - if feedback is not None: - # move progressbar forward - feedback.setProgress(int(i * (100.0 / Ta.__len__()))) - if feedback.isCanceled(): - feedback.setProgressText("Calculation cancelled") - break - else: - progress.update(1) - - # Daily water body temperature - if landcover == 1: - if ((dectime[i] - np.floor(dectime[i]))) == 0 or (i == 0): - Twater = np.mean(Ta[jday[0] == np.floor(dectime[i])]) - # Nocturnal cloudfraction from Offerle et al. 2003 - if (dectime[i] - np.floor(dectime[i])) == 0: - daylines = np.where(np.floor(dectime) == dectime[i]) - if daylines.__len__() > 1: - alt = altitude[0][daylines] - alt2 = np.where(alt > 1) - rise = alt2[0][0] - [_, CI, _, _, _] = clearnessindex_2013b( - zen[0, i + rise + 1], - jday[0, i + rise + 1], - Ta[i + rise + 1], - RH[i + rise + 1] / 100.0, - radG[i + rise + 1], - location, - P[i + rise + 1], - ) - if (CI > 1.0) or (CI == np.inf): - CI = 1.0 - else: - CI = 1.0 - - # Only if Kdir is derived from horizontal global shortwave and horizontal diffuse shortwave - # if altitude[0][i] > 0: - # radI[i] = radI[i]/np.sin(altitude[0][i] * np.pi/180) - # else: - # radG[i] = 0. - # radD[i] = 0. - # radI[i] = 0. - - ( - 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, - radIout, - radDout, - Lside, - Lsky_patch_characteristics, - CI_Tg, - CI_TgG, - KsideD, - dRad, - Kside, - steradians, - voxelTable, - ) = so.Solweig_2025a_calc( - i, - dsm, - scale, - rows, - cols, - svf, - svfN, - svfW, - svfE, - svfS, - svfveg, - svfNveg, - svfEveg, - svfSveg, - svfWveg, - svfaveg, - svfEaveg, - svfSaveg, - svfWaveg, - svfNaveg, - vegdsm, - vegdsm2, - albedo_b, - absK, - absL, - ewall, - Fside, - Fup, - Fcyl, - altitude[0][i], - azimuth[0][i], - zen[0][i], - jday[0][i], - usevegdem, - onlyglobal, - buildings, - location, - psi[0][i], - landcover, - lcgrid, - dectime[i], - altmax[0][i], - wallaspect, - wallheight, - cyl, - elvis, - Ta[i], - RH[i], - radG[i], - radD[i], - radI[i], - P[i], - 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[i], - wallScheme, - timeStep, - steradians, - walls_scheme, - dirwalls_scheme, - ) - - # Save I0 for I0 vs. Kdown output plot to check if UTC is off - # if i == (first_unique_day.shape[0] - 1): - # # Output I0 vs. Kglobal plot - # radG_for_plot = radG[DOY == first_unique_day[0]] - # # hours_for_plot = hours[DOY == first_unique_day[0]] - # dectime_for_plot = dectime[DOY == first_unique_day[0]] - # fig, ax = plt.subplots() - # ax.plot(dectime_for_plot, I0_array, label='I0') - # ax.plot(dectime_for_plot, radG_for_plot, label='Kglobal') - # ax.set_ylabel('Shortwave radiation [$Wm^{-2}$]') - # ax.set_xlabel('Decimal time') - # ax.set_title('UTC' + str(configDict['utc'])) - # ax.legend() - # fig.savefig(configDict['output_dir'] + '/metCheck.png', dpi=150) - # elif i < (first_unique_day.shape[0] - 1): - # I0_array[i] = I0 - - # Save I0 for I0 vs. Kdown output plot to check if UTC is off - if i < first_unique_day.shape[0]: - I0_array[i] = I0 - elif i == first_unique_day.shape[0]: - # Output I0 vs. Kglobal plot - radG_for_plot = radG[DOY == first_unique_day[0]] - # hours_for_plot = hours[DOY == first_unique_day[0]] - dectime_for_plot = dectime[DOY == first_unique_day[0]] - fig, ax = plt.subplots() - ax.plot(dectime_for_plot, I0_array, label="I0") - ax.plot(dectime_for_plot, radG_for_plot, label="Kglobal") - ax.set_ylabel("Shortwave radiation [$Wm^{-2}$]") - ax.set_xlabel("Decimal time") - ax.set_title("UTC" + str(configDict["utc"])) - ax.legend() - fig.savefig(configDict["output_dir"] + "/metCheck.png", dpi=150) - - tmrtplot = tmrtplot + Tmrt - - if altitude[0][i] > 0: - w = "D" - else: - w = "N" - - # Write to POIs - if not poisxy is None: - for k in range(0, poisxy.shape[0]): - poi_save = np.zeros((1, 41)) - poi_save[0, 0] = YYYY[0][i] - poi_save[0, 1] = jday[0][i] - poi_save[0, 2] = hours[i] - poi_save[0, 3] = minu[i] - poi_save[0, 4] = dectime[i] - poi_save[0, 5] = altitude[0][i] - poi_save[0, 6] = azimuth[0][i] - poi_save[0, 7] = radIout - poi_save[0, 8] = radDout - poi_save[0, 9] = radG[i] - poi_save[0, 10] = Kdown[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 11] = Kup[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 12] = Keast[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 13] = Ksouth[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 14] = Kwest[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 15] = Knorth[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 16] = Ldown[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 17] = Lup[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 18] = Least[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 19] = Lsouth[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 20] = Lwest[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 21] = Lnorth[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 22] = Ta[i] - poi_save[0, 23] = TgOut[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 24] = RH[i] - poi_save[0, 25] = esky - poi_save[0, 26] = Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 27] = I0 - poi_save[0, 28] = CI - poi_save[0, 29] = shadow[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 30] = svf[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 31] = svfbuveg[ - int(poisxy[k, 2]), int(poisxy[k, 1]) - ] - poi_save[0, 32] = KsideI[int(poisxy[k, 2]), int(poisxy[k, 1])] - # Recalculating wind speed based on powerlaw - WsPET = (1.1 / sensorheight) ** 0.2 * Ws[i] - WsUTCI = (10.0 / sensorheight) ** 0.2 * Ws[i] - resultPET = p._PET( - Ta[i], - RH[i], - Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])], - WsPET, - mbody, - age, - ht, - activity, - clo, - sex, - ) - poi_save[0, 33] = resultPET - resultUTCI = utci.utci_calculator( - Ta[i], - RH[i], - Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])], - WsUTCI, - ) - poi_save[0, 34] = resultUTCI - poi_save[0, 35] = CI_Tg - poi_save[0, 36] = CI_TgG - poi_save[0, 37] = KsideD[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 38] = Lside[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 39] = dRad[int(poisxy[k, 2]), int(poisxy[k, 1])] - poi_save[0, 40] = Kside[int(poisxy[k, 2]), int(poisxy[k, 1])] - data_out = ( - configDict["output_dir"] - + "/POI_" - + str(poiname[k]) - + ".txt" - ) - # f_handle = file(data_out, 'a') - f_handle = open(data_out, "ab") - np.savetxt(f_handle, poi_save, fmt=numformat) - f_handle.close() - - # If wall temperature parameterization scheme is in use - if ( - configDict["wallscheme"] == 1 - ): # folderWallScheme: TODO: Fix for standalone - # Store wall data for output - if not woisxy is None: - for k in range(0, woisxy.shape[0]): - temp_wall = voxelTable.loc[ - ( - (voxelTable["ypos"] == woisxy[k, 2]) - & (voxelTable["xpos"] == woisxy[k, 1]) - ), - "wallTemperature", - ].to_numpy() - K_in = voxelTable.loc[ - ( - (voxelTable["ypos"] == woisxy[k, 2]) - & (voxelTable["xpos"] == woisxy[k, 1]) - ), - "K_in", - ].to_numpy() - L_in = voxelTable.loc[ - ( - (voxelTable["ypos"] == woisxy[k, 2]) - & (voxelTable["xpos"] == woisxy[k, 1]) - ), - "L_in", - ].to_numpy() - wallShade = voxelTable.loc[ - ( - (voxelTable["ypos"] == woisxy[k, 2]) - & (voxelTable["xpos"] == woisxy[k, 1]) - ), - "wallShade", - ].to_numpy() - temp_all = np.concatenate( - [temp_wall, K_in, L_in, wallShade] - ) - # temp_all = np.concatenate([temp_wall]) - # wall_data = np.zeros((1, 7 + temp_wall.shape[0])) - wall_data = np.zeros((1, 7 + temp_all.shape[0])) - # Part of file name (wallid), i.e. WOI_wallid.txt - data_out = ( - configDict["output_dir"] - + "/WOI_" - + str(woiname[k]) - + ".txt" - ) - if i == 0: - # Output file header - # header = 'yyyy id it imin dectime Ta SVF Ts' - header = ( - "yyyy id it imin dectime Ta SVF" - + " Ts" * temp_wall.shape[0] - + " Kin" * K_in.shape[0] - + " Lin" * L_in.shape[0] - + " shade" * wallShade.shape[0] - ) - # Part of file name (wallid), i.e. WOI_wallid.txt - # woiname = voxelTable.loc[((voxelTable['ypos'] == woisxy[k, 2]) & (voxelTable['xpos'] == woisxy[k, 1])), 'wallId'].to_numpy()[0] - woi_save = [] # - np.savetxt( - data_out, - woi_save, - delimiter=" ", - header=header, - comments="", - ) - # Fill wall_data with variables - wall_data[0, 0] = YYYY[0][i] - wall_data[0, 1] = jday[0][i] - wall_data[0, 2] = hours[i] - wall_data[0, 3] = minu[i] - wall_data[0, 4] = dectime[i] - wall_data[0, 5] = Ta[i] - wall_data[0, 6] = svf[int(woisxy[k, 2]), int(woisxy[k, 1])] - wall_data[0, 7:] = temp_all - - # Num format for output file data - woi_numformat = ( - "%d %d %d %d %.5f %.2f %.2f" - + " %.2f" * temp_all.shape[0] - ) - # Open file, add data, save - f_handle = open(data_out, "ab") - np.savetxt(f_handle, wall_data, fmt=woi_numformat) - f_handle.close() - - # Save wall temperature/radiation as NetCDF TODO: fix for - # standAlone? - if configDict["wallnetcdf"] == 1: # wallNetCDF: - netcdf_output = configDict["output_dir"] + "/walls.nc" - walls_as_netcdf( - voxelTable, - rows, - cols, - met_for_xarray, - i, - dsm, - configDict["filepath_dsm"], - netcdf_output, - ) - - if hours[i] < 10: - XH = "0" - else: - XH = "" - if minu[i] < 10: - XM = "0" - else: - XM = "" - - time_code = ( - str(int(YYYY[0, i])) - + "_" - + str(int(DOY[i])) - + "_" - + XH - + str(int(hours[i])) - + XM - + str(int(minu[i])) - + w - ) - - if configDict["outputtmrt"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Tmrt_" + time_code + ".tif", - Tmrt, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Tmrt_" + time_code + ".tif", - Tmrt, - dsm_transf, - dsm_crs, - ) - if configDict["outputkup"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Kup_" + time_code + ".tif", - Kup, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Kup_" + time_code + ".tif", - Kup, - dsm_transf, - dsm_crs, - ) - if configDict["outputkdown"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Kdown_" + time_code + ".tif", - Kdown, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Kdown_" + time_code + ".tif", - Kdown, - dsm_transf, - dsm_crs, - ) - if configDict["outputlup"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Lup_" + time_code + ".tif", - Lup, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Lup_" + time_code + ".tif", - Lup, - dsm_transf, - dsm_crs, - ) - if configDict["outputldown"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Ldown_" + time_code + ".tif", - Ldown, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Ldown_" + time_code + ".tif", - Ldown, - dsm_transf, - dsm_crs, - ) - if configDict["outputsh"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Shadow_" + time_code + ".tif", - shadow, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Shadow_" + time_code + ".tif", - shadow, - dsm_transf, - dsm_crs, - ) - - if configDict["outputkdiff"] == 1: - if standAlone == 0: - saveraster( - gdal_dsm, - configDict["output_dir"] + "/Kdiff_" + time_code + ".tif", - dRad, - ) - else: - common.save_raster( - configDict["output_dir"] + "/Kdiff_" + time_code + ".tif", - dRad, - dsm_transf, - dsm_crs, - ) - - # Sky view image of patches - if (anisotropic_sky == 1) & (i == 0) & (not poisxy is None): - for k in range(poisxy.shape[0]): - Lsky_patch_characteristics[:, 2] = patch_characteristics[:, k] - skyviewimage_out = ( - configDict["output_dir"] - + "/POI_" - + str(poiname[k]) - + ".png" - ) - PolarBarPlot( - Lsky_patch_characteristics, - altitude[0][i], - azimuth[0][i], - "Hemisphere partitioning", - skyviewimage_out, - 0, - 5, - 0, - ) - - # Save files for Tree Planter - if configDict["outputtreeplanter"] == 1: # outputTreeplanter: - if feedback is not None: - feedback.setProgressText("Saving files for Tree Planter tool") - # Save DSM - copyfile( - configDict["filepath_dsm"], configDict["output_dir"] + "/DSM.tif" - ) - - # Save CDSM - if usevegdem == 1: - copyfile( - configDict["filepath_cdsm"], - configDict["output_dir"] + "/CDSM.tif", - ) - - albedo_g = param["Albedo"]["Effective"]["Value"]["Cobble_stone_2014a"] - eground = param["Emissivity"]["Value"]["Cobble_stone_2014a"] - - # Saving settings from SOLWEIG for SOLWEIG1D in TreePlanter - settingsHeader = "UTC, posture, onlyglobal, landcover, anisotropic, cylinder, albedo_walls, albedo_ground, emissivity_walls, emissivity_ground, absK, absL, elevation, patch_option" - settingsFmt = ( - "%i", - "%i", - "%i", - "%i", - "%i", - "%i", - "%1.2f", - "%1.2f", - "%1.2f", - "%1.2f", - "%1.2f", - "%1.2f", - "%1.2f", - "%i", - ) - settingsData = np.array( - [ - [ - int(configDict["utc"]), - pos, - onlyglobal, - landcover, - anisotropic_sky, - cyl, - albedo_b, - albedo_g, - ewall, - eground, - absK, - absL, - alt, - patch_option, - ] - ] - ) - # print(settingsData) - np.savetxt( - configDict["output_dir"] + "/treeplantersettings.txt", - settingsData, - fmt=settingsFmt, - header=settingsHeader, - delimiter=" ", - ) - - # Copying met file for SpatialTC - copyfile( - configDict["input_met"], configDict["output_dir"] + "/metforcing.txt" - ) - - tmrtplot = ( - tmrtplot / Ta.__len__() - ) # fix average Tmrt instead of sum, 20191022 - if standAlone == 0: - saveraster( - gdal_dsm, configDict["output_dir"] + "/Tmrt_average.tif", tmrtplot - ) - else: - common.save_raster( - configDict["output_dir"] + "/Tmrt_average.tif", - tmrtplot, - dsm_transf, - dsm_crs, - ) +# This is the main function of the SOLWEIG model +# 2025-Mar-21 +# Fredrik Lindberg, fredrikl@gvc.gu.se +# Goteborg Urban Climate Group +# Gothenburg University + +# sommon imports +from __future__ import absolute_import +from ...util.umep_solweig_export_component import read_solweig_config +from ...util.SEBESOLWEIGCommonFiles.Solweig_v2015_metdata_noload import ( + Solweig_2015a_metdata_noload, +) +from ...util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import ( + clearnessindex_2013b, +) + +from ...functions.SOLWEIGpython import ( + Solweig_2026a_calc_forprocessing as so, +) + +# from ...functions.SOLWEIGpython import WriteMetadataSOLWEIG # Not needed anymore? +from ...functions.SOLWEIGpython import PET_calculations as p +from ...functions.SOLWEIGpython import UTCI_calculations as utci +from ...functions.SOLWEIGpython.CirclePlotBar import PolarBarPlot +from ...functions.SOLWEIGpython.wall_surface_temperature import load_walls +from ...functions.SOLWEIGpython.wallOfInterest import pointOfInterest +from ...functions.SOLWEIGpython.patch_characteristics import hemispheric_image +from ...functions.SOLWEIGpython.wallsAsNetCDF import walls_as_netcdf +from ...functions.SOLWEIGpython.Tgmaps_v1 import Tgmaps_v1 +from ...functions import wallalgorithms as wa +from ...functions.SOLWEIGpython.ground_surface import ( + initiate_groundScheme, +) +import numpy as np +import json +import zipfile +import pandas as pd +import matplotlib.pylab as plt +from shutil import copyfile + +# imports from osgeo/qgis dependency +try: + from osgeo import gdal + from osgeo.gdalconst import * + from ...util.misc import saveraster, xy2latlon_fromraster + from qgis.core import QgsVectorLayer, QgsRasterLayer +except: + pass + +# imports for standalone +try: + from umep import common + from rasterio.transform import xy, rowcol + import pyproj + from tqdm import tqdm + import geopandas as gpd +except: + pass + + +# import numpy as np +# from .daylen import daylen +# from ...util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import clearnessindex_2013b +# from ...util.SEBESOLWEIGCommonFiles.diffusefraction import diffusefraction +# 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 .cylindric_wedge import cylindric_wedge +# from .TsWaveDelay_2015a import TsWaveDelay_2015a +# from .Kup_veg_2015a import Kup_veg_2015a +# # 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 + + +def solweig_run(configPath, feedback): + """ + Input: + configPath : config file including geodata paths and settings. + feedback : To communicate with qgis gui. Set to None if standalone + """ + + # Load config file + configDict = read_solweig_config(configPath) + + # Load parameters settings for SOLWEIG + with open(configDict["para_json_path"], "r") as jsn: + param = json.load(jsn) + + standAlone = int(configDict["standalone"]) + + # reading variables from config and parameters that is not yet presented + cyl = int(configDict["cyl"]) + albedo_b = param["Albedo"]["Effective"]["Value"]["Walls"] + ewall = param["Emissivity"]["Value"]["Walls"] + onlyglobal = int(configDict["onlyglobal"]) + elvis = 0.0 + absK = param["Tmrt_params"]["Value"]["absK"] + absL = param["Tmrt_params"]["Value"]["absL"] + + # Load DSM + if standAlone == 1: + dsm, dsm_transf, dsm_crs = common.load_raster( + configDict["filepath_dsm"], bbox=None + ) + scale = 1 / dsm_transf.a + # dsm_height, dsm_width = dsm.shape # y rows by x cols + # y is flipped - so return max for lower row + minx, miny = xy(dsm_transf, dsm.shape[0], 0) + # Define the source and target CRS + source_crs = pyproj.CRS(dsm_crs) + target_crs = pyproj.CRS(4326) # WGS 84 + # Create a transformer object + transformer = pyproj.Transformer.from_crs( + source_crs, target_crs, always_xy=True + ) + # Perform the transformation + lon, lat = transformer.transform(minx, miny) + nd = -9999 # TODO: extract nodatavalue from rasterio + else: + # dsmlayer = QgsRasterLayer(configDict['filepath_dsm']) + dsm_wkt = QgsRasterLayer(configDict["filepath_dsm"]).crs().toWkt() + gdal_dsm = gdal.Open(configDict["filepath_dsm"]) + lat, lon, scale, minx, miny = xy2latlon_fromraster(dsm_wkt, gdal_dsm) + dsm = gdal_dsm.ReadAsArray().astype(float) + nd = gdal_dsm.GetRasterBand(1).GetNoDataValue() + + rows = dsm.shape[0] + cols = dsm.shape[1] + + # response to issue #85 + dsm[dsm == nd] = 0.0 + if dsm.min() < 0: + dsmraise = np.abs(dsm.min()) + dsm = dsm + dsmraise + else: + dsmraise = 0 + + alt = np.median(dsm) + if alt < 0: + alt = 3 + + # Vegetation + transVeg = param["Tree_settings"]["Value"]["Transmissivity"] + trunkratio = param["Tree_settings"]["Value"]["Trunk_ratio"] + usevegdem = int(configDict["usevegdem"]) + if usevegdem == 1: + if standAlone == 0: + vegdsm = ( + gdal.Open(configDict["filepath_cdsm"]) + .ReadAsArray() + .astype(float) + ) + else: + vegdsm, _, _ = common.load_raster( + configDict["filepath_cdsm"], bbox=None + ) + if configDict["filepath_tdsm"] != "": + if standAlone == 0: + vegdsm2 = ( + gdal.Open(configDict["filepath_tdsm"]) + .ReadAsArray() + .astype(float) + ) + else: + vegdsm2, _, _ = common.load_raster( + configDict["filepath_tdsm"], bbox=None + ) + else: + vegdsm2 = vegdsm * trunkratio + else: + vegdsm = 0 + vegdsm2 = 0 + + # Land cover + landcover = int(configDict["landcover"]) + if landcover == 1: + if standAlone == 0: + lcgrid = ( + gdal.Open(configDict["filepath_lc"]) + .ReadAsArray() + .astype(float) + ) + else: + lcgrid, _, _ = common.load_raster( + configDict["filepath_lc"], bbox=None + ) + else: + lcgrid = np.ones_like(dsm) + + # DEM for buildings #TODO: fix nodata in standalone + demforbuild = int(configDict["demforbuild"]) + if demforbuild == 1: + if standAlone == 0: + gdal_dem = gdal.Open( + configDict["filepath_dem"] + ) # .ReadAsArray().astype(float) + dem = gdal_dem.ReadAsArray().astype(float) + nd = gdal_dem.GetRasterBand(1).GetNoDataValue() + else: + dem, _, _ = common.load_raster( + configDict["filepath_dem"], bbox=None + ) + nd = -9999 # TODO: standAlone nd exposure + + # response to issue and #230 + dem[dem == nd] = 0.0 + if dem.min() < 0: + demraise = np.abs(dem.min()) + dem = dem + demraise + else: + demraise = 0 + + # SVF + zip = zipfile.ZipFile(configDict["input_svf"], "r") + zip.extractall(configDict["working_dir"]) + zip.close() + + if standAlone == 0: + svf = ( + gdal.Open(configDict["working_dir"] + "/svf.tif") + .ReadAsArray() + .astype(float) + ) + svfN = ( + gdal.Open(configDict["working_dir"] + "/svfN.tif") + .ReadAsArray() + .astype(float) + ) + svfS = ( + gdal.Open(configDict["working_dir"] + "/svfS.tif") + .ReadAsArray() + .astype(float) + ) + svfE = ( + gdal.Open(configDict["working_dir"] + "/svfE.tif") + .ReadAsArray() + .astype(float) + ) + svfW = ( + gdal.Open(configDict["working_dir"] + "/svfW.tif") + .ReadAsArray() + .astype(float) + ) + else: + 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 + ) + + if usevegdem == 1: + if standAlone == 0: + svfveg = ( + gdal.Open(configDict["working_dir"] + "/svfveg.tif") + .ReadAsArray() + .astype(float) + ) + svfNveg = ( + gdal.Open(configDict["working_dir"] + "/svfNveg.tif") + .ReadAsArray() + .astype(float) + ) + svfSveg = ( + gdal.Open(configDict["working_dir"] + "/svfSveg.tif") + .ReadAsArray() + .astype(float) + ) + svfEveg = ( + gdal.Open(configDict["working_dir"] + "/svfEveg.tif") + .ReadAsArray() + .astype(float) + ) + svfWveg = ( + gdal.Open(configDict["working_dir"] + "/svfWveg.tif") + .ReadAsArray() + .astype(float) + ) + + svfaveg = ( + gdal.Open(configDict["working_dir"] + "/svfaveg.tif") + .ReadAsArray() + .astype(float) + ) + svfNaveg = ( + gdal.Open(configDict["working_dir"] + "/svfNaveg.tif") + .ReadAsArray() + .astype(float) + ) + svfSaveg = ( + gdal.Open(configDict["working_dir"] + "/svfSaveg.tif") + .ReadAsArray() + .astype(float) + ) + svfEaveg = ( + gdal.Open(configDict["working_dir"] + "/svfEaveg.tif") + .ReadAsArray() + .astype(float) + ) + svfWaveg = ( + gdal.Open(configDict["working_dir"] + "/svfWaveg.tif") + .ReadAsArray() + .astype(float) + ) + else: + 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 + ) + else: + svfveg = np.ones((rows, cols)) + svfNveg = np.ones((rows, cols)) + svfSveg = np.ones((rows, cols)) + svfEveg = np.ones((rows, cols)) + svfWveg = np.ones((rows, cols)) + svfaveg = np.ones((rows, cols)) + svfNaveg = np.ones((rows, cols)) + svfSaveg = np.ones((rows, cols)) + svfEaveg = np.ones((rows, cols)) + svfWaveg = np.ones((rows, cols)) + + tmp = svf + svfveg - 1.0 + tmp[tmp < 0.0] = 0.0 + # %matlab crazyness around 0 + svfalfa = np.arcsin(np.exp((np.log((1.0 - tmp)) / 2.0))) + + if standAlone == 0: + wallheight = ( + gdal.Open(configDict["filepath_wh"]).ReadAsArray().astype(float) + ) + wallaspect = ( + gdal.Open(configDict["filepath_wa"]).ReadAsArray().astype(float) + ) + else: + wallheight, _, _ = common.load_raster( + configDict["filepath_wh"], bbox=None + ) + wallaspect, _, _ = common.load_raster( + configDict["filepath_wa"], bbox=None + ) + + # Metdata + headernum = 1 + delim = " " + Twater = [] + + metdata = np.loadtxt( + configDict["input_met"], skiprows=headernum, delimiter=delim + ) + + location = {"longitude": lon, "latitude": lat, "altitude": alt} + YYYY, altitude, azimuth, zen, jday, leafon, dectime, altmax = ( + Solweig_2015a_metdata_noload(metdata, location, int(configDict["utc"])) + ) + + DOY = metdata[:, 1] + hours = metdata[:, 2] + minu = metdata[:, 3] + Ta = metdata[:, 11] + RH = metdata[:, 10] + radG = metdata[:, 14] + radD = metdata[:, 21] + radI = metdata[:, 22] + P = metdata[:, 12] + Ws = metdata[:, 9] + + # POIs check + if configDict["poi_file"] != "": # usePOI: + header = ( + "yyyy id it imin dectime altitude azimuth kdir kdiff kglobal kdown kup keast ksouth " + "kwest knorth ldown lup least lsouth lwest lnorth Ta Tg RH Esky Tmrt " + "I0 CI Shadow SVF_b SVF_bv KsideI PET UTCI CI_TgG KsideD Lside diffDown Kside " + ) + # poiname = [] + poi_field = configDict[ + "poi_field" + ] # self.parameterAsString(parameters, self.POI_FIELD, context) + if standAlone == 0: + # vlayer = QgsVectorLayer(configDict['poi_file'], 'point', 'ogr') + # idx = vlayer.fields().indexFromName(poi_field) + # numfeat = vlayer.featureCount() + # poisxy = np.zeros((numfeat, 3)) - 999 + # ind = 0 + # for f in vlayer.getFeatures(): # looping through each POI + # y = f.geometry().centroid().asPoint().y() + # x = f.geometry().centroid().asPoint().x() + # poiname.append(f.attributes()[idx]) + # poisxy[ind, 0] = ind + # poisxy[ind, 1] = np.round((x - minx) * scale) + # if miny >= 0: + # poisxy[ind, 2] = np.round((miny + rows * (1. / scale) - y) * scale) + # else: + # poisxy[ind, 2] = np.round((miny + rows * (1. / scale) - y) * scale) + # ind += 1 + + poi_field = configDict[ + "woi_field" + ] # self.parameterAsStrings(parameters, self.WOI_FIELD, context) + poisxy, poiname = pointOfInterest( + configDict["poi_file"], poi_field, scale, gdal_dsm + ) + + else: + pois_gdf = gpd.read_file(configDict["poi_file"]) + numfeat = pois_gdf.shape[0] + poisxy = np.zeros((numfeat, 3)) - 999 + for idx, row in pois_gdf.iterrows(): + y, x = rowcol( + dsm_transf, + row["geometry"].centroid.x, + row["geometry"].centroid.y, + ) # TODO: This produce different result since no standalone round coordinates + poiname.append(row[configDict["poi_field"]]) + poisxy[idx, 0] = idx + poisxy[idx, 1] = x + poisxy[idx, 2] = y + + for k in range(0, poisxy.shape[0]): + poi_save = [] # np.zeros((1, 33)) + data_out = ( + configDict["output_dir"] + "/POI_" + str(poiname[k]) + ".txt" + ) + np.savetxt( + data_out, poi_save, delimiter=" ", header=header, comments="" + ) + # print(poisxy) + # Num format for POI output + numformat = "%d %d %d %d %.5f " + "%.2f " * 35 + + # Other PET variables + sensorheight = param["Wind_Height"]["Value"]["magl"] + age = param["PET_settings"]["Value"]["Age"] + mbody = param["PET_settings"]["Value"]["Weight"] + ht = param["PET_settings"]["Value"]["Height"] + clo = param["PET_settings"]["Value"]["clo"] + activity = param["PET_settings"]["Value"]["Activity"] + sex = param["PET_settings"]["Value"]["Sex"] + else: + poisxy = None + + # Posture settings + if param["Tmrt_params"]["Value"]["posture"] == "Standing": + Fside = param["Posture"]["Standing"]["Value"]["Fside"] + Fup = param["Posture"]["Standing"]["Value"]["Fup"] + height = param["Posture"]["Standing"]["Value"]["height"] + Fcyl = param["Posture"]["Standing"]["Value"]["Fcyl"] + pos = 1 + else: + Fside = param["Posture"]["Sitting"]["Value"]["Fside"] + Fup = param["Posture"]["Sitting"]["Value"]["Fup"] + height = param["Posture"]["Sitting"]["Value"]["height"] + Fcyl = param["Posture"]["Sitting"]["Value"]["Fcyl"] + pos = 0 + + # Radiative surface influence, Rule of thumb by Schmid et al. (1990). + first = np.round(height) + if first == 0.0: + first = 1.0 + second = np.round((height * 20.0)) + + if usevegdem == 1: + # Conifer or deciduous + if configDict["conifer_bool"]: + leafon = np.ones((1, DOY.shape[0])) + else: + leafon = np.zeros((1, DOY.shape[0])) + if ( + param["Tree_settings"]["Value"]["First_day_leaf"] + > param["Tree_settings"]["Value"]["Last_day_leaf"] + ): + leaf_bool = ( + DOY > param["Tree_settings"]["Value"]["First_day_leaf"] + ) | (DOY < param["Tree_settings"]["Value"]["Last_day_leaf"]) + else: + leaf_bool = ( + DOY > param["Tree_settings"]["Value"]["First_day_leaf"] + ) & (DOY < param["Tree_settings"]["Value"]["Last_day_leaf"]) + leafon[0, leaf_bool] = 1 + + # % Vegetation transmittivity of shortwave radiation + psi = leafon * transVeg + psi[leafon == 0] = 0.5 + # amaxvalue + vegmax = vegdsm.max() + amaxvalue = dsm.max() - dsm.min() + amaxvalue = np.maximum(amaxvalue, vegmax) + + # Elevation vegdsms if buildingDEM includes ground heights + vegdsm = vegdsm + dsm + vegdsm[vegdsm == dsm] = 0 + vegdsm2 = vegdsm2 + dsm + vegdsm2[vegdsm2 == dsm] = 0 + + # % Bush separation + bush = np.logical_not((vegdsm2 * vegdsm)) * vegdsm + + svfbuveg = svf - (1.0 - svfveg) * ( + 1.0 - transVeg + ) # % major bug fixed 20141203 + else: + psi = leafon * 0.0 + 1.0 + svfbuveg = svf + bush = np.zeros([rows, cols]) + amaxvalue = 0 + + # Initialization of maps + Knight = np.zeros((rows, cols)) + Tgmap1 = np.zeros((rows, cols)) + Tgmap1E = np.zeros((rows, cols)) + Tgmap1S = np.zeros((rows, cols)) + Tgmap1W = np.zeros((rows, cols)) + Tgmap1N = np.zeros((rows, cols)) + + # Create building boolean raster from either land cover or height rasters + if demforbuild == 0: + buildings = np.copy(lcgrid) + buildings[buildings == 7] = 1 + buildings[buildings == 6] = 1 + buildings[buildings == 5] = 1 + buildings[buildings == 4] = 1 + buildings[buildings == 3] = 1 + buildings[buildings == 2] = 0 + else: + buildings = dsm - dem + buildings[buildings < 2.0] = 1.0 + buildings[buildings >= 2.0] = 0.0 + + if int(configDict["savebuild"]) == 1: + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/buildings.tif", + buildings, + ) + else: + common.save_raster( + configDict["output_dir"] + "/buildings.tif", + buildings, + dsm_transf, + dsm_crs, + ) + + # Import shadow matrices (Anisotropic sky) + anisotropic_sky = int(configDict["aniso"]) + if anisotropic_sky == 1: # UseAniso + data = np.load(configDict["input_aniso"]) + shmat = data["shadowmat"] + vegshmat = data["vegshadowmat"] + vbshvegshmat = data["vbshmat"] + if usevegdem == 1: + diffsh = np.zeros((rows, cols, shmat.shape[2])) + for i in range(0, shmat.shape[2]): + diffsh[:, :, i] = shmat[:, :, i] - (1 - vegshmat[:, :, i]) * ( + 1 - transVeg + ) # changes in psi not implemented yet + else: + diffsh = shmat + + # Estimate number of patches based on shadow matrices + if shmat.shape[2] == 145: + patch_option = 1 # patch_option = 1 # 145 patches + elif shmat.shape[2] == 153: + patch_option = 2 # patch_option = 2 # 153 patches + elif shmat.shape[2] == 306: + patch_option = 3 # patch_option = 3 # 306 patches + elif shmat.shape[2] == 612: + patch_option = 4 # patch_option = 4 # 612 patches + + # asvf to calculate sunlit and shaded patches + asvf = np.arccos(np.sqrt(svf)) + + # Empty array for steradians + steradians = np.zeros((shmat.shape[2])) + else: + # anisotropic_sky = 0 + diffsh = None + shmat = None + vegshmat = None + vbshvegshmat = None + asvf = None + patch_option = 0 + steradians = 0 + shadow = np.zeros_like(dsm) + + # % Ts parameterisation maps + if landcover == 1.0: + # Get land cover properties for Tg wave (land cover scheme based on Bogren et al. 2000, explained in Lindberg et al., 2008 and Lindberg, Onomura & Grimmond, 2016) + [ + TgK, + Tstart, + alb_grid, + emis_grid, + TgK_wall, + Tstart_wall, + TmaxLST, + TmaxLST_wall, + ] = Tgmaps_v1(lcgrid.copy(), param) + else: + TgK = Knight + param["Ts_deg"]["Value"]["Cobble_stone_2014a"] + Tstart = Knight - param["Tstart"]["Value"]["Cobble_stone_2014a"] + TmaxLST = param["TmaxLST"]["Value"]["Cobble_stone_2014a"] + alb_grid = ( + Knight + + param["Albedo"]["Effective"]["Value"]["Cobble_stone_2014a"] + ) + emis_grid = Knight + param["Emissivity"]["Value"]["Cobble_stone_2014a"] + TgK_wall = param["Ts_deg"]["Value"]["Walls"] + Tstart_wall = param["Tstart"]["Value"]["Walls"] + TmaxLST_wall = param["TmaxLST"]["Value"]["Walls"] + + # Parameterization for the 2026 ground scheme + groundSurface = int(configDict["groundmodel"]) + if groundSurface == 1: + # Initiate the maps if the surface temperature is available + if configDict["input_surf"] != "": + surfData = pd.read_csv(configDict["input_surf"]) + Tg = surfData["Tg"] + Tm = np.mean(surfData["Tg"]) + ( + _, + _, + Rn, + Rn_past, + G, + cap_grid, + diff_grid, + a1_grid, + a2_grid, + a3_grid, + ) = initiate_groundScheme( + lcgrid.copy(), param, DOY[0], Ta, location + ) + else: + ( + Tg, + Tm, + Rn, + Rn_past, + G, + cap_grid, + diff_grid, + a1_grid, + a2_grid, + a3_grid, + ) = initiate_groundScheme( + lcgrid.copy(), param, DOY[0], Ta, location + ) + else: + pass + + # Replace the ground view factors with integration of solid angles + outgoingLW = int(configDict["outgoingLW"]) + + # Import data for wall temperature parameterization TODO: fix for standalone + wallScheme = int(configDict["wallscheme"]) + if wallScheme == 1: + wallData = np.load(configDict["input_wall"]) + voxelMaps = wallData["voxelId"] + voxelTable = wallData["voxelTable"] + # Get wall type from standalone + if standAlone == 1: + wall_type_standalone = { + "Brick_wall": "100", + "Concrete_wall": "101", + "Wood_wall": "102", + } + wall_type = wall_type_standalone[configDict["walltype"]] + else: + # Get wall type set in GUI + wall_type = configDict[ + "walltype" + ] # str(100 + int(self.parameterAsString(parameters, self.WALL_TYPE, context))) #TODO + + # Calculate wall height for wall scheme, i.e. include corners (thicker walls) + walls_scheme = wa.findwalls_sp( + dsm, 2, np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]]) + ) + # Calculate wall aspect for wall scheme, i.e. include corners (thicker walls) + dirwalls_scheme = wa.filter1Goodwin_as_aspect_v3( + walls_scheme.copy(), scale, dsm, feedback, 100.0 / 180.0 + ) + + # Used in wall temperature parameterization scheme + first_timestep = ( + pd.to_datetime(YYYY[0][0], format="%Y") + + pd.to_timedelta(DOY[0] - 1, unit="d") + + pd.to_timedelta(hours[0], unit="h") + + pd.to_timedelta(minu[0], unit="m") + ) + second_timestep = ( + pd.to_datetime(YYYY[0][1], format="%Y") + + pd.to_timedelta(DOY[1] - 1, unit="d") + + pd.to_timedelta(hours[1], unit="h") + + pd.to_timedelta(minu[1], unit="m") + ) + + timeStep = (second_timestep - first_timestep).seconds + + # Load voxelTable as Pandas DataFrame + voxelTable, dirwalls_scheme = load_walls( + voxelTable, + param, + wall_type, + dirwalls_scheme, + Ta[0], + timeStep, + albedo_b, + ewall, + alb_grid, + landcover, + lcgrid, + dsm, + ) + + # Use wall of interest + woi_file = configDict["woi_file"] + if woi_file: + # (dsm_minx, dsm_x_size, dsm_x_rotation, dsm_miny, dsm_y_rotation, dsm_y_size) = gdal_dsm.GetGeoTransform() #TODO: fix for standalone + if standAlone == 0: + woi_field = configDict[ + "woi_field" + ] # self.parameterAsStrings(parameters, self.WOI_FIELD, context) + woisxy, woiname = pointOfInterest( + configDict["woi_file"], woi_field, scale, gdal_dsm + ) + else: + pois_gdf = gpd.read_file(configDict["poi_file"]) + numfeat = pois_gdf.shape[0] + poisxy = np.zeros((numfeat, 3)) - 999 + for idx, row in pois_gdf.iterrows(): + y, x = rowcol( + dsm_transf, + row["geometry"].centroid.x, + row["geometry"].centroid.y, + ) # TODO: This produce different result since no standalone round coordinates + poiname.append(row[configDict["poi_field"]]) + poisxy[idx, 0] = idx + poisxy[idx, 1] = x + poisxy[idx, 2] = y + + # Create pandas datetime object to be used when createing an xarray DataSet where wall temperatures/radiation is stored and eventually saved as a NetCDf + if configDict["wallnetcdf"] == 1: + met_for_xarray = ( + pd.to_datetime(YYYY[0][:], format="%Y") + + pd.to_timedelta(DOY - 1, unit="d") + + pd.to_timedelta(hours, unit="h") + + pd.to_timedelta(minu, unit="m") + ) + else: + wallScheme = 0 + voxelMaps = 0 + voxelTable = 0 + timeStep = 0 + # thermal_effusivity = 0 + walls_scheme = np.ones((rows, cols)) * 10.0 + dirwalls_scheme = np.ones((rows, cols)) * 10.0 + + # Initialisation of time related variables + if Ta.__len__() == 1: + timestepdec = 0 + else: + timestepdec = dectime[1] - dectime[0] + timeadd = 0.0 + # timeaddE = 0. + # timeaddS = 0. + # timeaddW = 0. + # timeaddN = 0. + firstdaytime = 1.0 + + # Save hemispheric image + if anisotropic_sky == 1: + if not poisxy is None: + patch_characteristics = hemispheric_image( + poisxy, shmat, vegshmat, vbshvegshmat, voxelMaps, wallScheme + ) + + # If metfile starts at night + CI = 1.0 + + # Main loop + tmrtplot = np.zeros((rows, cols)) + + # Initiate array for I0 values + if np.unique(DOY).shape[0] > 1: + unique_days = np.unique(DOY) + first_unique_day = DOY[DOY == unique_days[0]] + I0_array = np.zeros((first_unique_day.shape[0])) + else: + first_unique_day = DOY.copy() + I0_array = np.zeros((DOY.shape[0])) + + if standAlone == 1: + progress = tqdm(total=Ta.__len__()) + + for i in np.arange(0, Ta.__len__()): + if feedback is not None: + feedback.setProgress( + int(i * (100.0 / Ta.__len__())) + ) # move progressbar forward + if feedback.isCanceled(): + feedback.setProgressText("Calculation cancelled") + break + else: + progress.update(1) + + # Daily water body temperature + if landcover == 1: + if ((dectime[i] - np.floor(dectime[i]))) == 0 or (i == 0): + Twater = np.mean(Ta[jday[0] == np.floor(dectime[i])]) + # Nocturnal cloudfraction from Offerle et al. 2003 + if (dectime[i] - np.floor(dectime[i])) == 0: + daylines = np.where(np.floor(dectime) == dectime[i]) + if daylines.__len__() > 1: + alt = altitude[0][daylines] + alt2 = np.where(alt > 1) + rise = alt2[0][0] + [_, CI, _, _, _] = clearnessindex_2013b( + zen[0, i + rise + 1], + jday[0, i + rise + 1], + Ta[i + rise + 1], + RH[i + rise + 1] / 100.0, + radG[i + rise + 1], + location, + P[i + rise + 1], + ) + if (CI > 1.0) or (CI == np.inf): + CI = 1.0 + else: + CI = 1.0 + + # Only if Kdir is derived from horizontal global shortwave and horizontal diffuse shortwave + # if altitude[0][i] > 0: + # radI[i] = radI[i]/np.sin(altitude[0][i] * np.pi/180) + # else: + # radG[i] = 0. + # radD[i] = 0. + # radI[i] = 0. + + # Timestep of the simulation used in the ground scheme calculation + first_timestep = ( + pd.to_datetime(YYYY[0][0], format="%Y") + + pd.to_timedelta(DOY[0] - 1, unit="d") + + pd.to_timedelta(hours[0], unit="h") + + pd.to_timedelta(minu[0], unit="m") + ) + second_timestep = ( + pd.to_datetime(YYYY[0][1], format="%Y") + + pd.to_timedelta(DOY[1] - 1, unit="d") + + pd.to_timedelta(hours[1], unit="h") + + pd.to_timedelta(minu[1], unit="m") + ) + timeStep = (second_timestep - first_timestep).seconds + + ( + 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, + radIout, + radDout, + Lside, + Lsky_patch_characteristics, + CI_TgG, + KsideD, + dRad, + Kside, + steradians, + voxelTable, + Rn, + Rn_past, + Tm, + G, + ) = so.Solweig_2026a_calc( + i, + dsm, + scale, + rows, + cols, + svf, + svfN, + svfW, + svfE, + svfS, + svfveg, + svfNveg, + svfEveg, + svfSveg, + svfWveg, + svfaveg, + svfEaveg, + svfSaveg, + svfWaveg, + svfNaveg, + vegdsm, + vegdsm2, + albedo_b, + absK, + absL, + ewall, + Fside, + Fup, + Fcyl, + altitude[0][i], + azimuth[0][i], + zen[0][i], + jday[0][i], + usevegdem, + onlyglobal, + buildings, + location, + psi[0][i], + landcover, + lcgrid, + dectime[i], + altmax[0][i], + wallaspect, + wallheight, + cyl, + elvis, + Ta[i], + RH[i], + radG[i], + radD[i], + radI[i], + P[i], + 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, + diffsh, + shmat, + vegshmat, + vbshvegshmat, + anisotropic_sky, + asvf, + patch_option, + voxelMaps, + voxelTable, + Ws[i], + wallScheme, + timeStep, + steradians, + walls_scheme, + dirwalls_scheme, + groundScheme, + outgoingLW, + Tg, + Rn, + Rn_past, + G, + Tm, + cap_grid, + diff_grid, + a1_grid, + a2_grid, + a3_grid, + shadow, + ) + + # Save I0 for I0 vs. Kdown output plot to check if UTC is off + if i < first_unique_day.shape[0]: + I0_array[i] = I0 + elif i == first_unique_day.shape[0]: + # Output I0 vs. Kglobal plot + radG_for_plot = radG[DOY == first_unique_day[0]] + # hours_for_plot = hours[DOY == first_unique_day[0]] + dectime_for_plot = dectime[DOY == first_unique_day[0]] + fig, ax = plt.subplots() + ax.plot(dectime_for_plot, I0_array, label="I0") + ax.plot(dectime_for_plot, radG_for_plot, label="Kglobal") + ax.set_ylabel("Shortwave radiation [$Wm^{-2}$]") + ax.set_xlabel("Decimal time") + ax.set_title("UTC" + str(configDict["utc"])) + ax.legend() + fig.savefig(configDict["output_dir"] + "/metCheck.png", dpi=150) + + tmrtplot = tmrtplot + Tmrt + + if altitude[0][i] > 0: + w = "D" + else: + w = "N" + + # Write to POIs + if not poisxy is None: + for k in range(0, poisxy.shape[0]): + poi_save = np.zeros((1, 40)) + poi_save[0, 0] = YYYY[0][i] + poi_save[0, 1] = jday[0][i] + poi_save[0, 2] = hours[i] + poi_save[0, 3] = minu[i] + poi_save[0, 4] = dectime[i] + poi_save[0, 5] = altitude[0][i] + poi_save[0, 6] = azimuth[0][i] + poi_save[0, 7] = radIout + poi_save[0, 8] = radDout + poi_save[0, 9] = radG[i] + poi_save[0, 10] = Kdown[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 11] = Kup[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 12] = Keast[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 13] = Ksouth[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 14] = Kwest[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 15] = Knorth[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 16] = Ldown[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 17] = Lup[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 18] = Least[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 19] = Lsouth[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 20] = Lwest[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 21] = Lnorth[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 22] = Ta[i] + poi_save[0, 23] = Tg[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 24] = RH[i] + poi_save[0, 25] = esky + poi_save[0, 26] = Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 27] = I0 + poi_save[0, 28] = CI + poi_save[0, 29] = shadow[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 30] = svf[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 31] = svfbuveg[ + int(poisxy[k, 2]), int(poisxy[k, 1]) + ] + poi_save[0, 32] = KsideI[int(poisxy[k, 2]), int(poisxy[k, 1])] + # Recalculating wind speed based on powerlaw + WsPET = (1.1 / sensorheight) ** 0.2 * Ws[i] + WsUTCI = (10.0 / sensorheight) ** 0.2 * Ws[i] + resultPET = p._PET( + Ta[i], + RH[i], + Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])], + WsPET, + mbody, + age, + ht, + activity, + clo, + sex, + ) + poi_save[0, 33] = resultPET + resultUTCI = utci.utci_calculator( + Ta[i], + RH[i], + Tmrt[int(poisxy[k, 2]), int(poisxy[k, 1])], + WsUTCI, + ) + poi_save[0, 34] = resultUTCI + poi_save[0, 35] = CI_TgG + poi_save[0, 36] = KsideD[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 37] = Lside[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 38] = dRad[int(poisxy[k, 2]), int(poisxy[k, 1])] + poi_save[0, 39] = Kside[int(poisxy[k, 2]), int(poisxy[k, 1])] + data_out = ( + configDict["output_dir"] + + "/POI_" + + str(poiname[k]) + + ".txt" + ) + # f_handle = file(data_out, 'a') + f_handle = open(data_out, "ab") + np.savetxt(f_handle, poi_save, fmt=numformat) + f_handle.close() + + # If wall temperature parameterization scheme is in use + if ( + configDict["wallscheme"] == 1 + ): # folderWallScheme: TODO: Fix for standalone + # Store wall data for output + if not woisxy is None: + for k in range(0, woisxy.shape[0]): + temp_wall = voxelTable.loc[ + ( + (voxelTable["ypos"] == woisxy[k, 2]) + & (voxelTable["xpos"] == woisxy[k, 1]) + ), + "wallTemperature", + ].to_numpy() + K_in = voxelTable.loc[ + ( + (voxelTable["ypos"] == woisxy[k, 2]) + & (voxelTable["xpos"] == woisxy[k, 1]) + ), + "K_in", + ].to_numpy() + L_in = voxelTable.loc[ + ( + (voxelTable["ypos"] == woisxy[k, 2]) + & (voxelTable["xpos"] == woisxy[k, 1]) + ), + "L_in", + ].to_numpy() + wallShade = voxelTable.loc[ + ( + (voxelTable["ypos"] == woisxy[k, 2]) + & (voxelTable["xpos"] == woisxy[k, 1]) + ), + "wallShade", + ].to_numpy() + temp_all = np.concatenate( + [temp_wall, K_in, L_in, wallShade] + ) + # temp_all = np.concatenate([temp_wall]) + # wall_data = np.zeros((1, 7 + temp_wall.shape[0])) + wall_data = np.zeros((1, 7 + temp_all.shape[0])) + # Part of file name (wallid), i.e. WOI_wallid.txt + data_out = ( + configDict["output_dir"] + + "/WOI_" + + str(woiname[k]) + + ".txt" + ) + if i == 0: + # Output file header + # header = 'yyyy id it imin dectime Ta SVF Ts' + header = ( + "yyyy id it imin dectime Ta SVF" + + " Ts" * temp_wall.shape[0] + + " Kin" * K_in.shape[0] + + " Lin" * L_in.shape[0] + + " shade" * wallShade.shape[0] + ) + # Part of file name (wallid), i.e. WOI_wallid.txt + # woiname = voxelTable.loc[((voxelTable['ypos'] == woisxy[k, 2]) & (voxelTable['xpos'] == woisxy[k, 1])), 'wallId'].to_numpy()[0] + woi_save = [] # + np.savetxt( + data_out, + woi_save, + delimiter=" ", + header=header, + comments="", + ) + # Fill wall_data with variables + wall_data[0, 0] = YYYY[0][i] + wall_data[0, 1] = jday[0][i] + wall_data[0, 2] = hours[i] + wall_data[0, 3] = minu[i] + wall_data[0, 4] = dectime[i] + wall_data[0, 5] = Ta[i] + wall_data[0, 6] = svf[int(woisxy[k, 2]), int(woisxy[k, 1])] + wall_data[0, 7:] = temp_all + + # Num format for output file data + woi_numformat = ( + "%d %d %d %d %.5f %.2f %.2f" + + " %.2f" * temp_all.shape[0] + ) + # Open file, add data, save + f_handle = open(data_out, "ab") + np.savetxt(f_handle, wall_data, fmt=woi_numformat) + f_handle.close() + + # Save wall temperature/radiation as NetCDF TODO: fix for standAlone? + if configDict["wallnetcdf"] == "1": # wallNetCDF: + netcdf_output = configDict["outputDir"] + "/walls.nc" + walls_as_netcdf( + voxelTable, + rows, + cols, + met_for_xarray, + i, + dsm, + configDict["filepath_dsm"], + netcdf_output, + ) + + if hours[i] < 10: + XH = "0" + else: + XH = "" + if minu[i] < 10: + XM = "0" + else: + XM = "" + + time_code = ( + str(int(YYYY[0, i])) + + "_" + + str(int(DOY[i])) + + "_" + + XH + + str(int(hours[i])) + + XM + + str(int(minu[i])) + + w + ) + + if configDict["outputtmrt"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Tmrt_" + time_code + ".tif", + Tmrt, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Tmrt_" + time_code + ".tif", + Tmrt, + dsm_transf, + dsm_crs, + ) + if configDict["outputkup"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Kup_" + time_code + ".tif", + Kup, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Kup_" + time_code + ".tif", + Kup, + dsm_transf, + dsm_crs, + ) + if configDict["outputkdown"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Kdown_" + time_code + ".tif", + Kdown, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Kdown_" + time_code + ".tif", + Kdown, + dsm_transf, + dsm_crs, + ) + if configDict["outputlup"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Lup_" + time_code + ".tif", + Lup, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Lup_" + time_code + ".tif", + Lup, + dsm_transf, + dsm_crs, + ) + if configDict["outputldown"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Ldown_" + time_code + ".tif", + Ldown, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Ldown_" + time_code + ".tif", + Ldown, + dsm_transf, + dsm_crs, + ) + if configDict["outputsh"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Shadow_" + time_code + ".tif", + shadow, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Shadow_" + time_code + ".tif", + shadow, + dsm_transf, + dsm_crs, + ) + + if configDict["outputkdiff"] == "1": + if standAlone == 0: + saveraster( + gdal_dsm, + configDict["output_dir"] + "/Kdiff_" + time_code + ".tif", + dRad, + ) + else: + common.save_raster( + configDict["output_dir"] + "/Kdiff_" + time_code + ".tif", + dRad, + dsm_transf, + dsm_crs, + ) + + # Sky view image of patches + if (anisotropic_sky == 1) & (i == 0) & (not poisxy is None): + for k in range(poisxy.shape[0]): + Lsky_patch_characteristics[:, 2] = patch_characteristics[:, k] + skyviewimage_out = ( + configDict["output_dir"] + + "/POI_" + + str(poiname[k]) + + ".png" + ) + PolarBarPlot( + Lsky_patch_characteristics, + altitude[0][i], + azimuth[0][i], + "Hemisphere partitioning", + skyviewimage_out, + 0, + 5, + 0, + ) + + # Save files for Tree Planter + if configDict["outputtreeplanter"] == "1": # outputTreeplanter: + if feedback is not None: + feedback.setProgressText("Saving files for Tree Planter tool") + # Save DSM + copyfile( + configDict["filepath_dsm"], configDict["output_dir"] + "/DSM.tif" + ) + + # Save CDSM + if usevegdem == 1: + copyfile( + configDict["filepath_cdsm"], + configDict["output_dir"] + "/CDSM.tif", + ) + + albedo_g = param["Albedo"]["Effective"]["Value"]["Cobble_stone_2014a"] + eground = param["Emissivity"]["Value"]["Cobble_stone_2014a"] + + # Saving settings from SOLWEIG for SOLWEIG1D in TreePlanter + settingsHeader = "UTC, posture, onlyglobal, landcover, anisotropic, cylinder, albedo_walls, albedo_ground, emissivity_walls, emissivity_ground, absK, absL, elevation, patch_option" + settingsFmt = ( + "%i", + "%i", + "%i", + "%i", + "%i", + "%i", + "%1.2f", + "%1.2f", + "%1.2f", + "%1.2f", + "%1.2f", + "%1.2f", + "%1.2f", + "%i", + ) + settingsData = np.array( + [ + [ + int(configDict["utc"]), + pos, + onlyglobal, + landcover, + anisotropic_sky, + cyl, + albedo_b, + albedo_g, + ewall, + eground, + absK, + absL, + alt, + patch_option, + ] + ] + ) + # print(settingsData) + np.savetxt( + configDict["output_dir"] + "/treeplantersettings.txt", + settingsData, + fmt=settingsFmt, + header=settingsHeader, + delimiter=" ", + ) + + # Copying met file for SpatialTC + copyfile( + configDict["input_met"], configDict["output_dir"] + "/metforcing.txt" + ) + + tmrtplot = ( + tmrtplot / Ta.__len__() + ) # fix average Tmrt instead of sum, 20191022 + if standAlone == 0: + saveraster( + gdal_dsm, configDict["output_dir"] + "/Tmrt_average.tif", tmrtplot + ) + else: + common.save_raster( + configDict["output_dir"] + "/Tmrt_average.tif", + tmrtplot, + dsm_transf, + dsm_crs, + ) diff --git a/functions/SOLWEIGpython/ground_surface.py b/functions/SOLWEIGpython/ground_surface.py new file mode 100644 index 0000000..c0aaf53 --- /dev/null +++ b/functions/SOLWEIGpython/ground_surface.py @@ -0,0 +1,695 @@ +""" +@author: Eliott Bridoux, University of Gothenburg +""" + +import numpy as np +import math +import matplotlib.pyplot as plt + +# Stefan-Boltzmann s constant +SBC = 5.67e-8 + + +def saturated_vp(T): + """ + Used in the calculation of the surface temperature for the water bodies. + + :Parameters: + :T: the temperature just above the water + + :Return: + :qs: the saturated vapor pressure + :slope: the slope of the tangent to the saturated vapor pressure-temperature curve + """ + + L = 2.260e6 # Latent heat of vaporisation + R = 8.314 # Gas constant + + # August-Roche-Magnus approx. in Pa + qs = 6109.4 * np.exp(17.625 * T / (T + 243.04)) + + # Clausius-Clapeyron equation + slope = L * qs / R / (T + 273.15) ** 2 + + return slope, qs + + +def initiate_groundScheme(lc_grid, solweig_parameters, day, Ta, location): + """ + Setup the maps used in the ground scheme calculations depending on the landcover + + :Parameters: + :lc_grid: Array of landcover values + :solweig_param: Dict of physical parameters + :day: day of the year (int) + :Ta: air temperature (float) + :location: Dict containing latitude and longitude + + :Return: Array for the surface temperature, radiation flux and physical parameters + """ + + # Get the landcover data from lc_grid array + lc_grid[lc_grid >= 100] = 2 + id = np.unique(lc_grid) + id = id.astype(int) + + # Physical parameters grids + cap_grid = np.copy(lc_grid) # Heat capacity + diff_grid = np.copy(lc_grid) # Thermal diffusivity + a1_grid = np.copy(lc_grid) + a2_grid = np.copy(lc_grid) + a3_grid = np.copy(lc_grid) # The 3 OHM coefficients + + # Initial fluxes and temperature grids + Rn = np.zeros_like(lc_grid) # Net radiation + Rn_past = np.zeros_like(lc_grid) # Stored net radiation + G = np.zeros_like(lc_grid) # Ground heat flux + Tg = np.copy(lc_grid) # Surface temperature + Tm = np.copy(lc_grid) # Mean daily surface temperature + + for i in id: + cap_grid[cap_grid == i] = solweig_parameters["Heat capacity"]["Value"][ + solweig_parameters["Names"]["Value"][str(i)] + ] + diff_grid[diff_grid == i] = solweig_parameters["Thermal_diffusivity"][ + "Value" + ][solweig_parameters["Names"]["Value"][str(i)]] + + # Coefficients of the OHM per land cover + mean_a1 = solweig_parameters["OHM_coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][0] + phi_a1 = solweig_parameters["OHM_coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][1] + a1_grid[a1_grid == i] = mean_a1 * ( + 1 + + 0.33 + * np.sin(2 * np.pi / 365.25 * day + phi_a1) + * np.sign(location["latitude"]) + ) + a2_grid[a2_grid == i] = solweig_parameters["OHM_coefficients"][ + "Values" + ][solweig_parameters["Names"]["Value"][str(i)]][2] + a3_grid[a3_grid == i] = solweig_parameters["OHM_coefficients"][ + "Values" + ][solweig_parameters["Names"]["Value"][str(i)]][3] + + # Initial ground surface temperature parameters + offset_Tg = solweig_parameters["Tg_ini coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][0] + slope_Tg = solweig_parameters["Tg_ini coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][1] + ratio_Tg = solweig_parameters["Tg_ini coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][2] + phi_Tg = 1.6 + + # Correct the offset value given the latitude + offset_Tg += slope_Tg * location["latitude"] + + # Mean daily soil temperature parameters + ampl_Tm = solweig_parameters["Tm_ini coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][0] + phi_Tm = 1.7 + slope_Tm = solweig_parameters["Tm_ini coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][1] + offset_Tm = solweig_parameters["Tm_ini coefficients"]["Values"][ + solweig_parameters["Names"]["Value"][str(i)] + ][2] + + # Correct the offset value given the latitude + offset_Tm += slope_Tm * location["latitude"] + + if i == 0 or i == 1: + # For paved and asphalt landcover + Tg[Tg == i] = ( + Ta[0] + + offset_Tg + * ( + 1 + + ratio_Tg + * np.sin(2 * np.pi / 365.25 * day + phi_Tg) + * np.sign(location["latitude"]) + ) + + 4 + ) + Tm[Tm == i] = ( + np.mean(Ta) + + ampl_Tm + * np.sin(2 * np.pi / 365.25 * day + phi_Tm) + * np.sign(location["latitude"]) + + offset_Tm + + 4 + ) + + elif i == 2: + # For roofs + Tg[Tg == i] = ( + Ta[0] + + offset_Tg + * ( + 1 + + ratio_Tg + * np.sin(2 * np.pi / 365.25 * day + phi_Tg) + * np.sign(location["latitude"]) + ) + + 4 + ) + Tm[Tm == i] = np.mean(Ta) + offset_Tm + + elif i == 5: + # For grass surfaces + Tg[Tg == i] = Ta[0] + offset_Tg * ( + 1 + + ratio_Tg + * np.sin(2 * np.pi / 365.25 * day + phi_Tg) + * np.sign(location["latitude"]) + ) + Tm[Tm == i] = ( + np.mean(Ta) + + ampl_Tm + * np.sin(2 * np.pi / 365.25 * day + phi_Tm) + * np.sign(location["latitude"]) + + offset_Tm + ) + + elif i == 6: + # For bare soil landcover + Tg[Tg == i] = ( + Ta[0] + + offset_Tg + * ( + 1 + + ratio_Tg + * np.sin(2 * np.pi / 365.25 * day + phi_Tg) + * np.sign(location["latitude"]) + ) + + 2 + ) + Tm[Tm == i] = ( + np.mean(Ta) + + ampl_Tm + * np.sin(2 * np.pi / 365.25 * day + phi_Tm) + * np.sign(location["latitude"]) + + offset_Tm + + 2 + ) + + elif i == 7: + # For water bodies + Tg[Tg == i] = Ta[0] + + return ( + Tg, + Tm, + Rn, + Rn_past, + G, + cap_grid, + diff_grid, + a1_grid, + a2_grid, + a3_grid, + ) + + +def surfaceTemperature_calc( + Kdown, + Ldown, + Rn, + Rn_past, + G, + Tg, + Tm, + alb, + emis, + cap, + diff, + lc_grid, + a1, + a2, + a3, + timestep, + RH, + shadow, + shadow_past, +): + """ + Calculation of the ground surface temperature + + Based on the force restore method with the ground heat flux modeled according to the heat storage flux formulation depicted in the OHM (Grimmond 1991). + A simple model is implemented to assess the surface temperature for water landcover (lc==7) since the OHM fails to model the corresponding ground heat + flux. The temporal integration is done using the Runge-Kutta 2nd order scheme + + :Parameters: + :Kdown: the downwelling shortwave radiation + :Ldown: the downwelling longwave radiation + :Rn_past: net radiation stored to calculate the radiation rate + :G: ground heat flux stored from the previous timestep + :Tg: ground surface temperature grid + :Tm: Temperature of the deep soil + :alb: albedo grid of the ground surface + :emis: emissivity grid of the surface + :cap: heat capacity grid of the material + :diff: thermal diffusivity grid of the material + :lc_grid: landcover grid to identify the water bodies + :a: coefficient grids related to the OHM + :timestep: the timestep between 2 iteration of the simulation in min + + :Return: + :Tg: Surface temperature + :Rn: Net radiation for the current timestep + :Rn_past: Stored net radiation for upcoming radiation rate calc + :G: Ground heat flux for the current timestep + """ + + # Store the past ground surface temperature + Tg_stored = Tg + + # Damping depths of the daily surface temperature wave + D = np.sqrt((2 * diff) / (2 * math.pi / 86400)) + + ### Runge Kutta method for the surface temperature calc + # First estimate of the surface temperature and of the deep soil temperature given past ground heat flux + k1 = 2 * G / cap / D - 2 * math.pi / 86400 * (Tg - Tm) + Tg_temp = Tg + k1 * timestep + + ### Estimate k2 the surface temperature step based on updated heat fluxes + # The fluxes involved are calculated using the estimated Ts + Lup_temp = SBC * emis * (Tg_temp + 273.15) ** 4 + Ldown * ( + 1 - emis + ) # Temporary outgoing longwave rad (W.m-2) + Rn_temp = Kdown * (1 - alb) + Ldown - Lup_temp # Temporary net rad (W.m-2) + RnStar_temp = (Rn_temp - Rn) / 1 # Temporary radiation rate (W.m-2.h-1) + G_temp = ( + a1 * Rn_temp + a2 * RnStar_temp + a3 + ) # Temporary ground heat flux (W.m-2) + + # Damping of the ground heat flux if it increases (or drops) too quickly + deltaG = abs(G_temp - G) + radCriterion = abs( + a1 * (Rn_temp - Rn_past) + ) # Criterion regarding the radiation step + mask = np.logical_and( + deltaG > radCriterion, abs(shadow - shadow_past) > 0.5 + ) # Grid of the pixels where the ground heat flux spikes + G_temp[mask] = G[mask] + np.sign(G_temp - G)[mask] * radCriterion[mask] + + # Correction of the temperature estimates + k2 = 2 * G_temp / cap / D - 2 * math.pi / 86400 * (Tg_temp - Tm) + Tg += (k1 + k2) / 2 * timestep + + ### Finally calculation of the updated heat fluxes + Rn_past = Rn + G_past = G + Lup = SBC * emis * (Tg + 273.15) ** 4 + (1 - emis) * Ldown + Rn = (1 - alb) * Kdown + Ldown - Lup + Rn_star = (Rn - Rn_past) / 1 + G = a1 * Rn + a2 * Rn_star + a3 + + # Damping of the ground heat flux if it increases (or decreases) too quickly + deltaG = abs(G - G_past) + radCriterion = abs( + a1 * (Rn - Rn_past) + ) # Criterion regarding the radiation step + mask = np.logical_and( + deltaG > radCriterion, abs(shadow - shadow_past) > 0.5 + ) # Grid of the pixels where the ground heat flux spikes + G[mask] = G_past[mask] + np.sign(G - G_past)[mask] * radCriterion[mask] + + ### Water bodies surface temperature estimate + beta = 0.45 # Amount of shortwave rad absorbed by the water surface layer + thickness = 1 # Depth of the water layer + lamb = 2.260e6 # Latent heat of vaporisation + rho = 1000 # Density of water (kg.m-3) + Rn_water = ( + Kdown * (1 - alb) * (beta + (1 - beta) * (1 - np.exp(-1))) + + Ldown + - Lup + ) # Net radiation for the top water layer beta described the transmitted rad + _, es = saturated_vp(Tg) + E = 0.0858 * (es / 1000) * (1 - RH / 100) / 3600 / 1000 * rho * lamb + deltaTg = np.copy(lc_grid) + deltaTg = ( + timestep + / cap + / thickness + * (Rn_water - E - diff * cap / thickness * (Tg - Tm)) + ) + Tg[lc_grid == 7] = Tg_stored[lc_grid == 7] + deltaTg[lc_grid == 7] + + return Tg, Rn, Rn_past, G + + +def outgoingLongwave_calc( + Tg, + Tgwall, + Ta, + Ldown, + emis, + alb, + buildings, + shadow, + sunwall, + walls, + rows, + cols, + sizepx, +): + """ + Calculation of the outgoing longwave radiation from the ground, + + :Parameters: + :Tg: ground surface temperature grid + :Tgwall: wall surface temperature grid + :Ta: air temperature grid + :emis: emissivity grid of the surface + :alb: albedo grid of the surface + :emis_wall: emissivity of the wall (for now float = 0.9) + :buildings: boolean grid, 0 if the landcover is roof and 1 if there is no building + :shadow: boolean grid, 0 when the pixel is shadowed and 1 when sunlit + :sunwall: Grid where non zero values indicate a sunlit wall and its height + :walls: grid containing the heights of the walls + :aspect: grid containing the angles of the normal dir to walls + :rows: number of rows in the grids + :cols: number of columns in the grids + :sizepx: size of a pixel in m (1/scale) + + :Return: + """ + + # Assessment of the distance from a pixel at which most of the radiation are received (cf view factor Lambert) + factor = 0.99 # Percentage of radiation accounted for + zs = 1.1 # in m + r_max = zs * np.sqrt( + factor / (1 - factor) + ) # in m, maximum radius for the radiation calc + + # Emissivity of the wall + emis_wall = 0.9 + alb_wall = 0.2 + + # Copy of the sunlit wall grid and replacement of the wall height with 1 if sunlit + sunlitwall = sunwall + sunlitwall[sunlitwall > 0] = 1 + + # The alb grids only take into account the sunlit surfaces in the alb calculation albnosh calculate it for all the surfaces + albsunlit = alb * shadow + + # Boolean array 1 if the pixel is a wall, 0 if not + wallbol = (walls > 0) * 1 + + # step in meters between every iteration + step = 1 + + ### Initialize the ground view factor grids as np.zeros() + # Upwelling longwave radiation + gvfLup = np.zeros((rows, cols)) + gvfLupE = np.zeros((rows, cols)) + gvfLupS = np.zeros((rows, cols)) + gvfLupW = np.zeros((rows, cols)) + gvfLupN = np.zeros((rows, cols)) + + # Albedo of the sunlit surfaces + gvfalbsun = np.zeros((rows, cols)) + gvfalbsunE = np.zeros((rows, cols)) + gvfalbsunS = np.zeros((rows, cols)) + gvfalbsunW = np.zeros((rows, cols)) + gvfalbsunN = np.zeros((rows, cols)) + + # Albedo complete + gvfalbtot = np.zeros((rows, cols)) + gvfalbtotE = np.zeros((rows, cols)) + gvfalbtotS = np.zeros((rows, cols)) + gvfalbtotW = np.zeros((rows, cols)) + gvfalbtotN = np.zeros((rows, cols)) + + # Longwave radiation coming from the side + gvfLsideE = np.zeros((rows, cols)) + gvfLsideS = np.zeros((rows, cols)) + gvfLsideW = np.zeros((rows, cols)) + gvfLsideN = np.zeros((rows, cols)) + + # Add the radiation from the pixel directly below, only for the total gvf + # Do not take the roofs into account for now + view_factor = (sizepx / 2) ** 2 / ((sizepx / 2) ** 2 + zs**2) + gvfLup = ( + gvfLup + (SBC * emis * (Tg + 273.15) ** 4) * view_factor * buildings + ) + gvfalbsun = gvfalbsun + albsunlit * view_factor * buildings + gvfalbtot = gvfalbtot + alb * view_factor * buildings + + # Division of the 360° field of view in 20 and convert the array in radian + azimuths = np.linspace(18, 360, num=20, endpoint=True) + azimuths = azimuths * (np.pi / 180) + + ### Loop for the number of azimuth values + for azimuth in azimuths: + # Copy of the building grid + building_copy = buildings + + # Initialisation of the tables + # First the ones containing the translated rasters (temporary) + building_temp = np.zeros((rows, cols)) + Lup_temp = np.zeros((rows, cols)) + Lwall_temp = np.zeros((rows, cols)) + albsun_temp = np.zeros((rows, cols)) + albtot_temp = np.zeros((rows, cols)) + sunlitwall_temp = np.zeros((rows, cols)) + + # Then the tables containing the sum of the radiations (or albedo) for this azimuth + Lup_sum = np.zeros((rows, cols)) + LsideE_sum = np.zeros((rows, cols)) + LsideN_sum = np.zeros((rows, cols)) + LsideW_sum = np.zeros((rows, cols)) + LsideS_sum = np.zeros((rows, cols)) + albsun_sum = np.zeros((rows, cols)) + albtot_sum = np.zeros((rows, cols)) + + ### Shadow casting algorithm + # Translation ranges from 1/2 a pixel to the max radius r_max + for r in np.arange(sizepx / 2, r_max, step=step): + # Longwave radiation grids both at the ground level and from the walls + Lup = ( + SBC * emis * (Tg + 273.15) ** 4 + Ldown * (1 - emis) + ) * building_copy + Lwall = ( + SBC * emis_wall * (Tgwall + Ta + 273.15) ** 4 * building_copy + ) + + # Step of the raster translation + dx = -np.cos(azimuth) + dy = -np.sin(azimuth) + + # Scale so that the grid is at least translated from 1px + if abs(dx) > abs(dy): + dx = -r * np.sign(np.cos(azimuth)) + dy = -r * abs(np.tan(azimuth)) * np.sign(np.sin(azimuth)) + else: + dx = -r / abs(np.tan(azimuth)) * np.sign(np.cos(azimuth)) + dy = -r * np.sign(np.sin(azimuth)) + + # Select the interested part of the initial raster and the translated one from their four corners and + # translating toward the direction azimuth = 0° for dx > 0 + if dx > 0: + x_select_start = dx + x_select_end = rows + x_transl_start = 0 + x_transl_end = rows - dx + else: + x_select_start = 0 + x_select_end = rows + dx + x_transl_start = -dx + x_transl_end = rows + + # translating toward the direction azimuth = 90° for dy > 0 + if dy > 0: + y_select_start = dy + y_select_end = cols + y_transl_start = 0 + y_transl_end = cols - dy + else: + y_select_start = 0 + y_select_end = cols + dy + y_transl_start = -dy + y_transl_end = cols + + # Copy the initial rasters and input inside translated raster temporary and changing every iteration + # Building grid + building_temp[ + int(x_transl_start) : math.ceil(x_transl_end), + int(y_transl_start) : math.ceil(y_transl_end), + ] = buildings[ + int(x_select_start) : math.ceil(x_select_end), + int(y_select_start) : math.ceil(y_select_end), + ] + + # Ground longwave radiation grid + Lup_temp[ + int(x_transl_start) : math.ceil(x_transl_end), + int(y_transl_start) : math.ceil(y_transl_end), + ] = Lup[ + int(x_select_start) : math.ceil(x_select_end), + int(y_select_start) : math.ceil(y_select_end), + ] + + # Wall longwave radiation grid + Lwall_temp[ + int(x_transl_start) : math.ceil(x_transl_end), + int(y_transl_start) : math.ceil(y_transl_end), + ] = Lwall[ + int(x_select_start) : math.ceil(x_select_end), + int(y_select_start) : math.ceil(y_select_end), + ] + + # Albedo grid for the sunlit area + albsun_temp[ + int(x_transl_start) : math.ceil(x_transl_end), + int(y_transl_start) : math.ceil(y_transl_end), + ] = albsunlit[ + int(x_select_start) : math.ceil(x_select_end), + int(y_select_start) : math.ceil(y_select_end), + ] + + # Albedo grid for all the area + albtot_temp[ + int(x_transl_start) : math.ceil(x_transl_end), + int(y_transl_start) : math.ceil(y_transl_end), + ] = alb[ + int(x_select_start) : math.ceil(x_select_end), + int(y_select_start) : math.ceil(y_select_end), + ] + + # Sunlit wall grid + sunlitwall_temp[ + int(x_transl_start) : math.ceil(x_transl_end), + int(y_transl_start) : math.ceil(y_transl_end), + ] = sunlitwall[ + int(x_select_start) : math.ceil(x_select_end), + int(y_select_start) : math.ceil(y_select_end), + ] + + # Change the boolean building grid, if the px was already a building it remains one (px value = 0) + building_copy = np.min([building_copy, building_temp], axis=0) + + # For each pixel add the translated Lup to the received rad if there where there is no building + view_factor = ((r + step) ** 2 / (zs**2 + (r + step) ** 2)) - ( + r**2 / (zs**2 + r**2) + ) + Lup_sum += Lup_temp * view_factor * building_copy / 20 + albsun_sum += albsun_temp * view_factor * building_copy / 20 + albtot_sum += albtot_temp * view_factor * building_copy / 20 + + # Create a boolean grid to assert that the sunlit walls are not inside a building + wall_temp = wallbol * building_copy + onlysunwall_temp = sunlitwall_temp * building_copy + + # Then add the radiation incoming from those walls + Lup_sum += ( + wall_temp * Lwall_temp * zs**2 / ((r + step) ** 2 + zs**2) / 20 + ) + albsun_sum += ( + onlysunwall_temp + * alb_wall + * zs**2 + / ((r + step) ** 2 + zs**2) + / 20 + ) + albtot_sum += ( + wall_temp * alb_wall * zs**2 / ((r + step) ** 2 + zs**2) / 20 + ) + + # Finally add the radiation received from the side + dphi = np.arctan((r + step) / zs) - np.arctan(r / zs) + dtrigo = zs / np.sqrt(r**2 + zs**2) * r / np.sqrt( + r**2 + zs**2 + ) - zs / np.sqrt((r + step) ** 2 + zs**2) * (r + step) / np.sqrt( + (r + step) ** 2 + zs**2 + ) + + # Calculation of the solid angle for each of the cardinal points + steradiansW, steradiansS, steradiansE, steradiansN = 0, 0, 0, 0 + if (azimuth >= 0) and (azimuth < np.pi): + dthetaW = 2 * np.pi / 20 + steradiansW += dthetaW * (dphi + dtrigo) / 2 + + if (azimuth >= np.pi / 2) and (azimuth < 3 * np.pi / 2): + dthetaS = 2 * np.pi / 20 + steradiansS += dthetaS * (dphi + dtrigo) / 2 + + if (azimuth >= np.pi) and (azimuth < 2 * np.pi): + dthetaE = 2 * np.pi / 20 + steradiansE += dthetaE * (dphi + dtrigo) / 2 + + if (azimuth >= 3 * np.pi / 2) or (azimuth < np.pi / 2): + dthetaN = 2 * np.pi / 20 + steradiansN += dthetaN * (dphi + dtrigo) / 2 + + LsideW_sum += Lup_temp / np.pi * steradiansW * building_copy + LsideS_sum += Lup_temp / np.pi * steradiansS * building_copy + LsideE_sum += Lup_temp / np.pi * steradiansE * building_copy + LsideN_sum += Lup_temp / np.pi * steradiansN * building_copy + + ### Add the value for the computed part of the field of view + gvfLup += Lup_sum + gvfalbsun += albsun_sum + gvfalbtot += albtot_sum + + # Add the value if the azimuth correspond to the side of the compass + if (azimuth >= 0) and (azimuth < np.pi): + gvfLupW += Lup_sum + gvfalbsunW += albsun_sum + gvfalbtotW += albtot_sum + gvfLsideW += LsideW_sum + + if (azimuth >= np.pi / 2) and (azimuth < 3 * np.pi / 2): + gvfLupS += Lup_sum + gvfalbsunS += albsun_sum + gvfalbtotS += albtot_sum + gvfLsideS += LsideS_sum + + if (azimuth >= np.pi) and (azimuth < 2 * np.pi): + gvfLupE += Lup_sum + gvfalbsunE += albsun_sum + gvfalbtotE += albtot_sum + gvfLsideE += LsideE_sum + + if (azimuth >= 3 * np.pi / 2) or (azimuth < np.pi / 2): + gvfLupN += Lup_sum + gvfalbsunN += albsun_sum + gvfalbtotN += albtot_sum + gvfLsideN += LsideN_sum + + # If the px is associated with a roof landcover, for now Lup = 0 + # Here their Lup value is allocated to those px + gvfLup += (SBC * emis * (Tg + 273.15) ** 4) * (buildings * -1 + 1) + + # # Finally add the reflection from the downwelling longwave radiation + # gvfLup += Ldown * (1-emis) + + return ( + gvfLup, + gvfalbsun, + gvfalbtot, + gvfLupE, + gvfalbsunE, + gvfalbtotE, + gvfLupS, + gvfalbsunS, + gvfalbtotS, + gvfLupW, + gvfalbsunW, + gvfalbtotW, + gvfLupN, + gvfalbsunN, + gvfalbtotN, + gvfLsideW, + gvfLsideS, + gvfLsideE, + gvfLsideN, + ) diff --git a/parametersforsolweig.json b/parametersforsolweig.json new file mode 100644 index 0000000..167012f --- /dev/null +++ b/parametersforsolweig.json @@ -0,0 +1,362 @@ +{ + "Names": { + "Value": { + "0": "Cobble_stone_2014a", + "1": "Dark_asphalt", + "2": "Roofs(buildings)", + "5": "Grass_unmanaged", + "6": "Bare_soil", + "7": "Water", + "99": "Walls", + "100": "Brick_wall", + "101": "Concrete_wall", + "102": "Wood_wall" + }, + "Comment": "Name of each respective land cover class in land cover data." + }, + "Code": { + "Value": { + "Cobble_stone_2014a": 0, + "Dark_asphalt": 1, + "Roofs(buildings)": 2, + "Grass_unmanaged": 5, + "Bare_soil": 6, + "Water": 7, + "Walls": 99, + "Brick_wall": 100, + "Concrete_wall": 101, + "Wood_wall": 102 + }, + "Comment": "Code for each land cover class name." + }, + "Albedo": { + "Effective": { + "Value": { + "Cobble_stone_2014a": 0.2, + "Dark_asphalt": 0.18, + "Roofs(buildings)": 0.18, + "Grass_unmanaged": 0.16, + "Bare_soil": 0.25, + "Water": 0.05, + "Walls": 0.2 + }, + "Comment": "Effective albedos according to Lindberg et al., 2008; 2016." + }, + "Material": { + "Value": { + "Brick_wall": 0.2, + "Concrete_wall": 0.2, + "Wood_wall": 0.2 + }, + "Comment": "Material albedos according to Wallenberg et al., 2025." + } + }, + "Emissivity": { + "Value": { + "Cobble_stone_2014a": 0.95, + "Dark_asphalt": 0.95, + "Roofs(buildings)": 0.95, + "Grass_unmanaged": 0.94, + "Bare_soil": 0.94, + "Water": 0.98, + "Walls": 0.9, + "Brick_wall": 0.9, + "Concrete_wall": 0.9, + "Wood_wall": 0.9 + }, + "Comment": "Emissivity of each land cover class." + }, + "Specific_heat": { + "Value": { + "Cobble_stone_2014a": -9999.0, + "Dark_asphalt": -9999.0, + "Roofs(buildings)": -9999.0, + "Grass_unmanaged": -9999.0, + "Bare_soil": -9999.0, + "Water": -9999.0, + "Walls": -9999.0, + "Brick_wall": 800, + "Concrete_wall": 840, + "Wood_wall": 1880 + }, + "Comment": "Specific heat capacity, in units J kg-1 K-1, used for wall surface temperatures according to Wallenberg et al. 2025." + }, + "Heat capacity": { + "Value": { + "Cobble_stone_2014a": 2110000.0, + "Dark_asphalt": 1940000.0, + "Roofs(buildings)": 1250000.0, + "Grass_unmanaged": 1350000.0, + "Bare_soil": 2000000.0, + "Water": 4180000.0, + "Walls": 1250000.0, + "Brick_wall": -9999.0, + "Concrete_wall": -9999.0, + "Wood_wall": -9999.0 + }, + "Comment": "Heat capacity of each land cover class, in J m-3 K-1, equal to density times specific heat and used in the ground surface temperature scheme" + }, + "Thermal_conductivity": { + "Value": { + "Cobble_stone_2014a": 2, + "Dark_asphalt": 0.83, + "Roofs(buildings)": 1, + "Grass_unmanaged": 1.5, + "Bare_soil": 1.5, + "Water": 1, + "Walls": 1, + "Brick_wall": 0.84, + "Concrete_wall": 1.7, + "Wood_wall": 0.17 + }, + "Comment": "Thermal conductivity of each land cover class, in units W m-1 K-1, used for wall surface temperatures according to Wallenberg et al. 2025." + }, + "Thermal_diffusivity": { + "Value": { + "Cobble_stone_2014a": 0.00000072, + "Dark_asphalt": 0.00000038, + "Roofs(buildings)": 0.00000005, + "Grass_unmanaged": 0.00000021, + "Bare_soil": 0.00000030, + "Water": 0.0000001, + "Walls": 0.00000005, + "Brick_wall": -9999.0, + "Concrete_wall": -9999.0, + "Wood_wall": -9999.0 + }, + "Comment": "Thermal diffusivity of each land cover class, in , equal to thermal conductivity per density per specific heat and used in the ground surface temperature scheme" + }, + "Density": { + "Value": { + "Cobble_stone_2014a": -9999.0, + "Dark_asphalt": -9999.0, + "Roofs(buildings)": -9999.0, + "Grass_unmanaged": -9999.0, + "Bare_soil": -9999.0, + "Water": -9999.0, + "Walls": -9999.0, + "Brick_wall": 1700, + "Concrete_wall": 2200, + "Wood_wall": 700 + }, + "Comment": "Density of the material in units kg m-3, used for wall surface temperatures according to Wallenberg et al. 2025." + }, + "Wall_thickness": { + "Value": { + "Brick_wall": 0.1, + "Concrete_wall": 0.2, + "Wood_wall": 0.03 + }, + "Comment": "Wall thickness in units meters, used to calculate characteristic time for wall surface temperatures (Wallenberg et al., 2025)." + }, + "TmaxLST": { + "Value": { + "Cobble_stone_2014a": 15.0, + "Dark_asphalt": 15.0, + "Roofs(buildings)": 15.0, + "Grass_unmanaged": 14.0, + "Bare_soil": 14.0, + "Water": 12.0, + "Walls": 15.0, + "Brick_wall": -999.0, + "Concrete_wall": -999.0, + "Wood_wall": -999.0 + }, + "Comment": "TmaxLST used for ground surface temperatures and wall surface temperatures according to Lindberg et al. 2008; 2016." + }, + "Ts_deg": { + "Value": { + "Cobble_stone_2014a": 0.37, + "Dark_asphalt": 0.58, + "Roofs(buildings)": 0.58, + "Grass_unmanaged": 0.21, + "Bare_soil": 0.33, + "Water": 0.0, + "Walls": 0.37, + "Brick_wall": -999.0, + "Concrete_wall": -999.0, + "Wood_wall": -999.0 + }, + "Comment": "Ts_deg used for ground surface temperatures and wall surface temperatures according to Lindberg et al. 2008; 2016." + }, + "Tstart": { + "Value": { + "Cobble_stone_2014a": -3.41, + "Dark_asphalt": -9.78, + "Roofs(buildings)": -9.78, + "Grass_unmanaged": -3.38, + "Bare_soil": -3.01, + "Water": 0.0, + "Walls": -3.41, + "Brick_wall": -999.0, + "Concrete_wall": -999.0, + "Wood_wall": -999.0 + }, + "Comment": "Tstart used for ground surface temperatures and wall surface temperatures according to Lindberg et al. 2008; 2016." + }, + "OHM_coefficients": { + "Values": { + "Cobble_stone_2014a": [ + 0.61, + 1.39, + 0.28, + -23.9 + ], + "Dark_asphalt": [ + 0.50, + 1.39, + 0.28, + -31.45 + ], + "Roofs(buildings)": [ + 0.12, + 1.39, + 0.24, + -4.5 + ], + "Grass_unmanaged": [ + 0.27, + 1.39, + 0.33, + -21.75 + ], + "Bare_soil": [ + 0.3, + 1.39, + 0.44, + -24 + ], + "Water": [ + 0.1, + 1.39, + 0, + -10 + ] + }, + "Comment": "Used in the Objective Hysteresis Model leading to the ground surface temperature the mean of coef a1 (-), its phase (rad), a2 (h-1) and a3 (W m-2) " + }, + "Tg_ini coefficients": { + "Values": { + "Cobble_stone_2014a": [ + -1.7, + 1.5, + 1.61 + ], + "Dark_asphalt": [ + -2.35, + 0.9, + 1.60 + ], + "Roofs(buildings)": [ + 1, + -1.3, + 1.55 + ], + "Grass_unmanaged": [ + -1.4, + 1.1, + 1.58 + ], + "Bare_soil": [ + -1.8, + 0.75, + 1.58 + ], + "Water": [ + 0, + 0, + 0 + ] + }, + "Comment": "Offset in K, Ratio between amplitude and offset and phase used to model the seasonal pattern of the initial ground surface temperature when compared to the air temperature" + }, + "Tm_ini coefficients": { + "Values": { + "Cobble_stone_2014a": [ + -6.0, + 1.71, + 1.43 + ], + "Dark_asphalt": [ + -4.9, + 1.71, + -0.22 + ], + "Roofs(buildings)": [ + 0, + 0, + 0.3 + ], + "Grass_unmanaged": [ + -2.7, + 1.7, + -0.68 + ], + "Bare_soil": [ + -2.9, + 1.7, + -0.73 + ], + "Water": [ + 0, + 0, + 0 + ] + }, + "Comment": "Amplitude in K, phase and offset in K used to model the seasonal pattern of the initial ground surface temperature when compared to the air temperature" + }, + "Tmrt_params": { + "Value": { + "absK": 0.70, + "absL": 0.95, + "posture": "Standing" + }, + "Comment": "Absorption coefficients for mean radiant temperature (Tmrt) and posture. Posture is either standing or sitting." + }, + "PET_settings": { + "Value": { + "Age": 35, + "Weight": 75.0, + "Height": 180, + "Sex": "Male", + "Activity": 80.0, + "clo": 0.90 + }, + "Comment": "Settings to calculate Physiological Equivalent Temperature (PET). Sex is either Male or Female." + }, + "Wind_Height": { + "Value": { + "magl": 10.0 + }, + "Comment": "Height of wind sensor for PET and UTCI calculations." + }, + "Tree_settings":{ + "Value":{ + "Transmissivity": 0.03, + "Trunk_ratio": 0.25, + "First_day_leaf": 97, + "Last_day_leaf": 300 + }, + "Comment": "Settings for trees. Shortwave trasmissivity in %. Trunkratio as a fraction of total height" + }, + "Posture": { + "Standing" : { + "Value": { + "Fside" : 0.22, + "Fup" : 0.06, + "height" : 1.1, + "Fcyl" : 0.28 + }, + "Comment": "Standing posture of human body. Used in Tmrt calulations" + }, + "Sitting" : { + "Value":{ + "Fside" : 0.166666, + "Fup" : 0.166666, + "height" : 0.75, + "Fcyl" : 0.2 + }, + "Comment": "Sitting posture of human body. Used in Tmrt calulations" + } + } +} \ No newline at end of file diff --git a/processor/configsolweig.ini b/processor/configsolweig.ini index 45053bf..347305f 100644 --- a/processor/configsolweig.ini +++ b/processor/configsolweig.ini @@ -34,6 +34,8 @@ poi_file= poi_field= # Input file for wall temperture scheme (Wallenberg et al. 2025) input_wall = C:\Users\xlinfr\Desktop\SOLWEIGdata\blabla.npz +# Input file for surface temperature data +input_surf=C:\Users\xlinfr\Desktop\SOLWEIGdata\surftemp.txt #Point of Interest file for walls woi_file= woi_field= @@ -66,6 +68,10 @@ landcover=1 demforbuild=0 # use anisotrphic sky (Wallenberg et al. XXXX and Wallenberg et al. XXXX) aniso=0 +# use OHM for ground surface temperature modeling +groundmodel=1 +# compute the outgoing longwave rad using solid angles +outgoingLW=1 # use wall surface temperature scheme (Wallenberg et al. 2025, GMD) wallscheme=0 # If building materials is not included in lc, then this is used for all buildings (Wood, Brick or Concrete) diff --git a/processor/parametersforsolweig.json b/processor/parametersforsolweig.json index 17a8cbe..270e4c3 100644 --- a/processor/parametersforsolweig.json +++ b/processor/parametersforsolweig.json @@ -13,7 +13,7 @@ "102": "Wood_wall" }, "Comment": "Name of each respective land cover class in land cover data." - }, + }, "Code": { "Value": { "Cobble_stone_2014a": 0, @@ -28,7 +28,7 @@ "Wood_wall": 102 }, "Comment": "Code for each land cover class name." - }, + }, "Albedo": { "Effective": { "Value": { @@ -65,7 +65,7 @@ "Wood_wall": 0.9 }, "Comment": "Emissivity of each land cover class." - }, + }, "Specific_heat": { "Value": { "Cobble_stone_2014a": -9999.0, @@ -80,22 +80,52 @@ "Wood_wall": 1880 }, "Comment": "Specific heat capacity, in units J kg-1 K-1, used for wall surface temperatures according to Wallenberg et al. 2025." + }, + "Heat capacity": { + "Value": { + "Cobble_stone_2014a": 2110000.0, + "Dark_asphalt": 1940000.0, + "Roofs(buildings)": 1250000.0, + "Grass_unmanaged": 1350000.0, + "Bare_soil": 2000000.0, + "Water": 4180000.0, + "Walls": 1250000.0, + "Brick_wall": -9999.0, + "Concrete_wall": -9999.0, + "Wood_wall": -9999.0 }, + "Comment": "Heat capacity of each land cover class, in J m-3 K-1, equal to density times specific heat and used in the ground surface temperature scheme" + }, "Thermal_conductivity": { "Value": { - "Cobble_stone_2014a": -9999.0, - "Dark_asphalt": -9999.0, - "Roofs(buildings)": -9999.0, - "Grass_unmanaged": -9999.0, - "Bare_soil": -9999.0, - "Water": -9999.0, - "Walls": -9999.0, + "Cobble_stone_2014a": 2, + "Dark_asphalt": 0.83, + "Roofs(buildings)": 1, + "Grass_unmanaged": 1.5, + "Bare_soil": 1.5, + "Water": 1, + "Walls": 1, "Brick_wall": 0.84, "Concrete_wall": 1.7, "Wood_wall": 0.17 }, "Comment": "Thermal conductivity of each land cover class, in units W m-1 K-1, used for wall surface temperatures according to Wallenberg et al. 2025." + }, + "Thermal_diffusivity": { + "Value": { + "Cobble_stone_2014a": 0.00000072, + "Dark_asphalt": 0.00000038, + "Roofs(buildings)": 0.00000005, + "Grass_unmanaged": 0.00000021, + "Bare_soil": 0.00000030, + "Water": 0.0000001, + "Walls": 0.00000005, + "Brick_wall": -9999.0, + "Concrete_wall": -9999.0, + "Wood_wall": -9999.0 }, + "Comment": "Thermal diffusivity of each land cover class, in , equal to thermal conductivity per density per specific heat and used in the ground surface temperature scheme" + }, "Density": { "Value": { "Cobble_stone_2014a": -9999.0, @@ -133,7 +163,7 @@ "Wood_wall": -999.0 }, "Comment": "TmaxLST used for ground surface temperatures and wall surface temperatures according to Lindberg et al. 2008; 2016." - }, + }, "Ts_deg": { "Value": { "Cobble_stone_2014a": 0.37, @@ -148,7 +178,7 @@ "Wood_wall": -999.0 }, "Comment": "Ts_deg used for ground surface temperatures and wall surface temperatures according to Lindberg et al. 2008; 2016." - }, + }, "Tstart": { "Value": { "Cobble_stone_2014a": -3.41, @@ -163,7 +193,118 @@ "Wood_wall": -999.0 }, "Comment": "Tstart used for ground surface temperatures and wall surface temperatures according to Lindberg et al. 2008; 2016." + }, + "OHM_coefficients": { + "Values": { + "Cobble_stone_2014a": [ + 0.61, + 1.39, + 0.28, + -23.9 + ], + "Dark_asphalt": [ + 0.50, + 1.39, + 0.28, + -31.45 + ], + "Roofs(buildings)": [ + 0.12, + 1.39, + 0.24, + -4.5 + ], + "Grass_unmanaged": [ + 0.27, + 1.39, + 0.33, + -21.75 + ], + "Bare_soil": [ + 0.3, + 1.39, + 0.44, + -24 + ], + "Water": [ + 0.1, + 1.39, + 0, + -10 + ] + }, + "Comment": "Used in the Objective Hysteresis Model leading to the ground surface temperature the mean of coef a1 (-), its phase (rad), a2 (h-1) and a3 (W m-2) " + }, + "Tg_ini coefficients": { + "Values": { + "Cobble_stone_2014a": [ + -1.0, + -0.04, + 1.5 + ], + "Dark_asphalt": [ + -1.4, + -0.04, + 0.72 + ], + "Roofs(buildings)": [ + -2.9, + 0, + 0.43 + ], + "Grass_unmanaged": [ + -4.0, + -0.02, + 0.38 + ], + "Bare_soil": [ + -1.6, + -0.02, + 0.51 + ], + "Water": [ + 0, + 0, + 0 + ] }, + "Comment": "Offset in K, Slope of offset vs latitude, Ratio between amplitude and offset and phase used to model the seasonal pattern of the initial ground surface temperature when compared to the air temperature" + }, + "Tm_ini coefficients": { + "Values": { + "Cobble_stone_2014a": [ + -4.5, + -0.15, + 9 + ], + "Dark_asphalt": [ + -5.7, + -0.18, + 10 + ], + "Roofs(buildings)": [ + 0, + 0, + 6 + ], + "Grass_unmanaged": [ + -6.2, + -0.2, + 10 + ], + "Bare_soil": [ + -3.4, + -0.11, + 5 + ], + "Water": [ + 0, + 0, + 0 + ] + }, + "Comment": "Amplitude in K, slope of offset vs latitude in K/deg and offset in K used to model the seasonal pattern of the initial ground surface temperature when compared to the air temperature" + }, "Tmrt_params": { "Value": { "absK": 0.70, @@ -187,7 +328,7 @@ "Value": { "magl": 10.0 }, - "Comment": "Height of wind sensor for PET and UTCI calcualtions." + "Comment": "Height of wind sensor for PET and UTCI calculations." }, "Tree_settings":{ "Value":{ @@ -218,4 +359,4 @@ "Comment": "Sitting posture of human body. Used in Tmrt calulations" } } -} \ No newline at end of file +} diff --git a/processor/solweig_algorithm.py b/processor/solweig_algorithm.py index a8e3e6c..0ea103f 100644 --- a/processor/solweig_algorithm.py +++ b/processor/solweig_algorithm.py @@ -34,6 +34,7 @@ from qgis.core import ( QgsProcessing, QgsProcessingAlgorithm, + QgsProcessingParameterString, QgsProcessingParameterBoolean, QgsProcessingParameterNumber, QgsProcessingParameterFolderDestination, @@ -46,21 +47,52 @@ QgsProcessingParameterRasterLayer, ) +# from processing.gui.wrappers import WidgetWrapper +# from qgis.PyQt.QtWidgets import QDateEdit, QTimeEdit import numpy as np + +# import pandas as pd from osgeo import gdal +from osgeo.gdalconst import * import os from qgis.PyQt.QtGui import QIcon import inspect from pathlib import Path from ..util.misc import xy2latlon_fromraster import zipfile + +# from ..util.SEBESOLWEIGCommonFiles.Solweig_v2015_metdata_noload import Solweig_2015a_metdata_noload +# from ..util.SEBESOLWEIGCommonFiles import Solweig_v2015_metdata_noload as metload from ..util.umep_solweig_export_component import ( read_solweig_config, write_solweig_config, ) + +# from ..util.SEBESOLWEIGCommonFiles.clearnessindex_2013b import clearnessindex_2013b +# from ..functions.SOLWEIGpython.Tgmaps_v1 import Tgmaps_v1 +# from ..functions.SOLWEIGpython import Solweig_2022a_calc_forprocessing as so +# from ..functions.SOLWEIGpython import WriteMetadataSOLWEIG +# from ..functions.SOLWEIGpython import PET_calculations as p +# from ..functions.SOLWEIGpython import UTCI_calculations as utci +# from ..functions.SOLWEIGpython.CirclePlotBar import PolarBarPlot +# from ..functions.SOLWEIGpython.wall_surface_temperature import load_walls + +# from ..functions.SOLWEIGpython.wallOfInterest import wallOfInterest +# from ..functions.SOLWEIGpython.wallsAsNetCDF import walls_as_netcdf + from ..functions.SOLWEIGpython import Solweig_run as sr + +# import matplotlib.pyplot as plt import json +# For "Save necessary rasters for TreePlanter tool" +# from shutil import copyfile + +# import time + +# +# import datetime + class ProcessingSOLWEIGAlgorithm(QgsProcessingAlgorithm): """ @@ -88,6 +120,8 @@ class ProcessingSOLWEIGAlgorithm(QgsProcessingAlgorithm): INPUT_DEM = "INPUT_DEM" SAVE_BUILD = "SAVE_BUILD" INPUT_ANISO = "INPUT_ANISO" + USE_GROUNDSCHEME = "USE_GROUNDSCHEME" + INPUT_GROUNDSCHEME = "INPUT_GROUNDSCHEME" INPUT_WALLSCHEME = "INPUT_WALLSCHEME" WALLTEMP_NETCDF = "WALLTEMP_NETCDF" @@ -180,7 +214,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.TRANS_VEG, self.tr("Transmissivity of light through vegetation (%):"), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(3), True, minValue=0, @@ -194,7 +228,7 @@ def initAlgorithm(self, config): self.tr( "First day of year with leaves on trees (if deciduous)" ), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(97), False, minValue=0, @@ -208,7 +242,7 @@ def initAlgorithm(self, config): self.tr( "Last day of year with leaves on trees (if deciduous)" ), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(300), False, minValue=0, @@ -238,7 +272,7 @@ def initAlgorithm(self, config): self.tr( "Trunk zone height (percent of Canopy Height). Used if no Vegetation Trunk-zone DSM is loaded" ), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(25.0), optional=True, minValue=0.1, @@ -287,6 +321,22 @@ def initAlgorithm(self, config): optional=True, ) ) + self.addParameter( + QgsProcessingParameterBoolean( + self.USE_GROUNDSCHEME, + self.tr("Use ground cover scheme v2026a"), + defaultValue=False, + optional=True, + ) + ) + self.addParameter( + QgsProcessingParameterFile( + self.INPUT_GROUNDSCHEME, + self.tr("Ground surface temperature data (.txt)"), + extension="txt", + optional=True, + ) + ) self.addParameter( QgsProcessingParameterFile( self.INPUT_WALLSCHEME, @@ -307,6 +357,9 @@ def initAlgorithm(self, config): ) # Environmental parameters + # self.addParameter(QgsProcessingParameterNumber(self.EFFUS_WALL, + # self.tr('Wall type (only with wall scheme)'), QgsProcessingParameterNumber.Integer, + # QVariant(1065), False, minValue=0, maxValue=20000)) self.wallType = ( (self.tr("Brick"), "0"), (self.tr("Concrete"), "1"), @@ -324,7 +377,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.ALBEDO_WALLS, self.tr("Albedo (walls)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.20), False, minValue=0, @@ -335,7 +388,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.ALBEDO_GROUND, self.tr("Albedo (ground)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.15), False, minValue=0, @@ -346,7 +399,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.EMIS_WALLS, self.tr("Emissivity (walls)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.90), False, minValue=0, @@ -357,7 +410,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.EMIS_GROUND, self.tr("Emissivity (ground)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.95), False, minValue=0, @@ -370,7 +423,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.ABS_S, self.tr("Absorption of shortwave radiation of human body"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.70), False, minValue=0, @@ -381,7 +434,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.ABS_L, self.tr("Absorption of longwave radiation of human body"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.95), False, minValue=0, @@ -425,7 +478,7 @@ def initAlgorithm(self, config): QgsProcessingParameterNumber( self.UTC, self.tr("Coordinated Universal Time (UTC) "), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(0), False, minValue=-12, @@ -439,12 +492,11 @@ def initAlgorithm(self, config): self.tr( "Vector point file including Wall of Interest(s) for output with wall surface temperatures" ), - [QgsProcessing.SourceType.TypeVectorPoint], + [QgsProcessing.TypeVectorPoint], optional=True, ) woifile.setFlags( - woifile.flags() - | QgsProcessingParameterDefinition.Flag.FlagAdvanced + woifile.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(woifile) woi_field = QgsProcessingParameterField( @@ -452,27 +504,29 @@ def initAlgorithm(self, config): self.tr("Wall ID field"), "", self.WOI_FILE, - QgsProcessingParameterField.DataType.Numeric, + QgsProcessingParameterField.Numeric, optional=True, ) woi_field.setFlags( - woi_field.flags() - | QgsProcessingParameterDefinition.Flag.FlagAdvanced + woi_field.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(woi_field) # POIs for thermal comfort estimations + # poi = QgsProcessingParameterBoolean(self.POI, + # self.tr("Include Point of Interest(s) for thermal comfort calculations (PET and UTCI)"), defaultValue=False) + # poi.setFlags(poi.flags() | QgsProcessingParameterDefinition.FlagAdvanced) + # self.addParameter(poi) poifile = QgsProcessingParameterFeatureSource( self.POI_FILE, self.tr( "Vector point file including Point of Interest(s) for thermal comfort calculations (PET and UTCI)" ), - [QgsProcessing.SourceType.TypeVectorPoint], + [QgsProcessing.TypeVectorPoint], optional=True, ) poifile.setFlags( - poifile.flags() - | QgsProcessingParameterDefinition.Flag.FlagAdvanced + poifile.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(poifile) poi_field = QgsProcessingParameterField( @@ -480,12 +534,11 @@ def initAlgorithm(self, config): self.tr("ID field"), "", self.POI_FILE, - QgsProcessingParameterField.DataType.Numeric, + QgsProcessingParameterField.Numeric, optional=True, ) poi_field.setFlags( - poi_field.flags() - | QgsProcessingParameterDefinition.Flag.FlagAdvanced + poi_field.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(poi_field) @@ -493,66 +546,66 @@ def initAlgorithm(self, config): age = QgsProcessingParameterNumber( self.AGE, self.tr("Age (yy)"), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(35), optional=True, minValue=0, maxValue=120, ) age.setFlags( - age.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + age.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(age) act = QgsProcessingParameterNumber( self.ACTIVITY, self.tr("Activity (W)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(80), optional=True, minValue=0, maxValue=1000, ) act.setFlags( - act.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + act.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(act) clo = QgsProcessingParameterNumber( self.CLO, self.tr("Clothing (clo)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(0.9), optional=True, minValue=0, maxValue=10, ) clo.setFlags( - clo.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + clo.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(clo) wei = QgsProcessingParameterNumber( self.WEIGHT, self.tr("Weight (kg)"), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(75), optional=True, minValue=0, maxValue=500, ) wei.setFlags( - wei.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + wei.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(wei) hei = QgsProcessingParameterNumber( self.HEIGHT, self.tr("Height (cm)"), - QgsProcessingParameterNumber.Type.Integer, + QgsProcessingParameterNumber.Integer, QVariant(180), optional=True, minValue=0, maxValue=250, ) hei.setFlags( - hei.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + hei.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(hei) sex = QgsProcessingParameterEnum( @@ -563,20 +616,20 @@ def initAlgorithm(self, config): defaultValue=0, ) sex.setFlags( - sex.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + sex.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(sex) shei = QgsProcessingParameterNumber( self.SENSOR_HEIGHT, self.tr("Height of wind sensor (m agl)"), - QgsProcessingParameterNumber.Type.Double, + QgsProcessingParameterNumber.Double, QVariant(10), optional=True, minValue=0, maxValue=250, ) shei.setFlags( - shei.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced + shei.flags() | QgsProcessingParameterDefinition.FlagAdvanced ) self.addParameter(shei) @@ -701,6 +754,12 @@ def processAlgorithm(self, parameters, context, feedback): folderPathPerez = self.parameterAsString( parameters, self.INPUT_ANISO, context ) + useGroundScheme = self.parameterAsString( + parameters, self.USE_GROUNDSCHEME, context + ) + groundTempFile = self.parameterAsString( + parameters, self.INPUT_GROUNDSCHEME, context + ) folderWallScheme = self.parameterAsString( parameters, self.INPUT_WALLSCHEME, context ) @@ -724,6 +783,9 @@ def processAlgorithm(self, parameters, context, feedback): pos = self.parameterAsInt(parameters, self.POSTURE, context) cyl = self.parameterAsBool(parameters, self.CYL, context) + # cyl = 1 + # else: + # cyl = 0 if pos == 0: Fside = 0.22 @@ -772,8 +834,7 @@ def processAlgorithm(self, parameters, context, feedback): outputKdiff = False # outputSstr = False - # If "Save necessary rasters for TreePlanter tool" is ticked, save the - # following raster for TreePlanter or Spatial TC + # If "Save necessary rasters for TreePlanter tool" is ticked, save the following raster for TreePlanter or Spatial TC if outputTreeplanter: outputTmrt = True outputKup = True @@ -825,9 +886,10 @@ def processAlgorithm(self, parameters, context, feedback): # response to issue #85 nd = gdal_dsm.GetRasterBand(1).GetNoDataValue() + # dsm[dsm == nd] = 0. if dsm.min() < 0: dsmraise = np.abs(dsm.min()) - # dsm = dsm + dsmraise #moved to Solweig_run + # dsm = dsm + dsmraise feedback.setProgressText( "Digital Surface Model (DSM) included negative values. DSM raised with " + str(dsmraise) @@ -983,7 +1045,7 @@ def processAlgorithm(self, parameters, context, feedback): try: dataSet = gdal.Open(self.temp_dir + "/svf.tif") svf = dataSet.ReadAsArray().astype(float) - except BaseException: + except: raise QgsProcessingException( "SVF import error: The zipfile including the SVFs seems corrupt. Retry calcualting the SVFs in the Pre-processor or choose another file." ) @@ -1028,7 +1090,7 @@ def processAlgorithm(self, parameters, context, feedback): self.metdata = np.loadtxt( inputMet, skiprows=headernum, delimiter=delim ) - except BaseException: + except: raise QgsProcessingException( "Error: Make sure format of meteorological file is correct. You can" "prepare your data by using 'Prepare Existing Data' in " @@ -1119,6 +1181,16 @@ def processAlgorithm(self, parameters, context, feedback): feedback.setProgressText("Isotropic sky") anisotropic_sky = 0 + # Ground cover scheme + if useGroundScheme: + feedback.setProgressText( + "The ground cover scheme from v2026 is activated" + ) + else: + feedback.setProgressText( + "The ground cover scheme described in 2016 is activated" + ) + # % Ts parameterisation maps if landcover == 1.0: if folderWallScheme: @@ -1128,15 +1200,15 @@ def processAlgorithm(self, parameters, context, feedback): np.max(unique_landcover) > 7 or np.min(unique_landcover) < 1 ): - feedback.pushWarning( - "The land cover grid includes integer values higher (or lower) than standard UMEP-formatted. " - "Land cover grid (should be integer between 1 and 7). If other LC-classes should be included they also need to be included in landcoverclasses_2016a.txt" + raise QgsProcessingException( + "The land cover grid includes integer values higher (or lower) than UMEP-formatted " + "land cover grid (should be integer between 1 and 7). If other LC-classes should be included they also need to be included in landcoverclasses_2016a.txt" ) else: if np.max(lcgrid) > 7 or np.min(lcgrid) < 1: - feedback.pushWarning( - "The land cover grid includes integer values higher (or lower) than standard UMEP-formatted. " - "Land cover grid (should be integer between 1 and 7). If other LC-classes should be included they also need to be included in landcoverclasses_2016a.txt" + raise QgsProcessingException( + "The land cover grid includes integer values higher (or lower) than UMEP-formatted " + "land cover grid (should be integer between 1 and 7). If other LC-classes should be included they also need to be included in landcoverclasses_2016a.txt" ) if np.where(lcgrid) == 3 or np.where(lcgrid) == 4: raise QgsProcessingException( @@ -1171,22 +1243,29 @@ def processAlgorithm(self, parameters, context, feedback): "working_dir": self.temp_dir, "para_json_path": outputDir + "/solweig_parameters.json", "filepath_dsm": filepath_dsm, - "filepath_cdsm": filePath_cdsm, # vegdsm.dataProvider().dataSourceUri() if vegdsm.dataProvider().dataSourceUri() is not None else '', # if vegdsm is None: str(vegdsm.dataProvider().dataSourceUri()), - # str(vegdsm.dataProvider().dataSourceUri()) if vegdsm is None else - # '', #str(vegdsm2.dataProvider().dataSourceUri()), - "filepath_tdsm": filePath_tdsm, + "filepath_cdsm": ( + filePath_cdsm + ), # vegdsm.dataProvider().dataSourceUri() if vegdsm.dataProvider().dataSourceUri() is not None else '', # if vegdsm is None: str(vegdsm.dataProvider().dataSourceUri()), + "filepath_tdsm": ( + filePath_tdsm + ), # str(vegdsm.dataProvider().dataSourceUri()) if vegdsm is None else '', #str(vegdsm2.dataProvider().dataSourceUri()), "filepath_dem": filepath_dem, - # 'C:/Users/xlinfr/Documents/PythonScripts/SOLWEIG/SOLWEIGdata/landcover.tif', - "filepath_lc": filepath_lc, - # 'C:\\Users\\xlinfr\\Desktop\\SOLWEIGdata\\wallheight.tif', - "filepath_wh": filepath_wh, - # 'C:\\Users\\xlinfr\\Desktop\\SOLWEIGdata\\wallaspect.tif', - "filepath_wa": filepath_wa, + "filepath_lc": ( + filepath_lc + ), # 'C:/Users/xlinfr/Documents/PythonScripts/SOLWEIG/SOLWEIGdata/landcover.tif', + "filepath_wh": ( + filepath_wh + ), #'C:\\Users\\xlinfr\\Desktop\\SOLWEIGdata\\wallheight.tif', + "filepath_wa": ( + filepath_wa + ), #'C:\\Users\\xlinfr\\Desktop\\SOLWEIGdata\\wallaspect.tif', "input_svf": inputSVF, - # 'C:\\Users\\xlinfr\\Desktop\\SOLWEIGdata\\shadowmats.npz', - "input_aniso": folderPathPerez, + "input_aniso": ( + folderPathPerez + ), # 'C:\\Users\\xlinfr\\Desktop\\SOLWEIGdata\\shadowmats.npz', "poi_file": poi_file, "poi_field": poi_field, + "input_surf": groundTempFile, "input_wall": folderWallScheme, "woi_file": woi_file, "woi_field": woi_field, @@ -1206,7 +1285,7 @@ def processAlgorithm(self, parameters, context, feedback): "demforbuild": int(demforbuild), "aniso": int(anisotropic_sky), "wallscheme": wallScheme, - "walltype": wall_type, # 'Brick_wall', #:TODO + "walltype": wall_type, #'Brick_wall', #:TODO "outputtmrt": int(outputTmrt), "outputkup": int(outputKup), "outputkdown": int(outputKdown), @@ -1252,7 +1331,7 @@ def displayName(self): Returns the translated algorithm name, which should be used for any user-visible display of the algorithm name. """ - return self.tr("Outdoor Thermal Comfort: SOLWEIG v2025a") + return self.tr("Outdoor Thermal Comfort: SOLWEIG v2026a") def group(self): """ @@ -1273,22 +1352,15 @@ def groupId(self): def shortHelpString(self): return self.tr( - "SOLWEIG (v2025a) is a model which can be used to estimate spatial variations of 3D radiation fluxes and " - "mean radiant temperature (Tmrt) in complex urban settings. The SOLWEIG model originally followed the " + "SOLWEIG (v2026a) is a model which can be used to estimate spatial variations of 3D radiation fluxes and " + "mean radiant temperature (Tmrt) in complex urban settings. The SOLWEIG model follows the same " "approach commonly adopted to observe Tmrt, with shortwave and longwave radiation fluxes from " - "six directions being individually calculated to derive Tmrt. The model can also consider a person as a " - "standing cylinder. The model requires a limited number " - "of inputs, such as global shortwave radiation, air temperature, relative " + "six directions being individually calculated to derive Tmrt. The model requires a limited number " + "of inputs, such as direct, diffuse and global shortwave radiation, air temperature, relative " "humidity, urban geometry and geographical information (latitude, longitude and elevation). " "Additional vegetation and ground cover information can also be used to imporove the estimation of Tmrt.\n" "\n" - "Tools to generate sky view factors, wall height and aspect etc. is available in the pre-processing section in UMEP\n" - "\n" - "NOTE:\n" - "- Anisotrophic sky models for long- and diffuse shortwave are automatically activated when the *.npz file " - "for shadow maps are included.\n" - "- Wall scheme model is automatically activated when the .npz-file for wall voxels are included. This will " - "slow down the model consideraby as SOLWEIG become a near 3D-model.\n" + "Tools to generate sky view factors, wall height and aspect etc. is available in the pre-processing past in UMEP\n" "\n" "------------\n" "\n"