Skip to content

Commit

Permalink
Merge pull request #2064 from jthielen/refactor-unit-handling-bottlen…
Browse files Browse the repository at this point in the history 
…ecks

Improve some thermo calculation bottlenecks by refactoring how units are handled
  • Loading branch information
@dopplershift
dopplershift committed Jan 21, 2022
2 parents 10fea45 + 588f4e2 commit 94ef0e0
Show file tree
Hide file tree
Showing 8 changed files with 420 additions and 242 deletions.
2 changes: 1 addition & 1 deletion setup.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ per-file-ignores = examples/*.py: D MPY001 T003 T001
src/metpy/interpolate/*.py: RST306
src/metpy/io/*.py: RST306
src/metpy/future.py: RST307
src/metpy/constants.py: RST306
src/metpy/constants/*.py: RST306
docs/doc-server.py: T001
tests/*.py: MPY001
ci/filter_links.py: E731 T001
Expand Down
206 changes: 130 additions & 76 deletions src/metpy/calc/thermo.py
View file

Large diffs are not rendered by default.

61 changes: 6 additions & 55 deletions src/metpy/constants.py → src/metpy/constants/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Copyright (c) 2008,2015,2016,2018 MetPy Developers.
# Copyright (c) 2008,2015,2016,2018,2021 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause

Expand Down Expand Up @@ -70,59 +70,10 @@
.. [8] [Picard2008]_
""" # noqa: E501

from .package_tools import Exporter
from .units import units
from . import nounit # noqa: F401
from .default import * # noqa: F403
from ..package_tools import set_module

exporter = Exporter(globals())
__all__ = default.__all__[:] # pylint: disable=undefined-variable

# Export all the variables defined in this block
with exporter:
# Earth
earth_gravity = g = units.Quantity(9.80665, 'm / s^2')
Re = earth_avg_radius = units.Quantity(6371008.7714, 'm')
G = gravitational_constant = units.Quantity(6.67430e-11, 'm^3 / kg / s^2')
GM = geocentric_gravitational_constant = units.Quantity(3986005e8, 'm^3 / s^2')
omega = earth_avg_angular_vel = units.Quantity(7292115e-11, 'rad / s')
d = earth_sfc_avg_dist_sun = units.Quantity(149597870700., 'm')
S = earth_solar_irradiance = units.Quantity(1360.8, 'W / m^2')
delta = earth_max_declination = units.Quantity(23.45, 'degrees')
earth_orbit_eccentricity = units.Quantity(0.0167, 'dimensionless')
earth_mass = me = geocentric_gravitational_constant / gravitational_constant

# molar gas constant
R = units.Quantity(8.314462618, 'J / mol / K')

# Water
Mw = water_molecular_weight = units.Quantity(18.015268, 'g / mol')
Rv = water_gas_constant = R / Mw
rho_l = density_water = units.Quantity(999.97495, 'kg / m^3')
wv_specific_heat_ratio = units.Quantity(1.330, 'dimensionless')
Cp_v = wv_specific_heat_press = (
wv_specific_heat_ratio * Rv / (wv_specific_heat_ratio - 1)
)
Cv_v = wv_specific_heat_vol = Cp_v / wv_specific_heat_ratio
Cp_l = water_specific_heat = units.Quantity(4.2194, 'kJ / kg / K')
Lv = water_heat_vaporization = units.Quantity(2.50084e6, 'J / kg')
Lf = water_heat_fusion = units.Quantity(3.337e5, 'J / kg')
Cp_i = ice_specific_heat = units.Quantity(2090, 'J / kg / K')
rho_i = density_ice = units.Quantity(917, 'kg / m^3')

# Dry air
Md = dry_air_molecular_weight = units.Quantity(28.96546e-3, 'kg / mol')
Rd = dry_air_gas_constant = R / Md
dry_air_spec_heat_ratio = units.Quantity(1.4, 'dimensionless')
Cp_d = dry_air_spec_heat_press = (
dry_air_spec_heat_ratio * Rd / (dry_air_spec_heat_ratio - 1)
)
Cv_d = dry_air_spec_heat_vol = Cp_d / dry_air_spec_heat_ratio
rho_d = dry_air_density_stp = (
units.Quantity(1000., 'mbar') / (Rd * units.Quantity(273.15, 'K'))
).to('kg / m^3')

# General meteorology constants
P0 = pot_temp_ref_press = units.Quantity(1000., 'mbar')
kappa = poisson_exponent = (Rd / Cp_d).to('dimensionless')
gamma_d = dry_adiabatic_lapse_rate = g / Cp_d
epsilon = molecular_weight_ratio = (Mw / Md).to('dimensionless')

del Exporter
set_module(globals())
62 changes: 62 additions & 0 deletions src/metpy/constants/default.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
# Copyright (c) 2008,2015,2016,2018,2021 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Constant and thermophysical property values expressed as quantities."""

from ..package_tools import Exporter
from ..units import units

exporter = Exporter(globals())

# Export all the variables defined in this block
with exporter:
# Earth
earth_gravity = g = units.Quantity(9.80665, 'm / s^2')
Re = earth_avg_radius = units.Quantity(6371008.7714, 'm')
G = gravitational_constant = units.Quantity(6.67430e-11, 'm^3 / kg / s^2')
GM = geocentric_gravitational_constant = units.Quantity(3986005e8, 'm^3 / s^2')
omega = earth_avg_angular_vel = units.Quantity(7292115e-11, 'rad / s')
d = earth_sfc_avg_dist_sun = units.Quantity(149597870700., 'm')
S = earth_solar_irradiance = units.Quantity(1360.8, 'W / m^2')
delta = earth_max_declination = units.Quantity(23.45, 'degrees')
earth_orbit_eccentricity = units.Quantity(0.0167, 'dimensionless')
earth_mass = me = geocentric_gravitational_constant / gravitational_constant

# molar gas constant
R = units.Quantity(8.314462618, 'J / mol / K')

# Water
Mw = water_molecular_weight = units.Quantity(18.015268, 'g / mol')
Rv = water_gas_constant = R / Mw
rho_l = density_water = units.Quantity(999.97495, 'kg / m^3')
wv_specific_heat_ratio = units.Quantity(1.330, 'dimensionless')
Cp_v = wv_specific_heat_press = (
wv_specific_heat_ratio * Rv / (wv_specific_heat_ratio - 1)
)
Cv_v = wv_specific_heat_vol = Cp_v / wv_specific_heat_ratio
Cp_l = water_specific_heat = units.Quantity(4.2194, 'kJ / kg / K')
Lv = water_heat_vaporization = units.Quantity(2.50084e6, 'J / kg')
Lf = water_heat_fusion = units.Quantity(3.337e5, 'J / kg')
Cp_i = ice_specific_heat = units.Quantity(2090, 'J / kg / K')
rho_i = density_ice = units.Quantity(917, 'kg / m^3')
sat_pressure_0c = units.Quantity(6.112, 'millibar')

# Dry air
Md = dry_air_molecular_weight = units.Quantity(28.96546e-3, 'kg / mol')
Rd = dry_air_gas_constant = R / Md
dry_air_spec_heat_ratio = units.Quantity(1.4, 'dimensionless')
Cp_d = dry_air_spec_heat_press = (
dry_air_spec_heat_ratio * Rd / (dry_air_spec_heat_ratio - 1)
)
Cv_d = dry_air_spec_heat_vol = Cp_d / dry_air_spec_heat_ratio
rho_d = dry_air_density_stp = (
units.Quantity(1000., 'mbar') / (Rd * units.Quantity(273.15, 'K'))
).to('kg / m^3')

# General meteorology constants
P0 = pot_temp_ref_press = units.Quantity(1000., 'mbar')
kappa = poisson_exponent = (Rd / Cp_d).to('dimensionless')
gamma_d = dry_adiabatic_lapse_rate = g / Cp_d
epsilon = molecular_weight_ratio = (Mw / Md).to('dimensionless')

del Exporter
16 changes: 16 additions & 0 deletions src/metpy/constants/nounit.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
# Copyright (c) 2021 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""Subset of constant and thermophysical property values expressed as floats in base units."""

from . import default
from ..units import units

Rd = default.Rd.m_as('m**2 / K / s**2')
Lv = default.Lv.m_as('m**2 / s**2')
Cp_d = default.Cp_d.m_as('m**2 / K / s**2')
zero_degc = units.Quantity(0., 'degC').m_as('K')
sat_pressure_0c = default.sat_pressure_0c.m_as('Pa')
epsilon = default.epsilon.m_as('')
kappa = default.kappa.m_as('')
g = default.g.m_as('m / s**2')
113 changes: 54 additions & 59 deletions src/metpy/io/gempak.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
import contextlib
from datetime import datetime, timedelta
from enum import Enum
from itertools import product
from itertools import product, repeat
import logging
import math
import struct
Expand All @@ -21,12 +21,10 @@

from ._tools import IOBuffer, NamedStruct, open_as_needed
from .. import constants
from ..calc import (mixing_ratio_from_specific_humidity, scale_height,
specific_humidity_from_dewpoint, thickness_hydrostatic,
from ..calc import (scale_height, specific_humidity_from_dewpoint, thickness_hydrostatic,
virtual_temperature)
from ..package_tools import Exporter
from ..plots.mapping import CFProjection
from ..units import units

exporter = Exporter(globals())
log = logging.getLogger(__name__)
Expand Down Expand Up @@ -171,6 +169,11 @@ class DataSource(Enum):
])


def convert_degc_to_k(val, missing=-9999):
"""Convert scalar values from degC to K, handling missing values."""
return val + constants.nounit.zero_degc if val != missing else val


def _data_source(source):
"""Get data source from stored integer."""
try:
Expand Down Expand Up @@ -277,34 +280,31 @@ def _interp_logp_height(sounding, missing=-9999):

if maxlev < size - 1:
if maxlev > -1:
pb = units.Quantity(sounding['PRES'][maxlev], 'hPa')
zb = units.Quantity(sounding['HGHT'][maxlev], 'm')
tb = units.Quantity(sounding['TEMP'][maxlev], 'degC')
tdb = units.Quantity(sounding['DWPT'][maxlev], 'degC')
pb = sounding['PRES'][maxlev] * 100 # hPa to Pa
zb = sounding['HGHT'][maxlev] # m
tb = convert_degc_to_k(sounding['TEMP'][maxlev], missing)
tdb = convert_degc_to_k(sounding['DWPT'][maxlev], missing)
else:
pb = units.Quantity(missing, 'hPa')
zb = units.Quantity(missing, 'm')
tb = units.Quantity(missing, 'degC')
tdb = units.Quantity(missing, 'degC')
pb, zb, tb, tdb = repeat(missing, 4)

for i in range(maxlev + 1, size):
if sounding['HGHT'][i] == missing:
tt = units.Quantity(sounding['TEMP'][i], 'degC')
tdt = units.Quantity(sounding['DWPT'][i], 'degC')
pt = units.Quantity(sounding['PRES'][i], 'hPa')
tt = convert_degc_to_k(sounding['TEMP'][i], missing)
tdt = convert_degc_to_k(sounding['DWPT'][i], missing)
pt = sounding['PRES'][i] * 100 # hPa to Pa

pl = units.Quantity([pb.m, pt.m], 'hPa')
tl = units.Quantity([tb.m, tt.m], 'degC')
tdl = units.Quantity([tdb.m, tdt.m], 'degC')
pl = np.array([pb, pt])
tl = np.array([tb, tt])
tdl = np.array([tdb, tdt])

if missing in tdl.m:
rl = None
if missing in tdl:
tvl = tl
else:
ql = specific_humidity_from_dewpoint(pl, tdl)
rl = mixing_ratio_from_specific_humidity(ql)
ql = specific_humidity_from_dewpoint._nounit(pl, tdl)
tvl = virtual_temperature._nounit(tl, ql)

if missing not in [*tl.m, zb.m]:
sounding['HGHT'][i] = (zb + thickness_hydrostatic(pl, tl, rl)).m
if missing not in [*tvl, zb]:
sounding['HGHT'][i] = (zb + thickness_hydrostatic._nounit(pl, tvl))
else:
sounding['HGHT'][i] = missing

Expand Down Expand Up @@ -352,7 +352,7 @@ def _interp_moist_height(sounding, missing=-9999):
subroutine in GEMPAK. This the default behavior when
merging observed sounding data.
"""
hlist = units.Quantity(np.ones(len(sounding['PRES'])) * -9999, 'm')
hlist = np.ones(len(sounding['PRES'])) * -9999

ilev = -1
top = False
Expand All @@ -370,12 +370,12 @@ def _interp_moist_height(sounding, missing=-9999):
found = True

while not top:
plev = units.Quantity(sounding['PRES'][ilev], 'hPa')
pb = units.Quantity(sounding['PRES'][ilev], 'hPa')
tb = units.Quantity(sounding['TEMP'][ilev], 'degC')
tdb = units.Quantity(sounding['DWPT'][ilev], 'degC')
zb = units.Quantity(sounding['HGHT'][ilev], 'm')
zlev = units.Quantity(sounding['HGHT'][ilev], 'm')
plev = sounding['PRES'][ilev] * 100 # hPa to Pa
pb = sounding['PRES'][ilev] * 100 # hPa to Pa
tb = convert_degc_to_k(sounding['TEMP'][ilev], missing)
tdb = convert_degc_to_k(sounding['DWPT'][ilev], missing)
zb = sounding['HGHT'][ilev] # m
zlev = sounding['HGHT'][ilev] # m
jlev = ilev
klev = 0
mand = False
Expand All @@ -386,37 +386,33 @@ def _interp_moist_height(sounding, missing=-9999):
mand = True
top = True
else:
pt = units.Quantity(sounding['PRES'][jlev], 'hPa')
tt = units.Quantity(sounding['TEMP'][jlev], 'degC')
tdt = units.Quantity(sounding['DWPT'][jlev], 'degC')
zt = units.Quantity(sounding['HGHT'][jlev], 'm')
if (zt.m != missing
and tt != missing):
pt = sounding['PRES'][jlev] * 100 # hPa to Pa
tt = convert_degc_to_k(sounding['TEMP'][jlev], missing)
tdt = convert_degc_to_k(sounding['DWPT'][jlev], missing)
zt = sounding['HGHT'][jlev] # m
if (zt != missing and tt != missing):
mand = True
klev = jlev

pl = units.Quantity([pb.m, pt.m], 'hPa')
tl = units.Quantity([tb.m, tt.m], 'degC')
tdl = units.Quantity([tdb.m, tdt.m], 'degC')
pl = np.array([pb, pt])
tl = np.array([tb, tt])
tdl = np.array([tdb, tdt])

if missing in tdl.m:
rl = None
if missing in tdl:
tvl = tl
else:
ql = specific_humidity_from_dewpoint(pl, tdl)
rl = mixing_ratio_from_specific_humidity(ql)
tvl = virtual_temperature(tl, ql)
ql = specific_humidity_from_dewpoint._nounit(pl, tdl)
tvl = virtual_temperature._nounit(tl, ql)

if missing not in [*tl.m, zb.m]:
scale_z = scale_height(*tvl)
znew = zb + thickness_hydrostatic(pl, tl, rl)
if missing not in [*tl, zb]:
scale_z = scale_height._nounit(*tvl)
znew = zb + thickness_hydrostatic._nounit(pl, tvl)
tb = tt
tdb = tdt
pb = pt
zb = znew
else:
scale_z = units.Quantity(missing, 'm')
znew = units.Quantity(missing, 'm')
scale_z, znew = repeat(missing, 2)
hlist[jlev] = scale_z

if klev != 0:
Expand All @@ -427,17 +423,16 @@ def _interp_moist_height(sounding, missing=-9999):
hbb = zlev
pbb = plev
for ii in range(ilev + 1, jlev):
p = units.Quantity(sounding['PRES'][ii], 'hPa')
p = sounding['PRES'][ii] * 100 # hPa to Pa
scale_z = hlist[ii]
if missing not in [scale_z.m, hbb.m,
pbb.m, p.m]:
th = ((scale_z * constants.g) / constants.Rd).to('degC')
tbar = units.Quantity([th.m, th.m], 'degC')
pll = units.Quantity([pbb.m, p.m], 'hPa')
z = hbb + thickness_hydrostatic(pll, tbar)
if missing not in [scale_z, hbb, pbb, p]:
th = (scale_z * constants.nounit.g) / constants.nounit.Rd
tbar = np.array([th, th])
pll = np.array([pbb, p])
z = hbb + thickness_hydrostatic._nounit(pll, tbar)
else:
z = units.Quantity(missing, 'm')
sounding['HGHT'][ii] = z.m
z = missing
sounding['HGHT'][ii] = z
hbb = z
pbb = p

Expand Down
Loading

0 comments on commit 94ef0e0

Please sign in to comment.