From c30e79e19fc09ce2c06390ded9b34624ab755866 Mon Sep 17 00:00:00 2001 From: Ilay Falach Date: Sun, 17 May 2026 13:24:52 +0300 Subject: [PATCH 1/5] Migrate from unum to pint across simulations, evaporation, LSM, hydrodynamics, and risk assessment (#875) Replace all `from unum.units import *` star imports and bare Unum unit literals (m, kg, cm, s, etc.) with explicit `ureg.` equivalents. Replace `.asNumber()` calls with `.m_as()` and add `unumToPint()` at function boundaries for backwards compatibility with legacy Unum callers. Files migrated: - hera/simulations/gaussian/{gasCloud,toolkit,Meteorology,Sigma,FallingNonEvaporatingDroplets,MeshUtils,DropletCloud}.py - hera/simulations/evaporation/{monaghan,models}.py - hera/simulations/deposition/models.py - hera/simulations/LSM/{toolkit,template,singleSimulation}.py - hera/simulations/hydrodynamics/nearWallFlow.py - hera/riskassessment/agents/Agents.py - hera/riskassessment/agents/effects/InjuryLevel.py - hera/utils/matplotlibCountour.py Co-Authored-By: Claude Sonnet 4.6 --- hera/riskassessment/agents/Agents.py | 17 +- .../agents/effects/InjuryLevel.py | 3 +- hera/simulations/LSM/singleSimulation.py | 6 +- hera/simulations/LSM/template.py | 8 +- hera/simulations/LSM/toolkit.py | 10 +- hera/simulations/deposition/models.py | 17 +- hera/simulations/evaporation/models.py | 22 +- hera/simulations/evaporation/monaghan.py | 15 +- hera/simulations/gaussian/DropletCloud.py | 30 +-- .../gaussian/FallingNonEvaporatingDroplets.py | 127 +++++----- hera/simulations/gaussian/MeshUtils.py | 12 +- hera/simulations/gaussian/Meteorology.py | 55 +++-- hera/simulations/gaussian/Sigma.py | 28 +-- hera/simulations/gaussian/gasCloud.py | 222 +++++++++--------- hera/simulations/gaussian/toolkit.py | 70 +++--- .../simulations/hydrodynamics/nearWallFlow.py | 32 +-- hera/utils/matplotlibCountour.py | 4 +- 17 files changed, 338 insertions(+), 340 deletions(-) diff --git a/hera/riskassessment/agents/Agents.py b/hera/riskassessment/agents/Agents.py index 8689f7f8f..d06f9e830 100644 --- a/hera/riskassessment/agents/Agents.py +++ b/hera/riskassessment/agents/Agents.py @@ -1,4 +1,4 @@ -from hera.utils.unitHandler import ureg, Unum, Quantity, celsius, K +from hera.utils.unitHandler import ureg, unumToPint, Quantity, celsius, K from hera.riskassessment.agents.effects import injuryfactory import numpy import json @@ -254,10 +254,7 @@ def getVolatility(self,temperature): :return: The vapor saturation as Unum. """ - if isinstance(temperature, Unum): - temperature = temperature.asNumber(celsius) - elif isinstance(temperature, Quantity): - temperature = temperature.m_as(ureg.celsius) + temperature = unumToPint(temperature).m_as(ureg.celsius) MW = self.getMolecularWeight().to(ureg.g/ureg.mol) a,b,c,d = self._volatilityConst @@ -277,10 +274,7 @@ def getDensity(self, temperature): :return: The density as Unum """ - if isinstance(temperature, Unum): - temperature = temperature.asNumber(celsius) - elif isinstance(temperature, Quantity): - temperature = temperature.m_as(ureg.celsius) + temperature = unumToPint(temperature).m_as(ureg.celsius) a,b,c = self._densityConst return (a-b*(temperature-c))*ureg.g/ureg.cm**3 @@ -299,10 +293,7 @@ def vaporPressure(self, temperature): float Vapor pressure in bar. """ - if isinstance(temperature, Unum): - temperature = temperature.asNumber(K) - elif isinstance(temperature, Quantity): - temperature = temperature.m_as(ureg.kelvin) + temperature = unumToPint(temperature).m_as(ureg.kelvin) A = self._vaporConst["A"] B = self._vaporConst["B"] C = self._vaporConst["C"] diff --git a/hera/riskassessment/agents/effects/InjuryLevel.py b/hera/riskassessment/agents/effects/InjuryLevel.py index e59949a06..f4631a635 100755 --- a/hera/riskassessment/agents/effects/InjuryLevel.py +++ b/hera/riskassessment/agents/effects/InjuryLevel.py @@ -1,6 +1,7 @@ import geopandas import pandas from hera.utils import tounit, ureg +from hera.utils.unitHandler import unumToPint import numpy import json @@ -338,7 +339,7 @@ def _getGeopandas(self,concentrationField,x,y,**parameters): """ Return the correct geopandas of the Injury level """ - level = self.threshold.asNumber() + level = unumToPint(self.threshold).magnitude CS = plt.contour(concentrationField[x],concentrationField[y],concentrationField.squeeze(),levels=numpy.atleast_1d(level)) if numpy.max(CS.levels) < level: ret = geopandas.GeoDataFrame() diff --git a/hera/simulations/LSM/singleSimulation.py b/hera/simulations/LSM/singleSimulation.py index ef5084555..787c2d993 100644 --- a/hera/simulations/LSM/singleSimulation.py +++ b/hera/simulations/LSM/singleSimulation.py @@ -2,7 +2,7 @@ import xarray import numpy import os -from hera.utils.unitHandler import * +from hera.utils.unitHandler import ureg, unumToPint from ...utils import tounit,tonumber @@ -67,10 +67,6 @@ def getDosage(self, Q=1 * ureg.kg, time_units=ureg.min, q_units=ureg.mg): final_xarray.attrs['dt'] = dt_minutes.to(time_units) final_xarray.attrs['Q'] = Q.to(q_units) final_xarray.attrs['C'] = q_units/ ureg.m ** 3 - if not ret_pint: - final_xarray.attrs['dt'] = pintToUnum(final_xarray.attrs['dt']) - final_xarray.attrs['Q'] = pintToUnum(final_xarray.attrs['Q']) - final_xarray.attrs['C'] = pintToUnum(final_xarray.attrs['C']) Qfactor = (Q.to(q_units) * ureg.min / ureg.m ** 3).m_as(q_units * time_units / ureg.m ** 3) diff --git a/hera/simulations/LSM/template.py b/hera/simulations/LSM/template.py index 9d9db8c6d..7c444c4d2 100644 --- a/hera/simulations/LSM/template.py +++ b/hera/simulations/LSM/template.py @@ -9,7 +9,7 @@ from ..utils.inputForModelsCreation import InputForModelsCreator from hera.simulations.LSM.singleSimulation import SingleSimulation from hera.datalayer import datatypes -from hera.utils.unitHandler import Quantity, ureg, Unum +from hera.utils.unitHandler import Quantity, ureg, unumToPint from hera import toolkit from hera.utils.jsonutils import JSONToConfiguration, stripConfigurationUnits from hera.utils import dictToMongoQuery, get_classMethod_logger @@ -318,10 +318,8 @@ def prepareParams(desc, paramsToPrepare): if desc is not None and 'units' in desc: for key in desc["units"].keys(): param_item= paramsToPrepare[key] - if isinstance(param_item, Unum): - paramsToPrepare[key] = param_item.asNumber(eval(desc["units"][key])) - elif isinstance(param_item, Quantity): - paramsToPrepare[key] = param_item.m_as(desc["units"][key]) + if hasattr(param_item, 'asNumber') or hasattr(param_item, 'magnitude'): + paramsToPrepare[key] = unumToPint(param_item).m_as(ureg.parse_expression(desc["units"][key])) else: raise ValueError(f"parameters must use either pint or unum to specify units, currently type({param_item})={type(param_item)}") if 'duration' in paramsToPrepare and isinstance(paramsToPrepare['duration'], Quantity): diff --git a/hera/simulations/LSM/toolkit.py b/hera/simulations/LSM/toolkit.py index 250cecc5d..5d1ec4690 100644 --- a/hera/simulations/LSM/toolkit.py +++ b/hera/simulations/LSM/toolkit.py @@ -11,7 +11,7 @@ from ... import datalayer from ... import toolkit from .singleSimulation import SingleSimulation -from unum.units import * +from hera.utils.unitHandler import ureg, unumToPint from ..utils.coordinateHandler import coordinateHandler from pathlib import Path @@ -238,10 +238,10 @@ def getSimulations(self,simulationName=None,unitsTemplateVersion="v4-general", * for key in template._document['desc']["units"].keys(): if key in query.keys(): query_item= query[key] - if isinstance(query_item, Unum): - query[key] = query_item.asNumber(eval(template._document['desc']["units"][key])) - elif isinstance(query_item, Quantity): - query[key] = query_item.m_as(template._document['desc']["units"][key]) + if hasattr(query_item, 'asNumber') or hasattr(query_item, 'magnitude'): + pint_item = unumToPint(query_item) + unit_expr = template._document['desc']["units"][key] + query[key] = pint_item.m_as(ureg.parse_expression(unit_expr)) else: raise ValueError(f"query must use either pint or unum to specify units, currently type({query[key]})={type(query[key])}") else: diff --git a/hera/simulations/deposition/models.py b/hera/simulations/deposition/models.py index dff8c1a3a..a14b07df8 100644 --- a/hera/simulations/deposition/models.py +++ b/hera/simulations/deposition/models.py @@ -1,7 +1,7 @@ from ...datalayer import project from ...riskassessment import AgentHome from ..utils import tonumber, tounit -from unum.units import * +from hera.utils.unitHandler import ureg, unumToPint import numpy class depositionModels(object): @@ -71,16 +71,23 @@ def depositionModel(self,newModel): self._depositionModel=newModel def __init__(self, projectName = "deposition", surface={"name":"Desert","zrough":0.04}, - ustar=1, density=1500, diameter = 1E-6*m,heatFlux=0.1*W/(m**2), temperature=293*K, depositionModel="Petroff"): + ustar=1, density=1500, diameter=None, heatFlux=None, temperature=None, depositionModel="Petroff"): + + if diameter is None: + diameter = 1E-6*ureg.m + if heatFlux is None: + heatFlux = 0.1*ureg.W/ureg.m**2 + if temperature is None: + temperature = 293*ureg.K p = project.Project(projectName) self._depositionModel = depositionModel self._ustar = ustar self._density = density - self._temperature = tonumber(tounit(temperature,K),K) + self._temperature = tonumber(tounit(temperature, ureg.K), ureg.K) self._surface = surface if type(surface)==dict else p.getCacheDocuments(type="surface",surface=surface)[0].asDict()["desc"] - self._diameter = tonumber(tounit(diameter,m),m) - self._heatFlux = tonumber(tounit(heatFlux, W/(m**2)), W/(m**2)) + self._diameter = tonumber(tounit(diameter, ureg.m), ureg.m) + self._heatFlux = tonumber(tounit(heatFlux, ureg.W/ureg.m**2), ureg.W/ureg.m**2) def depositionRate(self): return getattr(self, f"depositionRate_{self._depositionModel}")() diff --git a/hera/simulations/evaporation/models.py b/hera/simulations/evaporation/models.py index 9a965d837..c763089ad 100644 --- a/hera/simulations/evaporation/models.py +++ b/hera/simulations/evaporation/models.py @@ -1,7 +1,7 @@ from ...datalayer import project from ...riskassessment import RiskToolkit from ..utils import tonumber, tounit -from unum.units import * +from hera.utils.unitHandler import ureg, unumToPint import numpy class evaporationModels(object): @@ -60,7 +60,7 @@ def Vair(self): @property def Magent(self): - return self._agent.physicalproperties.molecularWeight.asNumber(g/mol) + return unumToPint(self._agent.physicalproperties.molecularWeight).m_as(ureg.g/ureg.mol) @property def Vagent(self): @@ -75,14 +75,14 @@ def __init__(self,agent, evaporationModel="US",dinamicViscocityModel="powerLaw", self._dinamicViscocityModel = dinamicViscocityModel self._evaporationModel = evaporationModel try: - self._Vagent = self._agent.physicalproperties.molecularVolume.asNumber(cm ** 3 / mol) # agent.physicalproperties.molecularVolume.asNumber(ml/mol) + self._Vagent = unumToPint(self._agent.physicalproperties.molecularVolume).m_as(ureg.cm**3/ureg.mol) # agent.physicalproperties.molecularVolume.asNumber(ml/mol) except: print("Note that the molecular volume is not given, therefore, the molecular diffusion is set to EPA.") self._molecularDiffusionModel = "EPA" self._Vagent = None def molecularDiffusion(self,temperature): - temperature = tonumber(temperature, K) + temperature = tonumber(temperature, ureg.K) return getattr(self,f"molecularDiffusion_{self._molecularDiffusionModel}")(temperature) def molecularDiffusion_FSG(self,temperature): @@ -92,7 +92,7 @@ def molecularDiffusion_EPA(self,temperature): return 0.00000000409*(temperature**1.9)*numpy.sqrt(1/self.Mair+1/self.Magent)/numpy.cbrt(self.Magent) def dynamicViscocityAir(self,temperature): - temperature = tonumber(temperature, K) + temperature = tonumber(temperature, ureg.K) return getattr(self,f"dynamicViscocityAir_{self._dinamicViscocityModel}")(temperature) def dynamicViscocityAir_powerLaw(self,temperature): @@ -106,13 +106,15 @@ def Schmidt(self,temperature): density = self.Mair / (temperature * 8.205 * 10 ** (-5)) # (g/(m^3)) return self.dynamicViscocityAir(temperature=temperature)/(density*self.molecularDiffusion(temperature)) - def flux(self,diameter,velocity,temperature,units=g/(m**2*s)): - return tounit(getattr(self, f"flux_{self._evaporationModel}")(diameter,velocity,temperature)*g/(m**2*s),units) + def flux(self,diameter,velocity,temperature,units=None): + if units is None: + units = ureg.g/(ureg.m**2*ureg.s) + return tounit(getattr(self, f"flux_{self._evaporationModel}")(diameter,velocity,temperature)*ureg.g/(ureg.m**2*ureg.s), unumToPint(units)) def flux_US(self, diameter,velocity,temperature): - temperature = tonumber(temperature, K) - diameter = tonumber(diameter, m) - velocity = tonumber(velocity,m/s) + temperature = tonumber(temperature, ureg.K) + diameter = tonumber(diameter, ureg.m) + velocity = tonumber(velocity, ureg.m/ureg.s) Re = self.Reynolds(diameter=diameter,velocity=velocity,temperature=temperature) Sc = self.Schmidt(temperature=temperature) Km = 0.664*(Re**(-0.5))*(Sc**(-2/3))*velocity if Re<=20000 else 0.0366*(Re**(-0.2))*(Sc**(-2/3))*velocity diff --git a/hera/simulations/evaporation/monaghan.py b/hera/simulations/evaporation/monaghan.py index 38efeb9b1..5f022e99c 100644 --- a/hera/simulations/evaporation/monaghan.py +++ b/hera/simulations/evaporation/monaghan.py @@ -1,10 +1,9 @@ import pandas import numpy -from unum.units import * +from hera.utils.unitHandler import ureg, unumToPint, Quantity from ..gaussian.Meteorology import StandardMeteorolgyConstant_powerLaw from ..utils import tounit,tonumber from pyriskassessment.agents.Agents import Agent -from unum import Unum class MonaghanConstantConditions(object): """ @@ -34,7 +33,7 @@ def agent(self,value): @property def P0(self): - return 5.0/60*cm/s + return 5.0/60*ureg.cm/ureg.s @property def U2P1(self): @@ -46,7 +45,7 @@ def MASSRATIO(self): @property def WINDHEIGHT(self): - return 200*cm + return 200*ureg.cm @property @@ -105,13 +104,13 @@ def Q(self,timelist,dropletRadius): evaporation_nc[2] = (sorption_p+self.P0) / ((sorption_p+p1)*self.MASSRATIO) # equivalent to variable nc2 in the code of Hezi, but not multiplied by the droplet radius. evaporation_nc[3] = evaporation_nc[2]*evaporation_nc[1] # equivalent to variable nc3 in the code of Hezi - evaporation = [Unum.coercetounit(x).asNumber() for x in evaporation] - evaporation_nc = [Unum.coercetounit(x) for x in evaporation_nc] + evaporation = [x.magnitude if hasattr(x, 'magnitude') else float(x) for x in evaporation] + evaporation_nc = [unumToPint(x) if hasattr(x, 'asNumber') else x for x in evaporation_nc] ret = numpy.zeros(len(timelist)) - sstime = (2*tounit(dropletRadius,um)*evaporation_nc[3]).asNumber(s) - tau = (2 * tounit(dropletRadius, um) * evaporation_nc[1]).asNumber(s) + sstime = (2*tounit(dropletRadius,ureg.um)*evaporation_nc[3]).m_as(ureg.s) + tau = (2 * tounit(dropletRadius, ureg.um) * evaporation_nc[1]).m_as(ureg.s) sstime_timedelta = pandas.to_timedelta(sstime,"s") tau_timedelta = pandas.to_timedelta(tau, "s") diff --git a/hera/simulations/gaussian/DropletCloud.py b/hera/simulations/gaussian/DropletCloud.py index 64b48bba9..a4e18b60d 100644 --- a/hera/simulations/gaussian/DropletCloud.py +++ b/hera/simulations/gaussian/DropletCloud.py @@ -1,4 +1,4 @@ -from hera.utils.unitHandler import * +from hera.utils.unitHandler import ureg, unumToPint import pandas import numpy from scipy.stats import lognorm @@ -56,7 +56,7 @@ def getGround(self,T): print(f"solving droplets {i}/{total}" ,end="\r") res = droplet.solveToTime(T) res['N'] = droplet.N - res['DropletArea'] = droplet.AreaOnSurface.asNumber(um**2) + res['DropletArea'] = unumToPint(droplet.AreaOnSurface).m_as(ureg.um**2) retList.append(res.iloc[-1].to_frame().T) ret = pandas.concat(retList, ignore_index=True) @@ -65,22 +65,22 @@ def getGround(self,T): def _initDropletPosition(self, mmd, geometricstd, position, Q, clouds=30, meteorologyname="logNormal",met_kwargs={}, **kwargs): - rv = lognorm(numpy.log(geometricstd), scale=tonumber(mmd, m)) + rv = lognorm(numpy.log(geometricstd), scale=tonumber(mmd, ureg.m)) lower = rv.ppf(1e-4) upper = rv.ppf(1-1e-4) - met_params = dict(z0=10*cm) + met_params = dict(z0=10*ureg.cm) met_params.update(met_kwargs) interval = numpy.logspace(numpy.log(lower),numpy.log(upper),clouds,base=numpy.e) dh = numpy.diff(numpy.log(interval))[0] - massFractionVector = numpy.diff(rv.cdf(interval))*tonumber(Q,kg) + massFractionVector = numpy.diff(rv.cdf(interval))*tonumber(Q,ureg.kg) diameterVector = numpy.exp(numpy.log(interval[:-1])+dh/2.) # [m] for dropletDiam,dropletQ in zip(diameterVector,massFractionVector): - droplets = FallingNonEvaporatingDroplets(particleDiameter=dropletDiam*m, - Q=dropletQ*kg, + droplets = FallingNonEvaporatingDroplets(particleDiameter=dropletDiam*ureg.m, + Q=dropletQ*ureg.kg, position=position, meteorologyName=meteorologyname, met_kwargs=met_params, @@ -120,8 +120,8 @@ def __init__(self, mmd, geometricstd, position, Q, linelength, clouds=30,linepos self._dropletList = [] qCloud = Q/linepositions - for Ypos in numpy.linspace(0,tonumber(linelength,m),linepositions): - curpos = (position[0],tounit(Ypos,m),position[2]) + for Ypos in numpy.linspace(0,tonumber(linelength,ureg.m),linepositions): + curpos = (position[0],tounit(Ypos,ureg.m),position[2]) self._initDropletPosition(mmd, geometricstd, curpos, qCloud, clouds=clouds,meteorologyname=meteorologyname, **kwargs) @@ -158,9 +158,9 @@ def __init__(self,mmd,geometricstd,position,Q,clippedDiameter,clouds=30, meteoro def _initDropletPosition(self, mmd, geometricstd, position, Q, clouds=30, meteorologyname="StandardMeteorolgyConstant", **kwargs): - clippedDiameter = tonumber(self.clippedDiameter,m) + clippedDiameter = tonumber(self.clippedDiameter,ureg.m) - rv = lognorm(numpy.log(geometricstd), scale=tonumber(mmd, m)) + rv = lognorm(numpy.log(geometricstd), scale=tonumber(mmd, ureg.m)) lower = rv.ppf(1e-4) upper = clippedDiameter @@ -169,13 +169,13 @@ def _initDropletPosition(self, mmd, geometricstd, position, Q, clouds=30, meteor maxMass = rv.cdf(interval[-1]) - massFractionVector = numpy.diff(rv.cdf(interval))*tonumber(Q,kg) + massFractionVector = numpy.diff(rv.cdf(interval))*tonumber(Q,ureg.kg) diameterVector = numpy.exp(numpy.log(interval[:-1])+dh/2.) # [m] massFractionVector /= maxMass for dropletDiam,dropletQ in zip(diameterVector,massFractionVector): - droplets = FallingNonEvaporatingDroplets(particleDiameter=dropletDiam*m, Q=dropletQ*kg, position=position, meteorologyName=meteorologyname, **kwargs) + droplets = FallingNonEvaporatingDroplets(particleDiameter=dropletDiam*ureg.m, Q=dropletQ*ureg.kg, position=position, meteorologyName=meteorologyname, **kwargs) self._dropletList.append(droplets) @@ -184,7 +184,7 @@ class CirclePositionClippedDropletsCloud(FixedPointClippedDropletCloud): """ A line across the wind. """ - def __init__(self, mmd, geometricstd, position, Q, clippedDiameter, clouds=30,radius=10*m,circlepositions=4, meteorologyname="StandardMeteorolgyConstant", **kwargs): + def __init__(self, mmd, geometricstd, position, Q, clippedDiameter, clouds=30,radius=10*ureg.m,circlepositions=4, meteorologyname="StandardMeteorolgyConstant", **kwargs): """ Creates a list of clouds (discretization according to the number of clouds). @@ -212,7 +212,7 @@ def __init__(self, mmd, geometricstd, position, Q, clippedDiameter, clouds=30,ra self._dropletList = [] self.clippedDiameter = clippedDiameter - radius_m = tounit(radius,m) + radius_m = tounit(radius,ureg.m) qCloud = Q/circlepositions for deg in numpy.linspace(0,2*numpy.pi,circlepositions): diff --git a/hera/simulations/gaussian/FallingNonEvaporatingDroplets.py b/hera/simulations/gaussian/FallingNonEvaporatingDroplets.py index c04adf5ac..9ef7602e5 100644 --- a/hera/simulations/gaussian/FallingNonEvaporatingDroplets.py +++ b/hera/simulations/gaussian/FallingNonEvaporatingDroplets.py @@ -5,7 +5,7 @@ except ImportError: print("wrong scipy, can't solve the system") -from hera.utils.unitHandler import * +from hera.utils.unitHandler import ureg, unumToPint, tounit, tonumber from scipy.optimize import root from scipy.constants import g as gconst from hera.simulations.gaussian.Meteorology import MeteorologyFactory @@ -50,7 +50,7 @@ def position(self,value): except: raise ValueError("Must be a 3 tuple. ") - self._position=[tounit(x,m) for x in value] + self._position=[tounit(x,ureg.m) for x in value] @property def cloudSigma(self): @@ -64,7 +64,7 @@ def cloudSigma(self, value): except: raise ValueError("Must be a 3 tuple. ") - self._cloudSigma = [tounit(x, m) for x in value] + self._cloudSigma = [tounit(x, ureg.m) for x in value] @property def cloudSigmaX(self): @@ -83,7 +83,7 @@ def Q(self): return self._cloudQ @Q.setter def Q(self,value): - self._cloudQ = tounit(value,kg) + self._cloudQ = tounit(value,ureg.kg) @property def particleMass(self): @@ -111,7 +111,7 @@ def particleDiameter(self): @particleDiameter.setter def particleDiameter(self,value): - self._particleDiameter = tounit(value,um) + self._particleDiameter = tounit(value,ureg.um) self.__setParticleMass() @property @@ -120,47 +120,47 @@ def rho_p(self): @rho_p.setter def rho_p(self, value): - self._rho_p = tounit(value, g / cm ** 3) + self._rho_p = tounit(value, ureg.g / ureg.cm ** 3) self.__setParticleMass() - @property - def N(self): + @property + def N(self): """ - The total number of particles. + The total number of particles. """ - return (self.Q/self.particleMass).asNumber() + return (self.Q/self.particleMass).magnitude @property - def SpreadFactor(self): + def SpreadFactor(self): + """ + The spread factor of the droplets. Assuming oil and taken from the VX parameters. """ - The spread factor of the droplets. Assuming oil and taken from the VX parameters. - """ return 4.5 @property def g(self): - return gconst*m/s**2 + return gconst*ureg.m/ureg.s**2 @property - def AreaOnSurface(self): + def AreaOnSurface(self): """ - The total surface area that the droplet occupies after its spread. - Q 3 * Q * SpreadFactor - ------ * pi*(d/2)**2*SpreadFactor = --------------------- + The total surface area that the droplet occupies after its spread. + Q 3 * Q * SpreadFactor + ------ * pi*(d/2)**2*SpreadFactor = --------------------- Mass 2 * rho_l*d - + """ return (self.Q*self.SpreadFactor*numpy.pi*(self.particleDiameter/2.)**2)/self.particleMass def __setParticleMass(self): if (self.rho_p is None or self.particleDiameter is None): return - self._particleMass = (1 / 6. * self.rho_p * numpy.pi * self.particleDiameter ** 3).asUnit(kg) + self._particleMass = (1 / 6. * self.rho_p * numpy.pi * self.particleDiameter ** 3).to(ureg.kg) def __init__(self, particleDiameter, - Q=1*kg, - position=(0*m,0*m,0*m), + Q=1*ureg.kg, + position=(0*ureg.m,0*ureg.m,0*ureg.m), met_kwargs={}, meteorologyName="StandardMeteorolgyConstant", **kwargs): @@ -168,12 +168,12 @@ def __init__(self, Initializes the particle cloud. Currently initialized to point source. - See Aroesty - Atmospheric diffusion of droplet clouds for details. + See Aroesty - Atmospheric diffusion of droplet clouds for details. - Plume (Aroesty pg 19): Plume diffusion refers to circumstances where material release and smaplling times are long compared with travel time from the source. - Puff : material release and samplling times are short compared with travel time. + Plume (Aroesty pg 19): Plume diffusion refers to circumstances where material release and smaplling times are long compared with travel time from the source. + Puff : material release and samplling times are short compared with travel time. - Since the fall even from 500m is O(min) ~ O(release time) than we are always in the plume regime. + Since the fall even from 500m is O(min) ~ O(release time) than we are always in the plume regime. :param particleDiameter: @@ -186,21 +186,21 @@ def __init__(self, The initial cloud dimensions. :param position: - A 3 tuple of unum length that indicates the initial position of the cloud in space. + A 3 tuple of pint Quantity length that indicates the initial position of the cloud in space. default unit [m]. :param meteorology: The name meteorology class to use met_kwargs : the parameters or the meteorology. - :param dispersionType: + :param dispersionType: The dispersion is a plume or a puff. - - The differance is taken from Aerosty - Atmospheric diffusion of droplet clouds. - - - - According to Aerosty if it is a plume (release time + + The differance is taken from Aerosty - Atmospheric diffusion of droplet clouds. + + + + According to Aerosty if it is a plume (release time :param kwargs: * all Passed to the meteorology factory. @@ -215,7 +215,7 @@ def __init__(self, * rho_l : The density of the liquid. """ - cloudSigma = (0 * m, 0 * m, 0 * m) + cloudSigma = (0 * ureg.m, 0 * ureg.m, 0 * ureg.m) self._meteorology = MeteorologyFactory().getMeteorology(meteorologyName,**met_kwargs) dragCoeffFunc = kwargs.get("dragCoeffFunc","Haugen") @@ -225,14 +225,14 @@ def __init__(self, self._correctionfunc = getattr(self,"correctionCloud_%s" % cloudCorrectionFunc.title()) self.particleDiameter = particleDiameter - self.rho_p = kwargs.get("rho_l", 0.9 * g / cm ** 3) # oil + self.rho_p = kwargs.get("rho_l", 0.9 * ureg.g / ureg.cm ** 3) # oil self.Q = Q self.cloudSigma = cloudSigma self.position = position def getTerminalVelocity(self,height=None): Res = root(self._VTFunc,1,args=(height,)) - return Res.x[0]*m/s + return Res.x[0]*ureg.m/ureg.s def _DragCoefficient_Ik(self,Re): @@ -296,7 +296,7 @@ def _DragCoefficient_Fan(self,Re): :param Re: :return: - Drag coefficient + Drag coefficient. """ if Re < 2: Cd = 24 / Re @@ -322,30 +322,29 @@ def _DragCoefficient_Haugen(self,Re): def _VTFunc(self,Vt,height=None): - height = tounit(height if height is not None else self.z,m) - + height = tounit(height if height is not None else self.z, ureg.m) rho_air = self.meteorology.getAirDensity(height) nu_air = self.meteorology.getAirDynamicViscosity(height) - uVt = Vt[0]*m/s - Re = (uVt*self.particleDiameter/nu_air).asNumber() + uVt = Vt[0]*ureg.m/ureg.s + Re = (uVt*self.particleDiameter/nu_air).magnitude CD = self._dragfunc(Re) - GuessVt = 4 * gconst * (m/s**2) * (self.rho_p - rho_air) * self.particleDiameter ** 2 / (3 * nu_air * CD * Re) + GuessVt = 4 * gconst * (ureg.m/ureg.s**2) * (self.rho_p - rho_air) * self.particleDiameter ** 2 / (3 * nu_air * CD * Re) - return (GuessVt - uVt).asNumber(m/s) + return (GuessVt - uVt).m_as(ureg.m/ureg.s) def correctionCloud_None(self): return 1 def correctionCloud_Plume(self): - Ubar = numpy.mean([self.meteorology.getWindVelocity(z) for z in numpy.arange(0, self.z.asNumber(m))]) - return ((1 + (self.beta * self.getTerminalVelocity() / Ubar) ** 2) ** -0.25).asNumber() # taking beta=5 + Ubar = numpy.mean([self.meteorology.getWindVelocity(z) for z in numpy.arange(0, unumToPint(self.z).m_as(ureg.m))]) + return ((1 + (self.beta * self.getTerminalVelocity() / Ubar) ** 2) ** -0.25).magnitude # taking beta=5 def correctionCloud_Puff(self): - Ubar = numpy.mean([self.meteorology.getWindVelocity(z) for z in numpy.arange(0, self.z.asNumber(m))]) - return ((1 + (self.beta * self.getTerminalVelocity() / Ubar) ** 2) ** -0.5).asNumber() # taking beta=5 + Ubar = numpy.mean([self.meteorology.getWindVelocity(z) for z in numpy.arange(0, unumToPint(self.z).m_as(ureg.m))]) + return ((1 + (self.beta * self.getTerminalVelocity() / Ubar) ** 2) ** -0.5).magnitude # taking beta=5 def solveToTime(self, T): @@ -362,9 +361,9 @@ def solveToTime(self, T): - particle velocity """ # y0: x,y,z,dist,u,w - y0 = numpy.array([self.x.asNumber(m), - self.y.asNumber(m), - self.z.asNumber(m), + y0 = numpy.array([unumToPint(self.x).m_as(ureg.m), + unumToPint(self.y).m_as(ureg.m), + unumToPint(self.z).m_as(ureg.m), 0, 0, 0]) @@ -379,8 +378,8 @@ def solveToTime(self, T): correction = self._correctionfunc() ret = ret.assign(sigmaXCorrected=ret['sigmaX']*correction)\ .assign(sigmaZCorrected=ret['sigmaZ']*correction)\ - .assign(Q=self.Q.asNumber(kg))\ - .assign(diameter=self.particleDiameter.asNumber(um)) + .assign(Q=unumToPint(self.Q).m_as(ureg.kg))\ + .assign(diameter=unumToPint(self.particleDiameter).m_as(ureg.um)) return ret @@ -401,10 +400,10 @@ def _fallingParticle(self, t, y): x, yc, z, dist, u_p, w_p = y z = numpy.max(z,0) - z = tounit(z,m) + z = tounit(z, ureg.m) - u_p = tounit(u_p,m/s) - w_p = tounit(w_p, m / s) + u_p = tounit(u_p, ureg.m/ureg.s) + w_p = tounit(w_p, ureg.m / ureg.s) try: U = self.meteorology.getWindVelocity(z) @@ -420,7 +419,7 @@ def _fallingParticle(self, t, y): Uparticle = (u_p**2 + w_p**2)**0.5 # 2. Calculate Re - Re = (rho_air*Uabs * self.particleDiameter / nu_air).asNumber() + Re = (rho_air*Uabs * self.particleDiameter / nu_air).magnitude # 3. Calculate Cd Cd = self._dragfunc(Re) @@ -433,12 +432,12 @@ def _fallingParticle(self, t, y): new_w_p = -Cd*Coeff*Uabs*(w_p)- self.g newy = numpy.zeros(y.shape) - newy[0] = u_p.asNumber(m/s) # x - newy[1] = 0 # y - newy[2] = w_p.asNumber(m/s) # z - newy[3] = Uparticle.asNumber(m/s) # dist - distance, - newy[4] = new_u_p.asNumber() # u_p - newy[5] = new_w_p.asNumber() # w_p + newy[0] = u_p.m_as(ureg.m/ureg.s) # x + newy[1] = 0 # y + newy[2] = w_p.m_as(ureg.m/ureg.s) # z + newy[3] = Uparticle.m_as(ureg.m/ureg.s) # dist - distance, + newy[4] = new_u_p.magnitude # u_p + newy[5] = new_w_p.magnitude # w_p return newy class hit_ground(object): @@ -449,5 +448,3 @@ def __call__(self, t, y): #particleDiameter, Q = 1 * kg, cloudSigma = (1 * m, 1 * m, 1 * m), position = (0 * m, 0 * m,0 * m), meteorologyName = "StandardMeteorolgyConstant", ** kwargs): #particle = FallingNonEvaporatingDroplets(particleDiameter=3*mm,position = (0 * m, 0 * m,300 * m)) - - diff --git a/hera/simulations/gaussian/MeshUtils.py b/hera/simulations/gaussian/MeshUtils.py index b6cdfbd07..db867a78e 100644 --- a/hera/simulations/gaussian/MeshUtils.py +++ b/hera/simulations/gaussian/MeshUtils.py @@ -2,7 +2,7 @@ import xarray import numpy -from hera.utils.unitHandler import * +from hera.utils.unitHandler import ureg, unumToPint from scipy.special import erf class GaussianToMesh(object): @@ -95,7 +95,7 @@ def _defineCoordinates(self,gaussians): return xCoord,yCoord - def gaussianToMesh(self,gaussians,groupby=None,QName="Q",QUnits=kg): + def gaussianToMesh(self,gaussians,groupby=None,QName="Q",QUnits=ureg.kg): """ Spreads the list of gaussians on the regular mesh. if groupby is None then sum all the gaussians, @@ -125,7 +125,7 @@ def gaussianToMesh(self,gaussians,groupby=None,QName="Q",QUnits=kg): return ret - def _gaussianToMesh(self,gaussian,xCoords,yCoords,addTo=None,QName="Q",QUnits=kg): + def _gaussianToMesh(self,gaussian,xCoords,yCoords,addTo=None,QName="Q",QUnits=ureg.kg): """ Return a single gaussian that was spread on the mesh. @@ -154,7 +154,7 @@ def _gaussianToMesh(self,gaussian,xCoords,yCoords,addTo=None,QName="Q",QUnits=kg gaussY = xarray.DataArray(1/(twoPiFactor*sigy)*numpy.exp(-0.5*((yCoords-y)/sigy)**2), dims='y',coords={'y':yCoords}) fullX, fullY = xarray.broadcast(gaussX, gaussY) - ret = fullX*fullY*gaussian[QName] *((1*kg).asNumber(QUnits)) + ret = fullX*fullY*gaussian[QName] * (1*ureg.kg).m_as(unumToPint(QUnits)) if addTo is not None: ret = addTo+ ret @@ -166,7 +166,7 @@ class GaussianIntegrationToMesh(GaussianToMesh): This class calculates the concentration of meterial there is in the mesh by subtracting erf. This is more accurate than estimating the gaussian. """ - def _gaussianToMesh(self,gaussian,xCoords,yCoords,addTo=None,QName="Q",QUnits=kg): + def _gaussianToMesh(self,gaussian,xCoords,yCoords,addTo=None,QName="Q",QUnits=ureg.kg): """ Return the concentration a single gaussian that was spread on the mesh. @@ -198,7 +198,7 @@ def _gaussianToMesh(self,gaussian,xCoords,yCoords,addTo=None,QName="Q",QUnits=kg fullX, fullY = xarray.broadcast(gaussX, gaussY) - ret = fullX*fullY*gaussian[QName].asNumber(QUnits) + ret = fullX*fullY*unumToPint(gaussian[QName]).m_as(unumToPint(QUnits)) if addTo is not None: ret = addTo+ ret return ret diff --git a/hera/simulations/gaussian/Meteorology.py b/hera/simulations/gaussian/Meteorology.py index 55100a35a..731909acb 100644 --- a/hera/simulations/gaussian/Meteorology.py +++ b/hera/simulations/gaussian/Meteorology.py @@ -1,7 +1,7 @@ import pandas import numpy -from unum.units import * from hera.utils import * +from hera.utils.unitHandler import ureg, unumToPint # from hera.utils import tounit, tonumber @@ -31,7 +31,7 @@ def ustar(self): @ustar.setter def ustar(self, value): - self._ustar = tounit(value, m / s) + self._ustar = tounit(value, ureg.m / ureg.s) @property def refHeight(self): @@ -39,7 +39,7 @@ def refHeight(self): @refHeight.setter def refHeight(self, value): - self._refHeight = tounit(value, m) + self._refHeight = tounit(value, ureg.m) @property def wind_p(self): @@ -51,7 +51,7 @@ def u10(self): @u10.setter def u10(self, value): - self._u_refHeight = tounit(value, m / s) + self._u_refHeight = tounit(value, ureg.m / ureg.s) self._refHeight = 10 @property @@ -60,7 +60,7 @@ def u_refHeight(self): @u_refHeight.setter def u_refHeight(self, value): - self._u_refHeight = tounit(value, m / s) + self._u_refHeight = tounit(value, ureg.m / ureg.s) @property def stability(self): @@ -83,7 +83,7 @@ def z0(self): @z0.setter def z0(self, value): - self._z0 = tounit(value, m) + self._z0 = tounit(value, ureg.m) self._setPvalues() @property @@ -92,7 +92,7 @@ def temperature(self): @temperature.setter def temperature(self, value): - self._temperature = tounit(value, K) + self._temperature = tounit(value, ureg.K) @property def skinSurfaceTemperature(self): @@ -100,7 +100,7 @@ def skinSurfaceTemperature(self): @skinSurfaceTemperature.setter def skinSurfaceTemperature(self, value): - self._skinSurfaceTemperature = tounit(value, K) + self._skinSurfaceTemperature = tounit(value, ureg.K) # ================================ Wind profile calcluation ## Calculate the wind profile based on stability and z0. @@ -168,9 +168,9 @@ def getAirPressure(self, height): :return: The air pressure at mmHg units. """ - height = tounit(height, m) + height = unumToPint(height).to(ureg.m) - return 760. * numpy.exp(-1.186e-4 * height.asNumber(m)) * mmHg + return 760. * numpy.exp(-1.186e-4 * height.m_as(ureg.m)) * ureg.mmHg def getTKE(self, height): """ @@ -195,7 +195,7 @@ def getAirTemperature(self, height): :return: air temperature at C. """ - return self._temperature - 6.5e-3 * tonumber(height, m) * celsius + return self._temperature - 6.5e-3 * tonumber(height, ureg.m) * ureg.K def getAirDensity(self, height): """ @@ -214,11 +214,11 @@ def getAirDensity(self, height): :return: air density in kg/m**3 """ - P = self.getAirPressure(height).asNumber(mmHg) - T = self.getAirTemperature(height).asNumber(celsius) + P = unumToPint(self.getAirPressure(height)).m_as(ureg.mmHg) + T = unumToPint(self.getAirTemperature(height)).m_as(ureg.degC) - density = 1.701316e-6 * P / (1 + 0.00367 * T) * g / cm ** 3 - return density.asUnit(kg / m ** 3) + density = 1.701316e-6 * P / (1 + 0.00367 * T) * ureg.g / ureg.cm ** 3 + return density.to(ureg.kg / ureg.m ** 3) def getAirDynamicViscosity(self, height): """ @@ -233,8 +233,8 @@ def getAirDynamicViscosity(self, height): :return: The air viscosity in dyne/sec/m**2. """ - T = self.getAirTemperature(height).asNumber(celsius) - return (1e-6 * (170.27 + 0.911409 * T - 0.00786742 * T ** 2) * dyne * s / cm ** 2).asUnit(dyne * s / m ** 2) + T = unumToPint(self.getAirTemperature(height)).m_as(ureg.degC) + return (1e-6 * (170.27 + 0.911409 * T - 0.00786742 * T ** 2) * ureg.dyne * ureg.s / ureg.cm ** 2).to(ureg.dyne * ureg.s / ureg.m ** 2) def _setPvalues(self): """ @@ -247,7 +247,7 @@ def _setPvalues(self): if (self.z0 is None or self.stability is None): return pstab = self._pvalues[self.stability] - self._wind_p = numpy.interp(self.z0.asNumber(m), pstab.index, pstab) + self._wind_p = numpy.interp(unumToPint(self.z0).m_as(ureg.m), pstab.index, pstab) def getWindVelocity(self, height): """ @@ -263,8 +263,8 @@ def getWindVelocity(self, height): :return: The wind velocity at the requested height. """ - height = tonumber(height, m) - refHeight = tonumber(self.refHeight, m) + height = tonumber(height, ureg.m) + refHeight = tonumber(self.refHeight, ureg.m) height = numpy.min([numpy.max([height, 0]), 300]) return self.u_refHeight * (height / refHeight) ** self.wind_p @@ -287,14 +287,14 @@ def getWindVelocity(self, height): The wind velocity at the requested height. """ - z0 = tonumber(self.z0, m) - height = tonumber(height, m) + z0 = tonumber(self.z0, ureg.m) + height = tonumber(height, ureg.m) height = numpy.min([numpy.max([height, 0]), 300]) - refHeight = tonumber(self.refHeight, m) + refHeight = tonumber(self.refHeight, ureg.m) ustar_over_kappa = self.u_refHeight / numpy.log(refHeight / z0) if (height <= z0): - u = 0 * m / s + u = 0 * ureg.m / ureg.s else: u = ustar_over_kappa * numpy.log(height / z0) @@ -333,7 +333,7 @@ def __init__(self): log=StandardMeteorolgyConstant_log, uniformWind=StandardMeteorolgyConstant_uniformWind) - def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=20*celsius, stability="D", z0=0.1*m, ustar=0.3*m/s, skinSurfaceTemperature=35*celsius): + def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=20*ureg.degC, stability="D", z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC): """ Creating a meteorology object. @@ -348,8 +348,8 @@ def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", tempe return self.meteorology[verticalProfileType](u10=u10, inversion=inversion, temperature=temperature, stability=stability, z0=z0, ustar=ustar, skinSurfaceTemperature=skinSurfaceTemperature) - def getMeteorologyFromURefHeight(self, u, inversion, refHeight, verticalProfileType="log", temperature=20*celsius, stability="D", - z0=0.1*m, ustar=0.3*m/s, skinSurfaceTemperature=35*celsius): + def getMeteorologyFromURefHeight(self, u, inversion, refHeight, verticalProfileType="log", temperature=20*ureg.degC, stability="D", + z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC): """ Creating a meteorology object. @@ -365,4 +365,3 @@ def getMeteorologyFromURefHeight(self, u, inversion, refHeight, verticalProfileT stability=stability, z0=z0, ustar=ustar, skinSurfaceTemperature=skinSurfaceTemperature) - diff --git a/hera/simulations/gaussian/Sigma.py b/hera/simulations/gaussian/Sigma.py index 77389e42f..dceda27d0 100644 --- a/hera/simulations/gaussian/Sigma.py +++ b/hera/simulations/gaussian/Sigma.py @@ -1,7 +1,7 @@ import pandas import numpy from hera.utils import tounit,tonumber -from unum.units import * +from hera.utils.unitHandler import ureg, unumToPint # from sympy import solve, symbols from scipy import optimize @@ -32,17 +32,17 @@ def getVirtualDistance(self,sigma0,stability): stability=stability sx,sy,sz = sigma0 - sx = tonumber(sx, m) - sy = tonumber(sy, m) - sz = tonumber(sz, m) + sx = tonumber(sx, ureg.m) + sy = tonumber(sy, ureg.m) + sz = tonumber(sz, ureg.m) sigma_x = lambda x: self.getSigma(x=x, stability=stability, units=False)['sigmaX'][0] - sx sigma_y = lambda x: self.getSigma(x=x, stability=stability, units=False)['sigmaY'][0] - sy sigma_z = lambda x: self.getSigma(x=x, stability=stability, units=False)['sigmaZ'][0] - sz - Ix = optimize.newton(sigma_x, x0=0)*m - Iy = optimize.newton(sigma_y, x0=0) * m - Iz = optimize.newton(sigma_z, x0=0) * m + Ix = optimize.newton(sigma_x, x0=0)*ureg.m + Iy = optimize.newton(sigma_y, x0=0) * ureg.m + Iz = optimize.newton(sigma_z, x0=0) * ureg.m return Ix, Iy, Iz @@ -111,7 +111,7 @@ def getSigma(self, x, stability, sigma0=None, units=True): """ - x = numpy.array([tonumber(y,m) for y in numpy.atleast_1d(x)]) + x = numpy.array([tonumber(y,ureg.m) for y in numpy.atleast_1d(x)]) Ax, Bx, Cx = self._coeffX.loc[stability][['A','B','C']] Az, Bz, Cz = self._coeffZ.loc[stability][['A', 'B', 'C']] @@ -122,9 +122,9 @@ def getSigma(self, x, stability, sigma0=None, units=True): Ix,Iy,Iz = self.getVirtualDistance(sigma0,stability) - Ix = tonumber(Ix,m) - Iy = tonumber(Iy,m) - Iz = tonumber(Iz,m) + Ix = tonumber(Ix,ureg.m) + Iy = tonumber(Iy,ureg.m) + Iz = tonumber(Iz,ureg.m) # pandas_res = pandas.DataFrame({ # 'sigmaX' : [Ax*(x+Ix)*(1+Bx*(x+Ix))**Cx *m], @@ -134,9 +134,9 @@ def getSigma(self, x, stability, sigma0=None, units=True): if units: dict_res = { - 'sigmaX': Ax * (x + Ix) * (1 + Bx * (x + Ix)) ** Cx * m, - 'sigmaY': Ax * (x + Iy) * (1 + Bx * (x + Iy)) ** Cx * m, - 'sigmaZ': Az * (x + Iz) * (1 + Bz * (x + Iz)) ** Cz * m, 'distance': x * m + 'sigmaX': Ax * (x + Ix) * (1 + Bx * (x + Ix)) ** Cx * ureg.m, + 'sigmaY': Ax * (x + Iy) * (1 + Bx * (x + Iy)) ** Cx * ureg.m, + 'sigmaZ': Az * (x + Iz) * (1 + Bz * (x + Iz)) ** Cz * ureg.m, 'distance': x * ureg.m } else: diff --git a/hera/simulations/gaussian/gasCloud.py b/hera/simulations/gaussian/gasCloud.py index 9a9efee8b..51d9699b0 100644 --- a/hera/simulations/gaussian/gasCloud.py +++ b/hera/simulations/gaussian/gasCloud.py @@ -2,8 +2,8 @@ import numpy from numpy import matlib import math -from unum.units import * from hera.utils import * +from hera.utils.unitHandler import ureg, unumToPint import xarray from scipy import special @@ -16,16 +16,16 @@ def __init__(self, sourceQ, sourceHeight, initialCloudSize, sigmaType): Parameters ---------- - sourceQ : unum, method - If unum: + sourceQ : pint Quantity or unum + If Quantity: The unit determine the release time. [mass] - Instantaneous [mass/time] - Continuous else Continuous (not implementaed yet.) - sourceHeight : unum - initialCloudSize : 3-touple unum, the sigmas in each axis. + sourceHeight : pint Quantity + initialCloudSize : 3-touple pint Quantity, the sigmas in each axis. sigmaType : The sigma type, for example from Briggs, rural/urban. """ self.sourceHeight = sourceHeight @@ -41,29 +41,30 @@ def createGasCloud(sourceQ,sourceHeight,initialCloudSize,sigmaType): Return the type of the release based on the units of Q Parameters ---------- - sourceQ : unum, method - If unum: + sourceQ : pint Quantity or unum + If Quantity: The unit determine the release time. [mass] - Instantaneous [mass/time] - Continuous else Continuous (not implementaed yet.) - sourceHeight : unum - initialCloudSize : 3-touple unum, the sigmas in each axis. + sourceHeight : pint Quantity + initialCloudSize : 3-touple pint Quantity, the sigmas in each axis. Returns ------- """ + sourceQ = unumToPint(sourceQ) try: - sourceQ.asUnit(mg) + sourceQ.to(ureg.mg) instantaneous = True - except: + except Exception: try: - sourceQ.asUnit(mg/min) + sourceQ.to(ureg.mg/ureg.min) instantaneous = False - except: + except Exception: raise ValueError("Must be mass or mass per time!") returnCls = instantaneousReleaseGasCloud if instantaneous else continuousReleaseGasCloud @@ -75,10 +76,10 @@ def _getTXterm(self, stability, u, xcoordRange, tcoordRange): """ Parameters ---------- - initialCloudSize : 3-tuple of float/unum (default m) + initialCloudSize : 3-tuple of float/pint Quantity (default m) The initial cloud size (standard deviation) in the x,y and z dimensions. stability : the stability class - u: the wind speed /unum (default m/s) + u: the wind speed / pint Quantity (default m/s) xcoordRange : Tuple in numpy.arange format, unitless. tcoordRange : Tuple in numpy.arange format, unitless. @@ -99,7 +100,7 @@ def _getXYterm(self, stability, xcoordRange, ycoordRange): Parameters ---------- - initialCloudSize : 3-tuple of float/unum (default m) + initialCloudSize : 3-tuple of float/pint Quantity (default m) The initial cloud size (standard deviation) in the x,y and z dimensions. stability : the stability class xcoordRange : Tuple in numpy.arange format, unitless. @@ -123,10 +124,10 @@ def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflec Parameters ---------- - initialCloudSize : 3-tuple of float/unum (default m) + initialCloudSize : 3-tuple of float/pint Quantity (default m) The initial cloud size (standard deviation) in the x,y and z dimensions. stability : the stability class - inversion: the inversion height /unum (default m) + inversion: the inversion height / pint Quantity (default m) xcoordRange : Tuple in numpy.arange format, unitless. zcoordRange : Tuple in numpy.arange format, unitless. numOfReflections : The number of reflections of the summation of the Z component. @@ -135,7 +136,7 @@ def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflec ------- The Z component of the Gaussian concentration formula. """ - sourceHeight = tonumber(self.sourceHeight, m) + sourceHeight = tonumber(self.sourceHeight, ureg.m) X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] @@ -152,7 +153,7 @@ def _getXZterm(self, stability, inversion, xcoordRange, zcoordRange, numOfReflec - def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min): + def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min): """ This function generates an xarray of fractions of the mass of each isotope at every time-step. @@ -164,10 +165,10 @@ def fractions(self, fracVector, minx, miny, minz, maxx, maxy, maxz, timeSpan, dx :return: """ - xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy,m)) - ycoordRange = numpy.arange(tonumber(miny,m),tonumber(maxy,m),tonumber(dxdy,m)) - zcoordRange = numpy.arange(tonumber(minz, m), tonumber(maxz, m), tonumber(dz,m)) - tcoordRange = numpy.arange(0,tonumber(timeSpan,min),tonumber(dt,min)) + xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy,ureg.m)) + ycoordRange = numpy.arange(tonumber(miny,ureg.m),tonumber(maxy,ureg.m),tonumber(dxdy,ureg.m)) + zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz,ureg.m)) + tcoordRange = numpy.arange(0,tonumber(timeSpan,ureg.min),tonumber(dt,ureg.min)) frac, X = numpy.meshgrid(fracVector, xcoordRange, indexing='ij') XR_downwind = xarray.DataArray(frac, dims=("time", "x"), coords={"time": tcoordRange, "x": xcoordRange}) @@ -183,10 +184,10 @@ def _getTXDosage(self, stability, u, xcoordRange, tcoordRange): """ Parameters ---------- - initialCloudSize : 3-tuple of float/unum (default m) + initialCloudSize : 3-tuple of float/pint Quantity (default m) The initial cloud size (standard deviation) in the x,y and z dimensions. stability : the stability class - u: the wind speed /unum (default m/s) + u: the wind speed / pint Quantity (default m/s) xcoordRange : Tuple in numpy.arange format, unitless. tcoordRange : Tuple in numpy.arange format, unitless. @@ -290,11 +291,11 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): :return: The Depletion Factor """ - v = tonumber(0.003*m/s, m/min) #Deposition velocity. By default we take this value to be 0.003 [m/s] + v = tonumber(0.003*ureg.m/ureg.s, ureg.m/ureg.min) #Deposition velocity. By default we take this value to be 0.003 [m/s] X, Z = numpy.meshgrid(xcoordRange, zcoordRange, indexing='ij') sigmaZ = self.sigmaType.getSigma(x=X, stability=stability, sigma0=self.initialCloudSize, units=False)['sigmaZ'] - H = tonumber(self.sourceHeight, m) + H = tonumber(self.sourceHeight, ureg.m) dx = xcoordRange[1] - xcoordRange[0] DF = (numpy.e**(numpy.cumsum(1/(sigmaZ*numpy.e**(0.5*(H/sigmaZ)**2))*dx, axis=0)))**((-v/u)*numpy.sqrt(2/numpy.pi)) @@ -302,14 +303,14 @@ def _getDF(self, stability, u, xcoordRange, zcoordRange): - def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*m, dz=1*m, dt=1*min): + def getDF_noQ_xarray(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min): stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) + u = tonumber(meteorology.u10, ureg.m/ureg.min) - xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy, m)) - ycoordRange = numpy.arange(tonumber(miny, m), tonumber(maxy, m), tonumber(dxdy, m)) - zcoordRange = numpy.arange(tonumber(minz, m), tonumber(maxz, m), tonumber(dz, m)) - tcoordRange = numpy.arange(0, tonumber(timeSpan, min), tonumber(dt, min)) + xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy, ureg.m)) + ycoordRange = numpy.arange(tonumber(miny, ureg.m), tonumber(maxy, ureg.m), tonumber(dxdy, ureg.m)) + zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz, ureg.m)) + tcoordRange = numpy.arange(0, tonumber(timeSpan, ureg.min), tonumber(dt, ureg.min)) TX = self._getTXterm_ones(xcoordRange=xcoordRange, tcoordRange=tcoordRange) XY = self._getXYterm_ones(xcoordRange=xcoordRange, ycoordRange=ycoordRange) @@ -325,16 +326,16 @@ class instantaneousReleaseGasCloud(abstractGasCloud): def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3): - xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy,m)) - ycoordRange = numpy.arange(tonumber(miny,m),tonumber(maxy,m),tonumber(dxdy,m)) - zcoordRange = numpy.arange(tonumber(minz, m), tonumber(maxz, m), tonumber(dz,m)) - tcoordRange = numpy.arange(0,tonumber(timeSpan,min),tonumber(dt,min)) + xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy,ureg.m)) + ycoordRange = numpy.arange(tonumber(miny,ureg.m),tonumber(maxy,ureg.m),tonumber(dxdy,ureg.m)) + zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz,ureg.m)) + tcoordRange = numpy.arange(0,tonumber(timeSpan,ureg.min),tonumber(dt,ureg.min)) stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) - inversion = tonumber(meteorology.inversion, m) + u = tonumber(meteorology.u10, ureg.m/ureg.min) + inversion = tonumber(meteorology.inversion, ureg.m) TX = self._getTXterm(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) XY = self._getXYterm(stability=stability, xcoordRange=xcoordRange, ycoordRange=ycoordRange) @@ -342,21 +343,21 @@ def getConcentrationFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz numOfReflections=numOfReflections) ret = TX*XY*XZ - ret.attrs['Q'] = 1/m**3 + ret.attrs['Q'] = 1/ureg.m**3 return ret def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): stability = meteorology.stability - u = tonumber(meteorology.u10, m/min) - inversion = tonumber(meteorology.inversion, m) + u = tonumber(meteorology.u10, ureg.m/ureg.min) + inversion = tonumber(meteorology.inversion, ureg.m) - xcoordRange = numpy.arange(tonumber(minx, m), tonumber(maxx, m), tonumber(dxdy, m)) - ycoordRange = numpy.arange(tonumber(miny, m), tonumber(maxy, m), tonumber(dxdy, m)) - zcoordRange = numpy.arange(tonumber(minz, m), tonumber(maxz, m), tonumber(dz, m)) - tcoordRange = numpy.arange(0, tonumber(timeSpan, min), tonumber(dt, min)) + xcoordRange = numpy.arange(tonumber(minx, ureg.m), tonumber(maxx, ureg.m), tonumber(dxdy, ureg.m)) + ycoordRange = numpy.arange(tonumber(miny, ureg.m), tonumber(maxy, ureg.m), tonumber(dxdy, ureg.m)) + zcoordRange = numpy.arange(tonumber(minz, ureg.m), tonumber(maxz, ureg.m), tonumber(dz, ureg.m)) + tcoordRange = numpy.arange(0, tonumber(timeSpan, ureg.min), tonumber(dt, ureg.min)) TX = self._getTXDosage(stability=stability, u=u, xcoordRange=xcoordRange, tcoordRange=tcoordRange) XY = self._getXYterm(stability=meteorology.stability, xcoordRange=xcoordRange, ycoordRange=ycoordRange) @@ -367,12 +368,12 @@ def getDosageFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, D_F = self._getDF(stability=stability, u=u, xcoordRange=xcoordRange, zcoordRange=zcoordRange) ret = TX*XY*(XZ*D_F) - ret.attrs['Q'] = 1*min/m**3 + ret.attrs['Q'] = 1*ureg.min/ureg.m**3 return ret def getDosageFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, dz=dz, dt=dt,numOfReflections=numOfReflections) @@ -384,7 +385,7 @@ def getDosageFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt) ret = D_without_Q*D_F - ret.attrs['Q'] = 1*min/m**3 + ret.attrs['Q'] = 1*ureg.min/ureg.m**3 return ret @@ -392,37 +393,37 @@ def getDosageFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, def getConcentrationFromMinMaxRange_inst(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): C_without_Q = self.getConcentrationFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, dz=dz, dt=dt,numOfReflections=numOfReflections) - ret = tonumber(self.sourceQ, mg)*C_without_Q - ret.attrs['Q'] = 1*mg/m**3 + ret = tonumber(self.sourceQ, ureg.mg)*C_without_Q + ret.attrs['Q'] = 1*ureg.mg/ureg.m**3 return ret def getDosageFromMinMaxRange_inst(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - ret = tonumber(self.sourceQ, mg)*D_without_Q - ret.attrs['Q'] = 1*mg*min/m**3 + ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q + ret.attrs['Q'] = 1*ureg.mg*ureg.min/ureg.m**3 return ret def getDosageFromMinMaxRange_inst_NoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy,maxz=maxz, timeSpan=timeSpan,dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - ret = tonumber(self.sourceQ, mg)*D_without_Q - ret.attrs['Q'] = 1*mg*min/m**3 + ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q + ret.attrs['Q'] = 1*ureg.mg*ureg.min/ureg.m**3 return ret @@ -430,72 +431,78 @@ def getDosageFromMinMaxRange_inst_NoERF(self, meteorology, minx, miny, minz, max #---------------------------------------------Radiology--------------------------------------------- def concentrationConversion_mass_to_Bq(self, C, outputUnits, specificActivity): - factor = (C.attrs['Q'] * specificActivity).asNumber(outputUnits) + units = unumToPint(C.attrs['Q']) + out_units = unumToPint(outputUnits) + factor = (units * specificActivity).m_as(out_units) C_Bq = C * factor - C_Bq.attrs['Q'] = outputUnits + C_Bq.attrs['Q'] = out_units return C_Bq def getTIACFromMinMaxRange_inst(self,specifitActivity, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, outputUnits=Bq*s/m**3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False): D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - ret = tonumber(self.sourceQ, mg)*D_without_Q - factor = (1*mg*min/m**3 * specifitActivity).asNumber(outputUnits) + out_units = unumToPint(outputUnits) + ret = tonumber(self.sourceQ, ureg.mg)*D_without_Q + factor = (1*ureg.mg*ureg.min/ureg.m**3 * specifitActivity).m_as(out_units) ret *= factor - ret.attrs['Q'] = outputUnits + ret.attrs['Q'] = out_units return ret def getTIACFromMinMaxRange_inst_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, outputUnits=s/m**3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False): D_without_Q = self.getDosageFromMinMaxRange_inst_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - currentUnites = D_without_Q.attrs['Q'] - factor = (currentUnites).asNumber(outputUnits) + out_units = unumToPint(outputUnits) + currentUnites = unumToPint(D_without_Q.attrs['Q']) + factor = currentUnites.m_as(out_units) D_without_Q *= factor - D_without_Q.attrs['Q'] = outputUnits + D_without_Q.attrs['Q'] = out_units return D_without_Q def getTIACFromMinMaxRange_inst_NoERF(self, specifitActivity, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, outputUnits=Bq*s/m**3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False): D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - ret = tonumber(self.sourceQ, mg) * D_without_Q - factor = (1 * mg * min / m ** 3 * specifitActivity).asNumber(outputUnits) + out_units = unumToPint(outputUnits) + ret = tonumber(self.sourceQ, ureg.mg) * D_without_Q + factor = (1 * ureg.mg * ureg.min / ureg.m ** 3 * specifitActivity).m_as(out_units) ret *= factor - ret.attrs['Q'] = outputUnits + ret.attrs['Q'] = out_units return ret def getTIACFromMinMaxRange_inst_NoERF_noQ(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, outputUnits=s/m**3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, outputUnits=ureg.s/ureg.m**3, DF=False): D_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - currentUnites = D_without_Q.attrs['Q'] - factor = (currentUnites).asNumber(outputUnits) + out_units = unumToPint(outputUnits) + currentUnites = unumToPint(D_without_Q.attrs['Q']) + factor = currentUnites.m_as(out_units) D_without_Q *= factor - D_without_Q .attrs['Q'] = outputUnits + D_without_Q.attrs['Q'] = out_units return D_without_Q def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, outputUnits=Bq*s/m**3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, outputUnits=ureg.Bq*ureg.s/ureg.m**3, DF=False): """ :param C: xarray of concentrations in units of [mass/volume]. @@ -504,11 +511,12 @@ def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, meteorology, :return: TIAC (Time Integrated Air Concentration) in unites of [Bq*time/volume] """ - factor = (C.attrs['Q']).asNumber(mg/m**3) + out_units = unumToPint(outputUnits) + factor = unumToPint(C.attrs['Q']).m_as(ureg.mg/ureg.m**3) C_mg_m3 = C*factor - C_mg_m3.attrs['Q'] = mg/m**3 + C_mg_m3.attrs['Q'] = ureg.mg/ureg.m**3 - C_without_Q = C_mg_m3 / tonumber(self.sourceQ, mg) + C_without_Q = C_mg_m3 / tonumber(self.sourceQ, ureg.mg) D_without_Q = self.trapezoidal_integration(data = C_without_Q) D_F = 1 @@ -518,15 +526,15 @@ def getTIACFromConcentration_inst_NoERF(self, C, specifitActivity, meteorology, D_without_Q = D_without_Q * D_F - ret = tonumber(self.sourceQ, mg) * D_without_Q - factor = (1 * mg * min / m ** 3 * specifitActivity).asNumber(outputUnits) + ret = tonumber(self.sourceQ, ureg.mg) * D_without_Q + factor = (1 * ureg.mg * ureg.min / ureg.m ** 3 * specifitActivity).m_as(out_units) ret *= factor - ret.attrs['Q'] = outputUnits + ret.attrs['Q'] = out_units return ret def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, outputUnits=s/m**3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, outputUnits=ureg.s/ureg.m**3, DF=False): """ :param C_noQ: xarray of concentrations in units of [1/volume]. @@ -541,12 +549,13 @@ def getTIACFromConcentration_inst_NoERF_noQ(self, C_noQ, meteorology, minx, miny maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt) D_without_Q = D_without_Q*D_F - D_without_Q.attrs['Q'] = 1*min/m**3 #need to verify that C_noQ was generated with time steps in [min], not [s] + D_without_Q.attrs['Q'] = 1*ureg.min/ureg.m**3 #need to verify that C_noQ was generated with time steps in [min], not [s] - currentUnites = D_without_Q.attrs['Q'] - factor = (currentUnites).asNumber(outputUnits) + out_units = unumToPint(outputUnits) + currentUnites = unumToPint(D_without_Q.attrs['Q']) + factor = currentUnites.m_as(out_units) D_without_Q *= factor - D_without_Q .attrs['Q'] = outputUnits + D_without_Q.attrs['Q'] = out_units return D_without_Q @@ -578,7 +587,7 @@ def get_TIAC_from_dist(self, data, y, z, dist_list): class continuousReleaseGasCloud(abstractGasCloud): def getConcentrationFromMinMaxRange_cont(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): """ Returns ------- @@ -590,11 +599,11 @@ def getConcentrationFromMinMaxRange_cont(self, meteorology, minx, miny, minz, ma maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - return tonumber(self.sourceQ, mg/s)*C_without_Q + return tonumber(self.sourceQ, ureg.mg/ureg.s)*C_without_Q def getConcentrationFromMinMaxRange_cont_NoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): """ Returns ------- @@ -605,28 +614,28 @@ def getConcentrationFromMinMaxRange_cont_NoERF(self, meteorology, minx, miny, mi C_without_Q = self.getDosageFromMinMaxRange_inst_NoERF_noQ(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) - return tonumber(self.sourceQ, mg/s)*C_without_Q + return tonumber(self.sourceQ, ureg.mg/ureg.s)*C_without_Q def getDosageFromMinMaxRange_cont_NoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): C_without_Q = self.getConcentrationFromMinMaxRange_cont(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) D_without_Q = self.trapezoidal_integration(data=C_without_Q) - return tonumber(self.sourceQ, mg/min) * D_without_Q + return tonumber(self.sourceQ, ureg.mg/ureg.min) * D_without_Q def getDosageFromMinMaxRange_cont_doubleNoERF(self, meteorology, minx, miny, minz, maxx, maxy, maxz, timeSpan, - dxdy=10*m, dz=1*m, dt=1*min, numOfReflections=3, DF=False): + dxdy=10*ureg.m, dz=1*ureg.m, dt=1*ureg.min, numOfReflections=3, DF=False): C_without_Q = self.getConcentrationFromMinMaxRange_cont_NoERF(meteorology=meteorology, minx=minx, miny=miny, minz=minz, maxx=maxx, maxy=maxy, maxz=maxz, timeSpan=timeSpan, dxdy=dxdy, dz=dz, dt=dt, numOfReflections=numOfReflections, DF=DF) D_without_Q = self.trapezoidal_integration(data=C_without_Q) - return tonumber(self.sourceQ, mg/min) * D_without_Q + return tonumber(self.sourceQ, ureg.mg/ureg.min) * D_without_Q @@ -638,7 +647,7 @@ class Continuous(object): Timekernel = None _FullKernel = None - def __init__(self,dt,kernelsize,timetofinish=10*min): + def __init__(self,dt,kernelsize,timetofinish=10*ureg.min): """ Time to finish. the time (min) it take to reach 0.1. @@ -655,19 +664,19 @@ def __init__(self,dt,kernelsize,timetofinish=10*min): and that is what we should put in the kernel. """ - dt = tounum(dt,min) - self.dt = tonumber(dt,min) + dt = unumToPint(dt).to(ureg.min) + self.dt = dt.m_as(ureg.min) self.kernelsize = kernelsize - timetofinish = tounum(timetofinish,min) - alpha = numpy.log(0.1)/(-timetofinish.asNumber(min)) + timetofinish = unumToPint(timetofinish).to(ureg.min) + alpha = numpy.log(0.1)/(-timetofinish.m_as(ureg.min)) # build the kernel. - ts = numpy.arange(kernelsize,-1,-1)*dt.asNumber(min) + ts = numpy.arange(kernelsize,-1,-1)*dt.m_as(ureg.min) # the kernel - self.Timekernel = (numpy.exp(-alpha*ts[1:]) - numpy.exp(-alpha*ts[:-1]))/dt.asNumber(min) + self.Timekernel = (numpy.exp(-alpha*ts[1:]) - numpy.exp(-alpha*ts[:-1]))/dt.m_as(ureg.min) self.Timekernel = self.Timekernel.reshape([kernelsize,1,1,1]) @@ -685,4 +694,3 @@ def calc(self,data): - diff --git a/hera/simulations/gaussian/toolkit.py b/hera/simulations/gaussian/toolkit.py index 50c41787d..bdfb6e5a5 100644 --- a/hera/simulations/gaussian/toolkit.py +++ b/hera/simulations/gaussian/toolkit.py @@ -2,7 +2,7 @@ from hera.simulations.gaussian.Sigma import BriggsRural from hera.simulations.gaussian.gasCloud import abstractGasCloud from hera.simulations.gaussian.Meteorology import MeteorologyFactory -from unum.units import * +from hera.utils.unitHandler import ureg, unumToPint, celsius, K import matplotlib.pyplot as plt @@ -58,14 +58,14 @@ def listSigmaTypes(self): - def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=20*celsius, stability="D", - z0=0.1*m, ustar=0.3*m/s, skinSurfaceTemperature=35*celsius): + def getMeteorologyFromU10(self, u10, inversion, verticalProfileType="log", temperature=20*ureg.degC, stability="D", + z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC): return MeteorologyFactory().getMeteorologyFromU10(u10=u10, inversion=inversion, verticalProfileType=verticalProfileType, temperature=temperature, stability=stability, z0=z0, ustar=ustar, skinSurfaceTemperature=skinSurfaceTemperature) - def getMeteorologyFromURefHeight(self, u, refHeight, inversion, verticalProfileType="log", temperature=20*celsius, stability="D", - z0=0.1*m, ustar=0.3*m/s, skinSurfaceTemperature=35*celsius): + def getMeteorologyFromURefHeight(self, u, refHeight, inversion, verticalProfileType="log", temperature=20*ureg.degC, stability="D", + z0=0.1*ureg.m, ustar=0.3*ureg.m/ureg.s, skinSurfaceTemperature=35*ureg.degC): return MeteorologyFactory().getMeteorologyFromURefHeight(u=u, refHeight=refHeight, inversion=inversion, verticalProfileType=verticalProfileType, temperature=temperature, stability=stability, z0=z0, ustar=ustar,skinSurfaceTemperature=skinSurfaceTemperature) @@ -124,8 +124,8 @@ def plotMaxConcentrationOverTime(self,C, y, z, x_min=None,x_max=None): :param z: Z-coordinate of the line :return: """ - y = y.asNumber(m) - z = z.asNumber(m) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) conc_x_inst = C.sel(y=y, z=z, method="nearest").max(dim="time") x_array = conc_x_inst.squeeze().x @@ -139,8 +139,8 @@ def plotMaxConcentrationOverTime(self,C, y, z, x_min=None,x_max=None): plt.title(f"Maximum concentration over time. y={y}[m], z={z}[m]") plt.grid() if x_min is not None and x_max is not None: - x_min = x_min.asNumber(m) - x_max = x_max.asNumber(m) + x_min = unumToPint(x_min).m_as(ureg.m) + x_max = unumToPint(x_max).m_as(ureg.m) plt.xlim(x_min, x_max) plt.show() @@ -155,9 +155,9 @@ def plotFixedPointConcentrationOverTime(self, C, x, y, z, t_min=None,t_max=None) :return: """ - x = x.asNumber(m) - y = y.asNumber(m) - z = z.asNumber(m) + x = unumToPint(x).m_as(ureg.m) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) conc_inst_t = C.sel(x=x, y=y, z=z, method="nearest") time_array = conc_inst_t.squeeze().time @@ -171,8 +171,8 @@ def plotFixedPointConcentrationOverTime(self, C, x, y, z, t_min=None,t_max=None) plt.title(f"Receptor at x={x}[m], y={y}[m], z={z}[m].") plt.grid() if t_min is not None and t_max is not None: - t_min = t_min.asNumber(min) - t_max = t_max.asNumber(min) + t_min = unumToPint(t_min).m_as(ureg.min) + t_max = unumToPint(t_max).m_as(ureg.min) plt.xlim(t_min, t_max) plt.show() @@ -187,9 +187,9 @@ def plotDosagePerDistance(self,D, y, z, time): :param z: Z-coordinate of the line :return: """ - y = y.asNumber(m) - z = z.asNumber(m) - time = time.asNumber(min) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) + time = unumToPint(time).m_as(ureg.min) dos_x_inst = D.sel(y=y, z=z, time=time, method="nearest") x_array = dos_x_inst .squeeze().x @@ -215,9 +215,9 @@ def plotFixedPointDosageOverTime(self, D, x, y, z): :return: """ - x = x.asNumber(m) - y = y.asNumber(m) - z = z.asNumber(m) + x = unumToPint(x).m_as(ureg.m) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) dos_inst_t = D.sel(x=x, y=y, z=z, method="nearest") time_array = dos_inst_t .squeeze().time @@ -245,8 +245,8 @@ def plotMaxConcentrationOverTime_noQ(self,C, y, z, x_min=None,x_max=None): :param z: Z-coordinate of the line :return: """ - y = y.asNumber(m) - z = z.asNumber(m) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) conc_x_inst = C.sel(y=y, z=z, method="nearest").max(dim="time") x_array = conc_x_inst.squeeze().x @@ -257,8 +257,8 @@ def plotMaxConcentrationOverTime_noQ(self,C, y, z, x_min=None,x_max=None): plt.title(f"Maximum concentration over time. y={y}[m], z={z}[m]") plt.grid() if x_min is not None and x_max is not None: - x_min = x_min.asNumber(m) - x_max = x_max.asNumber(m) + x_min = unumToPint(x_min).m_as(ureg.m) + x_max = unumToPint(x_max).m_as(ureg.m) plt.xlim(x_min, x_max) plt.show() @@ -273,9 +273,9 @@ def plotFixedPointConcentrationOverTime_noQ(self, C, x, y, z, t_min=None,t_max=N :return: """ - x = x.asNumber(m) - y = y.asNumber(m) - z = z.asNumber(m) + x = unumToPint(x).m_as(ureg.m) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) conc_inst_t = C.sel(x=x, y=y, z=z, method="nearest") time_array = conc_inst_t.squeeze().time @@ -286,8 +286,8 @@ def plotFixedPointConcentrationOverTime_noQ(self, C, x, y, z, t_min=None,t_max=N plt.title(f"Receptor at x={x}[m], y={y}[m], z={z}[m].") plt.grid() if t_min is not None and t_max is not None: - t_min = t_min.asNumber(min) - t_max = t_max.asNumber(min) + t_min = unumToPint(t_min).m_as(ureg.min) + t_max = unumToPint(t_max).m_as(ureg.min) plt.xlim(t_min, t_max) plt.show() @@ -301,9 +301,9 @@ def plotTIACPerDistance_noQ(self,TIAC, y, z, time): :param z: Z-coordinate of the line :return: """ - y = y.asNumber(m) - z = z.asNumber(m) - time = time.asNumber(min) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) + time = unumToPint(time).m_as(ureg.min) dos_x_inst = TIAC.sel(y=y, z=z, time=time, method="nearest") x_array = dos_x_inst .squeeze().x @@ -326,9 +326,9 @@ def plotFixedPointTIACOverTime_noQ(self, TIAC, x, y, z): :return: """ - x = x.asNumber(m) - y = y.asNumber(m) - z = z.asNumber(m) + x = unumToPint(x).m_as(ureg.m) + y = unumToPint(y).m_as(ureg.m) + z = unumToPint(z).m_as(ureg.m) dos_inst_t = TIAC.sel(x=x, y=y, z=z, method="nearest") time_array = dos_inst_t .squeeze().time diff --git a/hera/simulations/hydrodynamics/nearWallFlow.py b/hera/simulations/hydrodynamics/nearWallFlow.py index 674291297..1ff6d3b4b 100644 --- a/hera/simulations/hydrodynamics/nearWallFlow.py +++ b/hera/simulations/hydrodynamics/nearWallFlow.py @@ -4,9 +4,9 @@ Taken from Boundary-Layer Theory 9th edition (schlichting) pg 537. """ -from ...utils.unum import tounit +from hera.utils import tounit, tonumber +from hera.utils.unitHandler import ureg, unumToPint from scipy.optimize import fsolve -from unum.units import * import numpy class functionG: @@ -57,7 +57,7 @@ def technical_roughness(self): @property def height(self): - return self._channelHeight.asNumber(m) + return self._channelHeight.m_as(ureg.m) @property def hydraulicHeight(self): @@ -69,7 +69,7 @@ def technical_roughness_plus(self,ustar): :param ustar: :return: """ - return (ustar*self.technical_roughness/self.kinematicViscosity).asNumber() + return (ustar*self.technical_roughness/self.kinematicViscosity).magnitude def Cplus(self,ustar): @@ -83,7 +83,7 @@ def Cplus(self,ustar): The C^+. """ kappa = 0.41 # von karman constant. - kplus = (ustar*self.technical_roughness/self.kinematicViscosity).asNumber() + kplus = (ustar*self.technical_roughness/self.kinematicViscosity).magnitude if kplus <= 5: ret = 5 @@ -114,8 +114,8 @@ def __init__(self,Ra,nu): nu: float/ unit The viscosity of the fluid. default unit [m^2/s] """ - self._nu = tounit(nu,m**2/s) - self._Ra = tounit(Ra,m) + self._nu = tounit(nu, ureg.m**2/ureg.s) + self._Ra = tounit(Ra, ureg.m) @@ -148,7 +148,7 @@ def __init__(self,Ra,nu,channelHeight): """ super().__init__(Ra,nu) - self._channelHeight = tounit(channelHeight,m) + self._channelHeight = tounit(channelHeight, ureg.m) self._functionG = functionG() self.C_bar_plus_C_barbar = -1.7 # Cbar + Cbar bar (equation 17.91). self.C_bar = 0.94 # Cbar + Cbar bar (equation 17.91). @@ -187,8 +187,8 @@ def Re_tau(self,ustar): :param ustar: :return: """ - ustar = tounit(ustar, m / s) - return ((ustar*self._channelHeight/2)/self.kinematicViscosity).asNumber() + ustar = tounit(ustar, ureg.m/ureg.s) + return ((ustar*self._channelHeight/2)/self.kinematicViscosity).magnitude def get_Umean_from_Ustar(self,ustar): @@ -207,7 +207,7 @@ def get_Umean_from_Ustar(self,ustar): """ kappa = 0.41 # von karman constant. - ustar = tounit(ustar,m/s) + ustar = tounit(ustar, ureg.m/ureg.s) Re_tau = self.Re_tau(ustar) @@ -235,7 +235,7 @@ def get_Ucenter_from_Ustar(self,ustar): """ kappa = 0.41 # von karman constant. - ustar = tounit(ustar,m/s) + ustar = tounit(ustar, ureg.m/ureg.s) Re_tau = self.Re_tau(ustar) @@ -275,7 +275,7 @@ def __init__(self, Ra, nu, channelHeight): """ super().__init__(Ra,nu) - self._channelHeight = tounit(channelHeight, m) + self._channelHeight = tounit(channelHeight, ureg.m) self._functionG = functionG() @@ -310,8 +310,8 @@ def Re_tau(self, ustar): :param ustar: :return: """ - ustar = tounit(ustar, m / s) - return ((ustar * self._channelHeight / 2) / self.kinematicViscosity).asNumber() + ustar = tounit(ustar, ureg.m/ureg.s) + return ((ustar * self._channelHeight / 2) / self.kinematicViscosity).magnitude def get_Um_from_Ustar(self, ustar): """ @@ -331,7 +331,7 @@ def get_Um_from_Ustar(self, ustar): """ kappa = 0.41 # von karman constant. - ustar = tounit(ustar, m / s) + ustar = tounit(ustar, ureg.m/ureg.s) Re_tau = self.Re_tau(ustar) diff --git a/hera/utils/matplotlibCountour.py b/hera/utils/matplotlibCountour.py index 032164895..654503a85 100644 --- a/hera/utils/matplotlibCountour.py +++ b/hera/utils/matplotlibCountour.py @@ -24,14 +24,14 @@ def toGeopandas(ContourData, inunits): A geopandas object with the contours as polygons and levels as attributes. """ - from hera.utils.unitHandler import ureg + from hera.utils.unitHandler import ureg, unumToPint try: from shapely import geometry import geopandas except ImportError: print("gis support not installed. ") - units_conversion = inunits.asNumber(ureg.m) + units_conversion = unumToPint(inunits).m_as(ureg.m) polyList = [] levelsList = [] for col, level in zip(ContourData.collections, ContourData.levels): From b093712a373e843c5bc4d2fe2a71efe8ed668f6b Mon Sep 17 00:00:00 2001 From: Ilay Falach Date: Sun, 17 May 2026 13:34:37 +0300 Subject: [PATCH 2/5] Migrate hera/utils/jsonutils.py from unum to pint (#875) Remove all `from unum import Unum` local imports used only for isinstance checks. Replace `isinstance(x, Unum)` with `hasattr(x, 'asNumber')` duck typing so the production codebase no longer depends on the unum library at all. The three affected functions (ConfigurationToJSON, JSONToConfiguration, stripConfigurationUnits) all retain backwards compatibility: Unum objects passed in are still handled via unumToPint() conversion. Co-Authored-By: Claude Sonnet 4.6 --- hera/utils/jsonutils.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/hera/utils/jsonutils.py b/hera/utils/jsonutils.py index 22bda7fd1..f645215bd 100644 --- a/hera/utils/jsonutils.py +++ b/hera/utils/jsonutils.py @@ -61,8 +61,6 @@ def ConfigurationToJSON(valueToProcess, standardize=False, splitUnits=False, kee dict with all the values as string """ - from unum import Unum - logger = logging.getLogger("hera.bin.ConfigurationToJSON") ret = {} logger.debug(f"Processing {valueToProcess}") @@ -79,7 +77,7 @@ def ConfigurationToJSON(valueToProcess, standardize=False, splitUnits=False, kee elif "'" in str(valueToProcess): logger.debug(f"\t {valueToProcess} is String, use as is") ret = valueToProcess - elif isinstance(valueToProcess,Unum): + elif hasattr(valueToProcess, 'asNumber'): # Unum object — convert via unumToPint for backwards compat logger.debug(f"\t Try to convert *{valueToProcess}* to string") @@ -142,7 +140,6 @@ def JSONToConfiguration(valueToProcess,returnUnum=False,returnStandardize=False) A dict with the values convected to unum. """ - from unum import Unum logger = logging.getLogger("hera.util.JSONToConfiguration") logger.debug(f"Processing {valueToProcess} of type {type(valueToProcess)}") @@ -162,23 +159,23 @@ def JSONToConfiguration(valueToProcess,returnUnum=False,returnStandardize=False) ret = [JSONToConfiguration(x,returnUnum=returnUnum,returnStandardize=returnStandardize) for x in valueToProcess] elif isinstance(valueToProcess,(int, float, bool)): ret = valueToProcess - elif isinstance(valueToProcess,Unum): + elif hasattr(valueToProcess, 'asNumber'): # Unum object — backwards compat + from hera.utils.unitHandler import unumToPint, pintToUnum ret = valueToProcess if returnUnum else unumToPint(valueToProcess) else: # we want to avoid importing pint when unnecessary since it takes a long time to load from pint import Quantity if isinstance(valueToProcess,Quantity): + from hera.utils.unitHandler import pintToUnum ret = pintToUnum(valueToProcess) if returnUnum else valueToProcess elif "'" in str(valueToProcess): logger.debug(f"\t {valueToProcess} is a String, use as is") ret = valueToProcess else: print(valueToProcess) - logger.debug(f"\t Try to convert {valueToProcess} to unum") + logger.debug(f"\t Try to convert {valueToProcess} to pint") try: - from unum import Unum - - from hera.utils.unitHandler import unumToPint, ureg, pintToUnum, UndefinedUnitError, DimensionalityError + from hera.utils.unitHandler import ureg, pintToUnum, UndefinedUnitError, DimensionalityError ret = ureg(valueToProcess) ret = pintToUnum(ret) if returnUnum else ret except (UndefinedUnitError, DimensionalityError, ValueError): @@ -202,7 +199,6 @@ def stripConfigurationUnits(valueToProcess,returnStandardize=False, ignoreStanda Same json with all units stripped leaving just magnitudes """ - from unum import Unum from pint import Quantity from hera.utils.unitHandler import unumToPint @@ -218,7 +214,7 @@ def stripConfigurationUnits(valueToProcess,returnStandardize=False, ignoreStanda logger.debug("\t list, transforming to str every item") ret = [stripConfigurationUnits(x,returnStandardize=returnStandardize, ignoreStandardization=ignoreStandardization) for x in valueToProcess] - elif isinstance(valueToProcess,Unum): + elif hasattr(valueToProcess, 'asNumber'): # Unum object — backwards compat logger.debug(f"\t Try to convert *{valueToProcess}* to string") origPintObj = unumToPint(valueToProcess) pintObj = origPintObj.to_base_units() if (returnStandardize) else origPintObj From 75d64f71eb90b64ba0328bee6931520b0c3808b6 Mon Sep 17 00:00:00 2001 From: yehudaa Date: Mon, 18 May 2026 10:20:20 +0300 Subject: [PATCH 3/5] Fixing the LSM toolkit. --- hera/simulations/LSM/singleSimulation.py | 2 +- hera/simulations/LSM/template.py | 41 +++++++++++++++++------- 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/hera/simulations/LSM/singleSimulation.py b/hera/simulations/LSM/singleSimulation.py index 787c2d993..89a1bcc54 100644 --- a/hera/simulations/LSM/singleSimulation.py +++ b/hera/simulations/LSM/singleSimulation.py @@ -97,7 +97,7 @@ def getConcentration(self, Q=1*ureg.kg, time_units=ureg.min, q_units=ureg.mg): finalxarray = self.getDosage(Q=Q, time_units=time_units, q_units=q_units) dDosage = finalxarray['Dosage'].diff('datetime').to_dataset().rename({'Dosage': 'dDosage'}) - dDosage['C'] = dDosage['dDosage'] / finalxarray.attrs['dt'] + dDosage['C'] = dDosage['dDosage'] / finalxarray.attrs['dt'].m_as(time_units) dDosage.attrs = finalxarray.attrs return dDosage diff --git a/hera/simulations/LSM/template.py b/hera/simulations/LSM/template.py index 7c444c4d2..d21183df6 100644 --- a/hera/simulations/LSM/template.py +++ b/hera/simulations/LSM/template.py @@ -160,6 +160,11 @@ def run(self,topography=None, stations=None,canopy=None,params=dict(),deposition f"{toolkit.TOOLKIT_SAVEMODE_FILEANDDB_REPLACE} in order to replace it.") if saveMode == toolkit.TOOLKIT_SAVEMODE_FILEANDDB_REPLACE: docList.delete() + + curr_counter = self.Toolkit.getCounterAndAdd(f"LSM_SIMULATION_{self.Toolkit.projectName}") + simulationName = simulationName if simulationName is not None else f"LSM_Simulation_{curr_counter}" + saveDir = os.path.join(saveDir, simulationName) + doc = self.Toolkit.addSimulationsDocument( type=self.doctype_simulation, resource='None', @@ -170,8 +175,7 @@ def run(self,topography=None, stations=None,canopy=None,params=dict(),deposition simulationName=simulationName, params=updated_params) ) - curr_counter = self.Toolkit.getCounterAndAdd(f"LSM_SIMULATION_{self.Toolkit.projectName}") - saveDir = os.path.join(saveDir, f"LSM_Simulation_{curr_counter}") + if self.to_xarray: doc['resource'] = os.path.join(saveDir, 'netcdf', '*') doc['dataFormat'] = datatypes.NETCDF_XARRAY @@ -315,15 +319,19 @@ def run(self,topography=None, stations=None,canopy=None,params=dict(),deposition @staticmethod def prepareParams(desc, paramsToPrepare): logger = get_logger(instance=None, name="hera.simulations.LSM.prepareParams") - if desc is not None and 'units' in desc: - for key in desc["units"].keys(): - param_item= paramsToPrepare[key] - if hasattr(param_item, 'asNumber') or hasattr(param_item, 'magnitude'): - paramsToPrepare[key] = unumToPint(param_item).m_as(ureg.parse_expression(desc["units"][key])) - else: - raise ValueError(f"parameters must use either pint or unum to specify units, currently type({param_item})={type(param_item)}") - if 'duration' in paramsToPrepare and isinstance(paramsToPrepare['duration'], Quantity): - paramsToPrepare['duration'] = paramsToPrepare['duration'].to(ureg.minutes) + try: + if desc is not None and 'units' in desc: + for key in desc["units"].keys(): + param_item= paramsToPrepare[key] + if hasattr(param_item, 'asNumber') or hasattr(param_item, 'magnitude'): + paramsToPrepare[key] = unumToPint(param_item).m_as(ureg.parse_expression(desc["units"][key])) + elif key=='duration': + paramsToPrepare[key] = param_item*ureg.minutes + else: + paramsToPrepare[key] = ureg.parse_expression(param_item).m_as(ureg.parse_expression(desc["units"][key])) + except: + raise ValueError(f"parameters must use either pint or unum to specify units, currently type({param_item})={type(param_item)}") + paramsToPrepare = stripConfigurationUnits(paramsToPrepare, returnStandardize=True, ignoreStandardization=["duration"]) for integer_field in ["TopoXn", "TopoYn"]: if not isinstance(paramsToPrepare[integer_field], int): @@ -435,6 +443,15 @@ def getSimulationByID(self,id): """ return SingleSimulation(self.getDocumentByID(id)) + def getSimulationByName(self,simulationName): + """ + get a simulation by document id + + :param id: + :return: + """ + return SingleSimulation(self.getSimulationsDocuments(type=self.doctype_simulation,templateName=self.templateName,simulationName=simulationName)) + def _getSimulationsList(self, **query): """ List the Simulation parameters that fulfil the query @@ -482,4 +499,4 @@ def getSimulationsTable(self,**query): return df except ValueError: - raise FileNotFoundError('No simulations.old found') + raise FileNotFoundError('No simulations found') From 40a13e7591c3da4b85b3624069d14232ce987046 Mon Sep 17 00:00:00 2001 From: yehudaa Date: Tue, 19 May 2026 10:21:16 +0300 Subject: [PATCH 4/5] Fixing the LSM toolkit. Fixing the shaply2 bug in the code --- .../agents/effects/Calculator.py | 23 ++++++++----------- .../presentation/casualtiesFigs.py | 7 +++--- hera/utils/matplotlibCountour.py | 4 +++- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/hera/riskassessment/agents/effects/Calculator.py b/hera/riskassessment/agents/effects/Calculator.py index bc1b808b6..e6d45c096 100755 --- a/hera/riskassessment/agents/effects/Calculator.py +++ b/hera/riskassessment/agents/effects/Calculator.py @@ -2,8 +2,10 @@ import pandas from hera.utils import ureg import xarray +from utils.unitHandler import unumToPint -class AbstractCalculator(object): + +class AbstractCalculator(object): """ Holds the abstract calculator. """ @@ -88,11 +90,10 @@ def calculate(self,concentrationField,breathingRate=10*ureg.L/ureg.min,time="dat """ + breathingRate = unumToPint(breathingRate) breathingRatio = (breathingRate/self.injuryBreathingRate).magnitude if inUnits is None: - if hasattr(concentrationField, "attrs"): - attrs_units = concentrationField.attrs.get("field", None) - inUnits = concentrationField.attrs[field] if attrs_units is not None else ureg.mg / ureg.m ** 3 + inUnits = concentrationField.attrs[field] if hasattr(concentrationField, "attrs") else 1*(ureg.mg / ureg.m ** 3) CunitConversion = inUnits.m_as(ureg.mg / ureg.m ** 3) @@ -165,12 +166,10 @@ def calculate(self,concentrationField,field,breathingRate=10*ureg.L/ureg.min,tim """ + breathingRate = unumToPint(breathingRate) breathingRatio = (breathingRate/self.injuryBreathingRate).magnitude - if inUnits is None: - if hasattr(concentrationField, "attrs"): - attrs_units = concentrationField.attrs.get("field", None) - inUnits = concentrationField.attrs[field] if attrs_units is not None else ureg.mg / ureg.m ** 3 + inUnits = concentrationField.attrs[field] if hasattr(concentrationField, "attrs") else 1*(ureg.mg / ureg.m ** 3) CunitConversion = inUnits.m_as(ureg.mg / ureg.m ** 3) if isinstance(concentrationField, xarray.Dataset): @@ -236,12 +235,10 @@ def calculate(self,concentrationField,field,breathingRate=10*ureg.L/ureg.min,tim """ - breathingRatio = (breathingRate / self.injuryBreathingRate).magnitude - + breathingRate = unumToPint(breathingRate) + breathingRatio = (breathingRate/self.injuryBreathingRate).magnitude if inUnits is None: - if hasattr(concentrationField, "attrs"): - attrs_units = concentrationField.attrs.get("field", None) - inUnits = attrs_units[field] if attrs_units is not None else ureg.mg / ureg.m ** 3 + inUnits = concentrationField.attrs[field] if hasattr(concentrationField, "attrs") else 1*(ureg.mg / ureg.m ** 3) CunitConversion = inUnits.m_as(ureg.mg / ureg.m ** 3) if isinstance(concentrationField, xarray.Dataset): diff --git a/hera/riskassessment/presentation/casualtiesFigs.py b/hera/riskassessment/presentation/casualtiesFigs.py index 9d738e72e..9f3ffaebd 100644 --- a/hera/riskassessment/presentation/casualtiesFigs.py +++ b/hera/riskassessment/presentation/casualtiesFigs.py @@ -202,9 +202,10 @@ def plotCasualtiesProjection(self, for severity,prop,lineprop in zip(severityList,cycler,boundarycycler): if severity not in projected.index: continue - if projected.loc[severity].geometry.filterType == 'GeometryCollection' or projected.loc[severity].geometry.filterType == 'MultiPolygon': - for pol in projected.loc[severity].geometry: - if pol.filterType == 'LineString': + geom =projected.loc[severity].geometry + if hasattr(geom,"geoms"): + for pol in geom.geoms: + if pol.geom_type == 'LineString': ax.plot(*pol.xy,**lineprop) else: ax.add_patch(PolygonPatch(pol,**prop) ) diff --git a/hera/utils/matplotlibCountour.py b/hera/utils/matplotlibCountour.py index 654503a85..0744dbc1c 100644 --- a/hera/utils/matplotlibCountour.py +++ b/hera/utils/matplotlibCountour.py @@ -9,7 +9,7 @@ def standardize_polygon(poly, units_conversion): return [(x * units_conversion, y * units_conversion) for (x, y) in zip(xs, ys)] -def toGeopandas(ContourData, inunits): +def toGeopandas(ContourData, inunits=None): """ Converts the contours of matplotlib to polygons. @@ -31,6 +31,8 @@ def toGeopandas(ContourData, inunits): import geopandas except ImportError: print("gis support not installed. ") + inunits = inunits if inunits is None else 1*ureg.m + units_conversion = unumToPint(inunits).m_as(ureg.m) polyList = [] levelsList = [] From be491a4b8d1ac80d49f03605a03ac90c2e4cc201 Mon Sep 17 00:00:00 2001 From: yehudaa Date: Tue, 19 May 2026 19:36:11 +0300 Subject: [PATCH 5/5] LSM : - fixed the plots of the rose and the projections. utils - added plots of polygons. --- hera/riskassessment/agents/effects/Injury.py | 9 ++- .../presentation/casualtiesFigs.py | 22 +++---- hera/utils/matplotlibCountour.py | 65 +++++++++++++++++++ 3 files changed, 78 insertions(+), 18 deletions(-) diff --git a/hera/riskassessment/agents/effects/Injury.py b/hera/riskassessment/agents/effects/Injury.py index 8bbcaed7f..eecd8e25f 100755 --- a/hera/riskassessment/agents/effects/Injury.py +++ b/hera/riskassessment/agents/effects/Injury.py @@ -195,7 +195,7 @@ def _postCalculatePointWise(self, retList): pass def calculate(self, concentrationField, field, time="datetime", x="x", y="y", breathingRate=10 * ureg.L / ureg.min, - **parameters): + sel={},isel={}): """Deprecated. Use ``calculateRegionOfInjured`` instead.""" import warnings warnings.warn("This function is obselete. Use calculateRegionOfInjured instead") @@ -205,9 +205,9 @@ def calculate(self, concentrationField, field, time="datetime", x="x", y="y", br x=x, y=y, breathingRate=breathingRate, - **parameters) + sel=sel,isel=isel) - def calculateRegionOfInjured(self,concentrationField,field,time="datetime",x="x",y="y",breathingRate=10*ureg.L / ureg.min,**parameters): + def calculateRegionOfInjured(self,concentrationField,field,time="datetime",x="x",y="y",breathingRate=10*ureg.L / ureg.min,sel={},isel={}): """ Calculates the fraction of the population that was effected in each point, and returns its contour. The levels are determined by the injury type. @@ -247,8 +247,7 @@ def calculateRegionOfInjured(self,concentrationField,field,time="datetime",x="x" toxicLoads = self.calculateToxicLoads(concentrationField=concentrationField, time=time, breathingRate=breathingRate, - field=field, - **parameters).compute() + field=field).sel(**sel).isel(**isel).compute() for lvl in self.levels: data = lvl.calculateContours(toxicLoads=toxicLoads,time=time,x=x,y=y) if data is not None: diff --git a/hera/riskassessment/presentation/casualtiesFigs.py b/hera/riskassessment/presentation/casualtiesFigs.py index 9f3ffaebd..4b6bcc04e 100644 --- a/hera/riskassessment/presentation/casualtiesFigs.py +++ b/hera/riskassessment/presentation/casualtiesFigs.py @@ -60,7 +60,7 @@ def plotCasualtiesRose(self, projectedData.append(injuryareas) projectedData = pandas.concat(projectedData) - pivotedData = projectedData.pivot("angle",'severity',effectedPopulation).reset_index().fillna(0) + pivotedData = projectedData.pivot(index="angle",columns='severity',values=effectedPopulation).reset_index().fillna(0) if (ax is None): fig = plt.gcf() ax = fig.add_subplot(111,polar=True) @@ -83,15 +83,15 @@ def plotCasualtiesRose(self, cycler += plt.cycler(width=[weights]*len(severityList)) bottom = numpy.zeros(pivotedData.shape[0]) + pivotedData["total"] = 0 for severity,plotprops in zip(severityList,cycler): if severity not in pivotedData.columns: continue ax.bar(pivotedData["angle"],pivotedData[severity],label=severity,bottom=bottom,**plotprops) bottom += pivotedData[severity] - pivotedData["total"] = 0 - for severity in severityList: pivotedData["total"] += pivotedData[severity] + if windDistribution is not None: windDist = windDistribution.copy() windTicks = [25,50,75,100] if windTicks is None else sorted(windTicks) @@ -175,7 +175,9 @@ def plotCasualtiesProjection(self, :cycler: a property cycler for the polygons. :boundarycycler: a property cycler for the overall polygons. """ - if (ax is None): + from hera.utils.matplotlibCountour import plot_polygons + + if (ax is None): fig = plt.gcf() ax = fig.add_subplot(111) elif isinstance(ax,list): @@ -196,21 +198,15 @@ def plotCasualtiesProjection(self, projected = retProj.dissolve("severity") boundarycycler = plt.cycler(color=plt.rcParams['axes.prop_cycle'].by_key()['color']) if boundarycycler is None else boundarycycler - cycler = plt.cycler(fc=plt.rcParams['axes.prop_cycle'].by_key()['color'])*plt.cycler(ec=['None']) if cycler is None else cycler + cycler = plt.cycler(facecolor=plt.rcParams['axes.prop_cycle'].by_key()['color'])*plt.cycler(edgecolor=['None']) if cycler is None else cycler patchList = [] for severity,prop,lineprop in zip(severityList,cycler,boundarycycler): if severity not in projected.index: continue geom =projected.loc[severity].geometry - if hasattr(geom,"geoms"): - for pol in geom.geoms: - if pol.geom_type == 'LineString': - ax.plot(*pol.xy,**lineprop) - else: - ax.add_patch(PolygonPatch(pol,**prop) ) - else: - ax.add_patch(PolygonPatch(projected.loc[severity].geometry,**prop) ) + plot_polygons(ax,geom,**prop) + for ((severity,severitydata),prop) in zip(results.shiftLocationAndAngle(loc,mathematical_angle=rotate_angle,geometry="TotalPolygon") .query("severity in %s" % numpy.atleast_1d(plumSeverity)) diff --git a/hera/utils/matplotlibCountour.py b/hera/utils/matplotlibCountour.py index 0744dbc1c..994ff5c8f 100644 --- a/hera/utils/matplotlibCountour.py +++ b/hera/utils/matplotlibCountour.py @@ -56,4 +56,69 @@ def toGeopandas(ContourData, inunits=None): ret = geopandas.GeoDataFrame({"Level": levelsList, "contour": polyList}, geometry="contour") return ret +def plot_polygons(ax, geom, facecolor="none", edgecolor="black", linewidth=1): + import matplotlib.pyplot as plt + from matplotlib.patches import Polygon as MplPolygon + from matplotlib.collections import PatchCollection + from shapely.geometry import ( + Polygon, MultiPolygon, + LineString, MultiLineString, + Point, MultiPoint, + GeometryCollection + ) + + if geom.is_empty: + return + + if isinstance(geom, Polygon): + # exterior + patch = MplPolygon( + list(geom.exterior.coords), + closed=True, + facecolor=facecolor, + edgecolor=edgecolor, + linewidth=linewidth, + ) + ax.add_patch(patch) + + # holes + for interior in geom.interiors: + x, y = interior.xy + ax.plot(x, y, color=edgecolor, linewidth=linewidth) + + elif isinstance(geom, MultiPolygon): + for part in geom.geoms: + plot_geom(ax, part, facecolor, edgecolor, linewidth) + + elif isinstance(geom, LineString): + x, y = geom.xy + ax.plot(x, y, color=edgecolor, linewidth=linewidth) + + elif isinstance(geom, MultiLineString): + for part in geom.geoms: + plot_geom(ax, part, facecolor, edgecolor, linewidth) + + elif isinstance(geom, Point): + ax.plot(geom.x, geom.y, "o", color=edgecolor) + + elif isinstance(geom, MultiPoint): + for part in geom.geoms: + plot_geom(ax, part, facecolor, edgecolor, linewidth) + + elif isinstance(geom, GeometryCollection): + for part in geom.geoms: + plot_geom(ax, part, facecolor, edgecolor, linewidth) + + else: + raise TypeError(f"Unsupported geometry type: {geom.geom_type}") + + +# usage +fig, ax = plt.subplots() + +plot_geom(ax, my_geometry, facecolor="lightblue", edgecolor="black") + +ax.set_aspect("equal") +ax.autoscale() +plt.show()