source: sasview/src/sas/sascalc/data_util/nxsunit.py @ 574adc7

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 574adc7 was 574adc7, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

convert sascalc to python 2/3 syntax

  • Property mode set to 100644
File size: 8.1 KB
Line 
1# This program is public domain
2# Author: Paul Kienzle
3"""
4Define unit conversion support for NeXus style units.
5
6The unit format is somewhat complicated.  There are variant spellings
7and incorrect capitalization to worry about, as well as forms such as
8"mili*metre" and "1e-7 seconds".
9
10This is a minimal implementation of units including only what I happen to
11need now.  It does not support the complete dimensional analysis provided
12by the package udunits on which NeXus is based, or even the units used
13in the NeXus definition files.
14
15Unlike other units packages, this package does not carry the units along with
16the value but merely provides a conversion function for transforming values.
17
18Usage example::
19
20    import nxsunit
21    u = nxsunit.Converter('mili*metre')  # Units stored in mm
22    v = u(3000,'m')  # Convert the value 3000 mm into meters
23
24NeXus example::
25
26    # Load sample orientation in radians regardless of how it is stored.
27    # 1. Open the path
28    file.openpath('/entry1/sample/sample_orientation')
29    # 2. scan the attributes, retrieving 'units'
30    units = [for attr,value in file.attrs() if attr == 'units']
31    # 3. set up the converter (assumes that units actually exists)
32    u = nxsunit.Converter(units[0])
33    # 4. read the data and convert to the correct units
34    v = u(file.read(),'radians')
35
36This is a standalone module, not relying on either DANSE or NeXus, and
37can be used for other unit conversion tasks.
38
39Note: minutes are used for angle and seconds are used for time.  We
40cannot tell what the correct interpretation is without knowing something
41about the fields themselves.  If this becomes an issue, we will need to
42allow the application to set the dimension for the unit rather than
43inferring the dimension from an example unit.
44"""
45
46# TODO: Add udunits to NAPI rather than reimplementing it in python
47# TODO: Alternatively, parse the udunits database directly
48# UDUnits:
49#  http://www.unidata.ucar.edu/software/udunits/udunits-1/udunits.txt
50
51from __future__ import division
52import math
53
54__all__ = ['Converter']
55
56# Limited form of units for returning objects of a specific type.
57# Maybe want to do full units handling with e.g., pyre's
58# unit class. For now lets keep it simple.  Note that
59def _build_metric_units(unit,abbr):
60    """
61    Construct standard SI names for the given unit.
62    Builds e.g.,
63        s, ns
64        second, nanosecond, nano*second
65        seconds, nanoseconds
66    Includes prefixes for femto through peta.
67
68    Ack! Allows, e.g., Coulomb and coulomb even though Coulomb is not
69    a unit because some NeXus files store it that way!
70
71    Returns a dictionary of names and scales.
72    """
73    prefix = dict(peta=1e15,tera=1e12,giga=1e9,mega=1e6,kilo=1e3,
74                  deci=1e-1,centi=1e-2,milli=1e-3,mili=1e-3,micro=1e-6,
75                  nano=1e-9,pico=1e-12,femto=1e-15)
76    short_prefix = dict(P=1e15,T=1e12,G=1e9,M=1e6,k=1e3,
77                        d=1e-1,c=1e-2,m=1e-3,u=1e-6,
78                        n=1e-9,p=1e-12,f=1e-15)
79    map = {abbr:1}
80    map.update([(P+abbr,scale) for (P,scale) in short_prefix.items()])
81    for name in [unit,unit.capitalize()]:
82        map.update({name:1,name+'s':1})
83        map.update([(P+name,scale) for (P,scale) in prefix.items()])
84        map.update([(P+'*'+name,scale) for (P,scale) in prefix.items()])
85        map.update([(P+name+'s',scale) for (P,scale) in prefix.items()])
86    return map
87
88def _build_plural_units(**kw):
89    """
90    Construct names for the given units.  Builds singular and plural form.
91    """
92    map = {}
93    map.update([(name,scale) for name,scale in kw.items()])
94    map.update([(name+'s',scale) for name,scale in kw.items()])
95    return map
96
97def _caret_optional(s):
98    """
99    Strip '^' from unit names.
100
101    * WARNING * this will incorrect transform 10^3 to 103.
102    """
103    s.update((k.replace('^',''),v)
104             for k, v in list(s.items())
105             if '^' in k)
106
107def _build_all_units():
108    distance = _build_metric_units('meter','m')
109    distance.update(_build_metric_units('metre','m'))
110    distance.update(_build_plural_units(micron=1e-6, Angstrom=1e-10))
111    distance.update({'A':1e-10, 'Ang':1e-10})
112
113    # Note: minutes are used for angle
114    time = _build_metric_units('second','s')
115    time.update(_build_plural_units(hour=3600,day=24*3600,week=7*24*3600))
116
117    # Note: seconds are used for time
118    angle = _build_plural_units(degree=1, minute=1/60.,
119                  arcminute=1/60., arcsecond=1/3600., radian=180/math.pi)
120    angle.update(deg=1, arcmin=1/60., arcsec=1/3600., rad=180/math.pi)
121
122    frequency = _build_metric_units('hertz','Hz')
123    frequency.update(_build_metric_units('Hertz','Hz'))
124    frequency.update(_build_plural_units(rpm=1/60.))
125
126    # Note: degrees are used for angle
127    # Note: temperature needs an offset as well as a scale
128    temperature = _build_metric_units('kelvin','K')
129    temperature.update(_build_metric_units('Kelvin','K'))
130    temperature.update(_build_metric_units('Celcius', 'C'))
131    temperature.update(_build_metric_units('celcius', 'C'))
132
133    charge = _build_metric_units('coulomb','C')
134    charge.update({'microAmp*hour':0.0036})
135
136    sld = { '10^-6 Angstrom^-2': 1e-6, 'Angstrom^-2': 1 }
137    Q = { 'invA': 1, 'invAng': 1, 'invAngstroms': 1, '1/A': 1,
138          '10^-3 Angstrom^-1': 1e-3, '1/cm': 1e-8, '1/m': 1e-10,
139          'nm^-1': 0.1, '1/nm': 0.1, 'n_m^-1': 0.1 }
140
141    _caret_optional(sld)
142    _caret_optional(Q)
143
144    dims = [distance, time, angle, frequency, temperature, charge, sld, Q]
145    return dims
146
147class Converter(object):
148    """
149    Unit converter for NeXus style units.
150    """
151    # Define the units, using both American and European spelling.
152    scalemap = None
153    scalebase = 1
154    dims = _build_all_units()
155
156    # Note: a.u. stands for arbitrary units, which should return the default
157    # units for that particular dimension.
158    # Note: don't have support for dimensionless units.
159    unknown = {None:1, '???':1, '': 1, 'a.u.': 1}
160
161    def __init__(self, name):
162        self.base = name
163        for map in self.dims:
164            if name in map:
165                self.scalemap = map
166                self.scalebase = self.scalemap[name]
167                return
168        if name in self.unknown:
169            return # default scalemap and scalebase correspond to unknown
170        else:
171            raise KeyError("Unknown unit %s"%name)
172
173    def scale(self, units=""):
174        if units == "" or self.scalemap is None: return 1
175        return self.scalebase/self.scalemap[units]
176
177    def __call__(self, value, units=""):
178        # Note: calculating a*1 rather than simply returning a would produce
179        # an unnecessary copy of the array, which in the case of the raw
180        # counts array would be bad.  Sometimes copying and other times
181        # not copying is also bad, but copy on modify semantics isn't
182        # supported.
183        if units == "" or self.scalemap is None: return value
184        try:
185            return value * (self.scalebase/self.scalemap[units])
186        except KeyError:
187            possible_units = ", ".join(str(k) for k in self.scalemap.keys())
188            raise KeyError("%s not in %s"%(units,possible_units))
189
190def _check(expect,get):
191    if expect != get:
192        raise ValueError("Expected %s but got %s"%(expect, get))
193     #print expect,"==",get
194
195def test():
196    _check(1,Converter('n_m^-1')(10,'invA')) # 10 nm^-1 = 1 inv Angstroms
197    _check(2,Converter('mm')(2000,'m')) # 2000 mm -> 2 m
198    _check(2.011e10,Converter('1/A')(2.011,"1/m")) # 2.011 1/A -> 2.011 * 10^10 1/m
199    _check(0.003,Converter('microseconds')(3,units='ms')) # 3 us -> 0.003 ms
200    _check(45,Converter('nanokelvin')(45))  # 45 nK -> 45 nK
201    _check(0.5,Converter('seconds')(1800,units='hours')) # 1800 s -> 0.5 hr
202    _check(123,Converter('a.u.')(123,units='mm')) # arbitrary units always returns the same value
203    _check(123,Converter('a.u.')(123,units='s')) # arbitrary units always returns the same value
204    _check(123,Converter('a.u.')(123,units='')) # arbitrary units always returns the same value
205    try:
206        Converter('help')
207    except KeyError:
208        pass
209    else:
210        raise Exception("unknown unit did not raise an error")
211
212    # TODO: more tests
213
214if __name__ == "__main__":
215    test()
Note: See TracBrowser for help on using the repository browser.