"""
A python module to parse VOUnit strings and compute conversion factors.
We believe this implements the full VOUnit specification.
To use this, you must have the gavo utils package installed. For
details, see http://soft.g-vo.org. The changes required to make this
work without gavo.utils are minute, though. If you want that, talk to
the authors.
Unit tests for this are at
http://svn.ari.uni-heidelberg.de/svn/gavo/python/trunk/tests/unitconvtest.py
"""
#c Copyright 2008-2023, the GAVO project <gavo@ari.uni-heidelberg.de>
#c
#c This program is free software, covered by the GNU GPL. See the
#c COPYING file in the source distribution.
import functools
import math
from gavo import utils
from gavo.utils import parsetricks
[docs]class IncompatibleUnits(utils.Error):
pass
[docs]class BadUnit(utils.Error):
pass
PLANCK_H = 6.62607015e-34
LIGHT_C = 299792458.0
BOLTZMANN_K = 1.380649e-23
# We can't yet restructure the tree, so we don't do SI-base-casting for
# compound units. Also, we don't change anything when a change of exponent
# would be necessary
PLAIN_UNITS = units = {
"a": (3600*24*365.25, "s"), # This is the SI/julian year!
"A": (1, "A"), # *Ampere* not Angstrom
"adu": (1, "adu"),
"\xc5": (1e-10, "m"),
"Angstrom": (1e-10, "m"),
"angstrom": (1e-10, "m"),
"AU": (1.495978701e11, "m"), # as per ESA/SP-1200
"arcmin": (math.pi/180./60., "rad"),
"arcsec": (math.pi/180./3600., "rad"),
"barn": (1, "barn"), # 1e-28 m2
"beam": (1, "beam"),
"bit": (1, "bit"),
"bin": (1, "bin"),
"byte": (1, "byte"), # I don't think we ever want to unify bit and byte
"C": (1, "C"), # A.s
"cd": (1, "cd"),
"ct": (1, "ct"),
"count": (1, "ct"),
"chan": (1, "chan"),
"D": (1e-19/3., "D"), # C.m
"d": (3600*24, "s"),
"deg": (math.pi/180., "rad"),
"erg": (1e-7, "J"),
"eV": (1.6021766208e-19, "J"),
"F": (1, "F"), # C/V
"g": (1e-3, "kg"),
"G": (1e-4, "T"),
"h": (3600., "s"),
"H": (1, "H"), # Wb/A
"Hz": (1, "Hz"), # s-1
"J": (1, "J"), # kg m^2/s^2
"Jy": (1, "Jy"), # 1e-26 W/m2/Hz
"K": (1, "K"),
"lm": (1, "lm"),
"lx": (1, "lx"),
"lyr": (2627980686828.0, "m"),
"m": (1, "m"),
"mag": (1, "mag"), # correlate that with, erm, lux?
"mas": (math.pi/180./3.6e6, "rad"),
"min": (60, "s"),
"mol": (1, "mol"),
"N": (1, "N"), # kg.m/s2
"Ohm": (1, "Ohm"), # V/A
"Pa": (1, "Pa"), # N/m^2
"pc": (3.0857e16, "m"),
"ph": (1, "ph"),
"photon": (1, "ph"),
"pix": (1, "pix"),
"pixel": (1, "pix"),
"R": (1, "R"), # Rayleigh
"rad": (1, "rad"),
"Ry": (2.17989e-18, "J"),
"s": (1, "s"),
"S": (1, "S"), # A/V
"solLum": (3.826e26, "W"),
"solMass": (1.989e30, "kg"),
"solRad": (6.9559e8, "m"),
"sr": (1, "sr"),
"T": (1, "T"), # V.s/m2
"u": (1.66053886e-27, "kg"),
"V": (1, "V"),
"voxel": (1, "voxel"),
"W": (1, "W"), # kg.m2/s3 or A.V -- that's going to be a tough one
"Wb": (1, "Wb"),
"yr": (3600*24*365.25, "s"), # This is the SI/julian year!
}
EXCEPTIONAL_CONVERSIONS = {
# tuple-form of units: conversions that are still linear but
# somehow idiosyncratic.
# wave number to energy
((('m', -1.0),), (('J', 1),)): PLANCK_H*LIGHT_C,
}
# These are the keys from PLAIN_UNITS that cannot take SI prefixes
NON_PREFIXABLE = frozenset([
"AU", "au", "D", "Ry", "arcmin", "beam", "bin", "chan", "d", "h", "mas",
"min", "ph", "photon", "pix", "pixel", "solLum", "solMass", "solRad",
"voxel"])
PREFIXES = prefixes = {"d": 1e-1, "c": 1e-2, "m":1e-3, "u":1e-6,
"n":1e-9, "p":1e-12, "f":1e-15, "a":1e-18, "z":1e-21, "y":1e-24,
"da": 10., "h":1e2, "k":1e3, "M":1e6, "G":1e9, "T":1e12, "P":1e15,
"E":1e18, "Z":1e21, "Y":1e24, "": 1}
# these map VOUnit Table 5 "magic" units to UnitNode constructor arguments
AMBIGUOUS_STRINGS = {
'Pa': ('Pa', ''),
'ha': ('yr', 'h'),
'cd': ('cd', ''),
'dB': ('Bel', 'd'),
'B': ('Byte', ''),
'au': ('AU', ''),
'dadu': ('adu', 'd'),
}
class _Node(object):
"""the abstract base for a node in a Unit tree.
All these must implement fromToks methods that can be called
as pyparsing callbacks.
"""
@classmethod
def fromToks(cls, s, p, toks):
raise NotImplementedError("%s objects have no implementation of"
" fromToks (i.e., they are broken)."%cls.__name__)
[docs]class UnitNode(_Node):
"""a terminal node containing a unit, possibly with a prefix
This is normally constructed through the fromToks constructor.
Check out the unit, prefix, and isUnknown attributes. The prefix
can be turned to a factor using the PREFIXES dictionary. An empty
prefix (factor 1) is represented by an empty string.
There is also isValid; that's isUnknown except if an unknown unit is
quoted.
"""
def __init__(self, unit, prefix=""):
self.unit = unit
self.isUnknown = self.unit not in PLAIN_UNITS
self.isValid = not self.isUnknown
self.prefix = prefix
[docs] @classmethod
def fromToks(cls, s, p, toks):
assert len(toks)==1
unitString = toks[0]
if unitString in AMBIGUOUS_STRINGS:
return cls(*AMBIGUOUS_STRINGS[unitString])
elif unitString in PLAIN_UNITS:
return cls(unitString)
elif len(unitString)>2 and unitString[:2]=="da":
return cls(unitString[2:], "da")
elif len(unitString)>1 and unitString[0] in PREFIXES:
if unitString[1:] in NON_PREFIXABLE:
raise BadUnit("No Prefixes allowed on %s"%unitString[1:])
return cls(unitString[1:], unitString[0])
else:
return cls(unitString, '')
def __eq__(self, other):
return (isinstance(other, UnitNode)
and self.unit==other.unit
and self.prefix==other.prefix)
def __str__(self):
return "%s%s"%(self.prefix, self.unit)
def __repr__(self):
if self.isUnknown:
return "U?(%s ; %s)"%(self.prefix, repr(self.unit))
else:
return "U(%s ; %s)"%(self.prefix, repr(self.unit))
[docs] def getSI(self):
"""returns a pair of factor and basic unit.
Basic units are what's in the defining pairs in the PLAIN_UNITS dict.
"""
if self.isUnknown:
# these don't hurt if they're on both sides of an equation
return PREFIXES[self.prefix], {self.unit: 1}
factor, basic = PLAIN_UNITS[self.unit]
return PREFIXES[self.prefix]*factor, {basic: 1}
[docs]class FunctionApplication(_Node):
"""A function applied to a term.
"""
_pythonFunc = {
"ln": math.log,
"log": math.log10,
"exp": math.exp,
"sqrt": math.sqrt,
}
def __init__(self, funcName, term):
self.funcName = funcName
self.term = term
self.isUnknown = (self.funcName not in self._pythonFunc
) or self.term.isUnknown
self.isValid = (self.funcName in self._pythonFunc
) and self.term.isValid
[docs] @classmethod
def fromToks(cls, s, p, toks):
assert len(toks)==2
return cls(toks[0], toks[1])
def __str__(self):
return "%s(%s)"%(self.funcName, self.term)
def __repr__(self):
return "A(%s ; %s)"%(repr(self.funcName), repr(self.term))
[docs] def getFactorWithUnit(self, unit):
"""see this method in Term.
We don't simplify any units here, as that's a bit hard in
the current logic -- though there's something to be won, in particular
with sqrt. Ah well.
"""
return None, None
[docs] def getSI(self):
factor, powers = self.term.getSI()
if self.funcName=="sqrt":
powers = dict((key, value/2.)
for key, value in powers.items())
else:
powers = dict(((self.funcName, key), value)
for key, value in powers.items())
return self._pythonFunc[self.funcName](factor), powers
[docs]class QuotedUnitNode(_Node):
"""a quoted ("defined unknown") unit.
"""
def __init__(self, unit):
self.unit = unit
self.isUnknown = True
self.isValid = True
[docs] @classmethod
def fromToks(cls, s, p, toks):
return cls(toks[1])
def __str__(self):
return "'%s'"%(self.unit)
def __repr__(self):
return "U?('%s')"%(repr(self.unit))
[docs] def getSI(self):
# These units stand for themselves. They don't really hurt if
# they're in both terms of a conversion
return 1, {self.unit: 1}
[docs]class Factor(_Node):
"""A UnitNode with a power.
"""
def __init__(self, unit, power):
self.unit, self.power = unit, power
self.isUnknown = self.unit.isUnknown
self.isValid = self.unit.isValid
[docs] @classmethod
def fromToks(cls, s, p, toks):
if len(toks)==2:
return cls(toks[0], toks[1])
elif len(toks)==1:
return cls(toks[0], 1)
else:
raise Exception("This cannot happen")
def __str__(self):
if self.power==0:
return "1"
powerLit = repr(self.power).rstrip("0").rstrip(".")
if "." in powerLit:
# see if we can come up with a nice fraction
for denom in range(2, 8):
if abs(int(self.power*denom)-self.power*denom)<1e-13:
powerLit = "(%d/%d)"%(round(self.power*denom), denom)
break
if powerLit=="1":
powerLit = ""
else:
powerLit = "**"+powerLit
return "%s%s"%(self.unit, powerLit)
def __repr__(self):
return "F(%s ; %s)"%(repr(self.unit), repr(self.power))
[docs] def getFactorWithUnit(self, unit):
"""see Term.getFactorWithUnit.
"""
if self.unit==unit:
return self, 1
else:
return None, None
[docs] def getSI(self):
factor, powers = self.unit.getSI()
powers[list(powers.keys())[0]] = self.power
return factor**self.power, powers
[docs]class Term(_Node):
"""A Node containing two factors and an operator.
The operator here is either . (product) or / (division).
"""
_expoMapping = {
'.': 1,
'/': -1,
}
def __init__(self, op1, operator, op2):
assert operator in {'.', '/'}
self.op1, self.operator, self.op2 = op1, operator, op2
self.isUnknown = self.op1.isUnknown or self.op2.isUnknown
self.isValid = self.op1.isValid and self.op2.isValid
[docs] @classmethod
def fromToks(cls, s, p, toks):
if len(toks)==1:
# avoid to many internal nodes: a 1-operand expression is the operand
return toks[0]
op1, operator, op2 = toks
if isinstance(op2, Factor):
# for now we only simplify with simple factors in the second
# operand. If there are complete expressions there, I'd
# have to do a lot more complex operations (e.g., take out
# factors from there).
existingFactorWithUnit, ef = op1.getFactorWithUnit(op2.unit)
if existingFactorWithUnit:
existingFactorWithUnit.power += \
op2.power*cls._expoMapping[operator]*ef
return op1
return cls(op1, operator, op2)
def __str__(self):
op1Lit, op2Lit = str(self.op1), str(self.op2)
if self.operator=='/' and isinstance(self.op2, Term):
op2Lit = "(%s)"%op2Lit
return "%s%s%s"%(op1Lit, self.operator, op2Lit)
def __repr__(self):
return "T(%s ; %s ; %s)"%(repr(self.op1),
repr(self.operator), repr(self.op2))
[docs] def getFactorWithUnit(self, unit):
"""returns a factor within self with unit if it exists.
Actually, it will return that factor and the exponentiation factor,
which is -1 if that factor is in the denominator, 1 otherwise.
This returns None, None if no such factor exists.
"""
from1, expo1 = self.op1.getFactorWithUnit(unit)
from2, expo2 = self.op2.getFactorWithUnit(unit)
if from1:
return from1, expo1
elif from2:
return from2, -1*expo2 if self.operator=='/' else 1
else:
return None, None
[docs] def getSI(self):
factor1, powers1 = self.op1.getSI()
factor2, powers2 = self.op2.getSI()
newPowers = powers1
if self.operator==".":
for si, power in powers2.items():
newPowers[si] = newPowers.get(si, 0)+power
return factor1*factor2, newPowers
elif self.operator=="/":
for si, power in powers2.items():
newPowers[si] = newPowers.get(si, 0)-power
return factor1/factor2, newPowers
else:
raise Exception("This can't happen")
[docs]class Expression(_Node):
"""The root node of an expression tree.
This contains a term and optionally a scale factor.
"""
def __init__(self, term, scaleFactor):
self.term, self.scaleFactor = term, scaleFactor
self.isUnknown = self.term.isUnknown
self.isValid = self.term.isValid
[docs] @classmethod
def fromToks(cls, s, p, toks):
if len(toks)==2:
return cls(toks[1], float(toks[0]))
elif len(toks)==1:
return cls(toks[0], 1)
else:
raise Exception("This can't happen")
def __str__(self):
if self.scaleFactor==1:
return str(self.term)
else:
return "%s %s"%(formatScaleFactor(self.scaleFactor), self.term)
def __repr__(self):
return "R(%s ; %s)"%(repr(self.scaleFactor), repr(self.term))
[docs] def getSI(self):
"""returns a pair of a numeric factor and a dict mapping SI units to
their powers.
"""
factor, siPowers = self.term.getSI()
return factor*self.scaleFactor, siPowers
def _buildTerm(s, pos, toks):
"""a parseAction for terms, making trees out of parse lists
left-associatively.
This will also normalise products of the same units to powers.
"""
toks = list(toks)
curOperand = toks.pop(0)
while len(toks)>1:
curOperand = Term.fromToks(s, pos, [curOperand, toks[0], toks[1]])
del toks[:2]
return curOperand
[docs]def evalAll(s, p, toks):
"""a parse action evaluating the whole match as a python expression.
Obviously, this should only be applied to carefully screened nonterminals.
"""
return eval("".join(str(tok) for tok in toks))
[docs]@functools.cache
def getUnitGrammar(debugging=False):
"""the grammar to parse VOUnits.
After initialization, the class has a "symbols" dictionary containing
the individual nonterminals.
"""
from gavo.utils.parsetricks import (Word, Literal, Regex,
Optional, ZeroOrMore, alphas,
Suppress, Forward, pyparsingWhitechars, ParserElement)
with pyparsingWhitechars(''):
unit_atom = Word(alphas).addParseAction(UnitNode.fromToks)
unit_atom.setName("atomic unit")
quoted_unit_atom = ("'" + Word(alphas) + "'"
).addParseAction(QuotedUnitNode.fromToks)
quoted_unit_atom.setName("quoted atomic unit")
OPEN_P = Literal('(')
CLOSE_P = Literal(')')
SIGN = Literal('+') | Literal('-')
FUNCTION_NAME = Word(alphas)
UNSIGNED_INTEGER = Word("01234567890")
SIGNED_INTEGER = SIGN + UNSIGNED_INTEGER
FLOAT = Regex(r"[+-]?([0-9]+(\.[0-9]*)?)")
VOFLOAT = Regex(r"0.[0-9]+([eE][+-]?[0-9]+)?"
"|[1-9][0-9]*(\.[0-9]+)?([eE][+-]?[0-9]+)?")
integer = (SIGNED_INTEGER | UNSIGNED_INTEGER).setName("integer")
power_operator = Literal('**')
multiplication_operator = Literal(".")
division_operator = Literal("/")
numeric_power = (integer
| OPEN_P + integer + CLOSE_P
| OPEN_P + FLOAT + CLOSE_P
| OPEN_P + integer + '/'
+ UNSIGNED_INTEGER.addParseAction(lambda s, p, t: t[0]+".")
+ CLOSE_P)
numeric_power.setParseAction(evalAll)
pow_10 = Literal("10") + power_operator + numeric_power
scale_factor = (pow_10 | VOFLOAT
).setName("scale_factor"
).setParseAction(evalAll)
any_unit_atom = unit_atom | quoted_unit_atom
factor = (any_unit_atom
+ Optional( Suppress(power_operator) + numeric_power )
).setName("factor"
).addParseAction(Factor.fromToks)
complete_expression = Forward()
function_application = (FUNCTION_NAME +
Suppress(OPEN_P) + complete_expression + Suppress(CLOSE_P)
).setName("function application")
function_application.addParseAction(FunctionApplication.fromToks)
unit_expression = (
Suppress(OPEN_P) + complete_expression + Suppress(CLOSE_P)
| ( factor
^ function_application )).setName("unit expression")
product_of_units = (unit_expression
+ ZeroOrMore(multiplication_operator + unit_expression)
).setName("units term"
).setParseAction(_buildTerm)
complete_expression << (
product_of_units + Optional(division_operator + unit_expression) )
complete_expression.setParseAction(Term.fromToks)
input = (Optional(scale_factor)
+ complete_expression
).setName("unit expression"
).setParseAction(Expression.fromToks)
if debugging:
for name, sym in locals().items():
if isinstance(sym, ParserElement):
sym.setDebug(True)
sym.setName(name)
del debugging
getUnitGrammar.symbols = locals()
return input
[docs]def asSequence(unitDict):
"""returns a sorted tuple of (si, power) for a result of getSI().
These should be more suitable for dimensional analysis.
"""
return tuple(sorted(unitDict.items()))
[docs]def parseUnit(unitStr, unitGrammar=getUnitGrammar()):
try:
return utils.pyparseString(unitGrammar, unitStr, parseAll=True)[0]
except parsetricks.ParseException as msg:
raise utils.logOldExc(
BadUnit("%s at col. %d"%(repr(unitStr), msg.column)))
[docs]def computeConversionFactor(unitStr1, unitStr2):
"""returns the factor needed to get from quantities given in unitStr1
to unitStr2.
Both must be given in VOUnits form.
This function may raise a BadUnit if one of the strings are
malformed, or an IncompatibleUnit exception if the units don't have
the same SI base.
If the function is successful, unitStr1 = result*unitStr2
"""
if unitStr1==unitStr2:
return 1
factor1, powers1 = parseUnit(unitStr1).getSI()
factor2, powers2 = parseUnit(unitStr2).getSI()
powers1, powers2 = asSequence(powers1), asSequence(powers2)
if (powers1, powers2) in EXCEPTIONAL_CONVERSIONS:
return factor1/factor2*EXCEPTIONAL_CONVERSIONS[powers1, powers2]
if (powers2, powers1) in EXCEPTIONAL_CONVERSIONS:
return factor1/factor2/EXCEPTIONAL_CONVERSIONS[powers2, powers1]
if powers1!=powers2:
raise IncompatibleUnits("%s and %s do not have the same SI base"%(
unitStr1, unitStr2))
# tuples as keys in powers come from non-polynomial function
# applications; in such cases, multiplication is not good enough
# for conversions, and thus we give up.
for u, _ in powers1:
if isinstance(u, tuple):
raise IncompatibleUnits("%s has a non-polynomial function. No"
" conversion by multiplication possible"%(unitStr1))
for u, _ in powers2:
if isinstance(u, tuple):
raise IncompatibleUnits("%s has a non-polynomial function. No"
" conversion by multiplication possible"%(unitStr2))
return factor1/factor2
[docs]def computeColumnConversions(newColumns, oldColumns):
"""returns a dict of conversion factors between newColumns and oldColumns.
Both arguments are iterables of columns.
For every column in newColumn, the function sees if the units of
newColumn and oldColumn match. If they don't, compute a conversion
factor to be multiplied to oldColumns values to make them newColumns
values and add it to the result dict.
The function raises a DataError if a column in newColumns has no
match in oldColumns.
"""
res = {}
for newCol in newColumns:
if not newCol.name in oldColumns:
raise utils.DataError(
"Request for column %s from %s cannot be satisfied in %s"%(
newCol.name, newColumns, oldColumns))
oldCol = oldColumns.getColumnByName(newCol.name)
try:
if newCol.unit!=oldCol.unit:
res[newCol.name] = computeConversionFactor(oldCol.unit, newCol.unit)
except BadUnit: # we ignore bad units, assume they'll be handled by
# valuemappers.
pass
return res
if __name__=="__main__":
g = getUnitGrammar(debugging=True)
res = g.parseString("10**28s", parseAll=True)[0]
print(res)