[SciPy-user] Units in SciPy

[Robert Kern]

Duncan Child wrote:

> 3) does anyone have input or suggestions regarding the temperature issue?

I've thought about this a little more, and I would like to describe here 
an outline of the system that I would like to see with examples drawn 
from my current research.

One of the many things that Konrad got right was the name: the system is 
dealing with physical quantities. The most general form would consist of 
"numerical quantities tagged by their measurement system." A large 
subset can be described as "numerical quantities with units." What we 
need is a system that can perform operations and conversions between 
measurement systems, some of which will be user-defined and arbitrarily 
complicated.

The methods for dealing with the unit-subset are sufficiently covered by 
the currently available software. With the exception of certain, very 
standard conversions like temperature scales, most of the 
non-unit-subset must be user-defined.

Allow me to describe some of the many measurement systems in my current 
project in more detail than you probably care; feel free to skip ahead. 
I am a geophysicist. I am modelling the distribution of slip over a 
fault surface during and after an earthquake. I am using measurements 
from a local GPS network and Interferometric Synthetic Aperture Radar 
(InSAR: I get images of the very small displacements from each pixel 
between two flyovers of a sattelite; it's kind of like having a very 
dense GPS network that doesn't get sampled very often). My measurement 
systems include:

Location:
  * Latitude and longitude of various entities: GPS stations, corners of 
fault patches, the corners of my InSAR images
  * Local northings, eastings, up (essentially the meters north, east, 
and up from a given local reference point)
  * X, Y, Z referenced to the fault geometry
  * Pixels in the InSAR images
  * Distance along-slip and along-dip on the fault (a location on the 
plane of the fault place itself)

Displacements:
  * Differential ground displacements of the GPS given a reference station
  * Differential ground displacements from the InSAR images given a 
reference pixel projected onto the line of sight of the satellite
  * Absolute ground displacements given a slip model

Time:
  * GPS time: seconds from an arbitrary start point, no leap seconds
  * Satellite time: seconds from Jan 1 00:00 of the launch year of the 
satellite, no leap seconds
  * Calendrical time: leap seconds!

Miscellaneous:
  * Slip, rake: magnitude and direction of the slip for each fault patch
  * Dip-slip, strike-slip: orthogonal slip components for each fault patch
  * Green functions: observable displacements (GPS and InSAR) as a 
function of slip across each fault patch
  * Smoothing kernel: for regularizing my inverse problem

... and probably some more. Now, I don't claim that even given my dream 
physical quantities system that I would commit all of these entities to 
it. Some, like time, are more suited to the strategy of converting 
everything as soon as it enters the picture. Others don't need 
conversions between different measurement systems.

Skip to here if you are bored by geophysics.

I have found that I have been duplicating snippets of code, throwing 
around little one-to-one conversion functions, or worst, just not 
bothering with some task because it's too annoying to switch between 
these systems.

Here is my sketch of a physical quantities system:

Each quantity has a value, which can be anything reasonably number-like, 
and measurement system, which probably should be a full-fledged object, 
but perhaps could also be a string for standard shortcuts.

   In[1]: x = quantity(10.5, apples)
   In[2]: y = quantity(12, oranges)

There will exist a function for converting quantities from one system to 
another.

   In[3]: z = convert(x, oranges)

In the case of real, honest-to-goodness units like "m" and "degC" (the 
differential kind, not absolute systems, not even Kelvin or Rankine), 
converting falls down to the well-known algorithms.

Failing that, convert() will look for an appropriate conversion path 
from all of the measurement system objects registered with it. This 
lookup algorithm can be stolen from PyProtocols[1]. Operations that can 
be performed on quantities of a given measurement system will be defined 
by the measurement system object.

For example, quantities in the (lat, lon) system can be converted to 
(northings, eastings). To a (lat, lon) quantity, I can add or subtract a 
unitted quantity of dimensions (length, length). I cannot add another 
(lat, lon) quantity, but I can subtract it yielding a (length, length), 
which, given a suitable context, can be converted to (range, azimuth). 
Multiplication does not make sense, nor does division, but trigonometric 
functions might depending on your use cases.

Of course, temperature scales fit neatly into this system. Since the 
"shift-scale-shift" type of conversion happens reasonably often, a bit 
of general-purpose scaffolding can be constructed for it and which the 
temperature scales can use.

[1] http://peak.telecommunity.com/PyProtocols.html

-- 
Robert Kern
rkern at ucsd.edu

"In the fields of hell where the grass grows high
  Are the graves of dreams allowed to die."
   -- Richard Harter