Report uncertainty: given a mean and the standard error, show only significant figures

The intent is to show the result of several observations without unnecessary digits i.e., to display a value with the number of significant digits that is consistent with a given uncertainty.

For example, if computed mean=123.45 and err=0.0012345 then the expected output may look like 123450 ± 1.2 (× 10^-3) where the following rules are used:

the error always has one or two significant digits. Two if the first digit is 1 (ignoring the leading zeros)
the mean value is rounded to drop uncertain digits except for the last one ("stop the mean at the same decade as that of the first significant (non-zero) digit in the SEM"). Trailing zeros are added to display the precision corresponding to the error if necessary.

How it might be used in Python:

import statistics

mean = statistics.mean(measurements)
err = statistics.stdev(measurements, mean) / len(measurements) ** 0.5
print("{} ± {} (×10<sup>{}</sup>)".format(*round_to_uncertainty(mean, err)))

The question is how to implement the round_to_uncertainty(value, uncertainty) function expressing the rules 1 and 2 above.

Note: the terms error, uncertainty are used loosely in the question. See the Guide to the Expression of Uncertainty in Measurement (GUM). Here's a related question for R.

Solution

decimal module could be used to manipulate the decimal representation of the numbers conveniently:

from decimal import Decimal

def round_to_uncertainty(value, uncertainty):
    # round the uncertainty to 1-2 significant digits
    u = Decimal(uncertainty).normalize()
    exponent = u.adjusted()  # find position of the most significant digit
    precision = (u.as_tuple().digits[0] == 1)  # is the first digit 1?
    u = u.scaleb(-exponent).quantize(Decimal(10)**-precision)

    # round the value to remove excess digits
    return round(Decimal(value).scaleb(-exponent).quantize(u)), u, exponent

Example:

for mean, err in [
    (123.45, 0.0012345),    # 123450 ± 1.2 (×10<sup>-3</sup>)
    (8165.666, 338.9741),   # 82 ± 3 (×10<sup>2</sup>)
]: 
    print("{} ± {} (×10<sup>{}</sup>)".format(*round_to_uncertainty(mean, err)))

The input 123.45, 0.0012345 is reported as 123450 ± 1.2 (×10^-3). And 8165.666, 338.9741 translates to 82 ± 3 (×10²) according to the rules from the current question.