Repository: IAMconsortium/pyam
Pull request #541 (2021-05-31 17:01:05)

Support units in binary operations
# Please confirm that this PR has done the following:

- [x] Tests Added
- [x] Documentation Added
- ~Name of contributors Added to AUTHORS.rst~
- [x] Description in RELEASE_NOTES.md Added

# Description of PR

This PR changes the implementation of the binary ops to support units (via iam-units and pint).

Default behaviour:
 - Cast values and units to a `pint.Quantity` using the **iam-units** registry and perform operations with pint unit handling (`ignore_units=False`)
    **Warning**: By casting units to pints and then performing operations, the return-format may look different.
   In some of the tests, `EJ/yr` is transformed to `exajoule / year` (in full format) or `EJ / a` (in compact format, implemented here). This is over-ridden for some methods (addition, multiplication, division) if all units are identical.
 - Passing a pint quantity as argument will work, e.g, `df.add("Primary Energy", Quantity(2, "EJ/yr"), "new")`.
 - **pint** will raise appropriate errors if it doesn't know how to handle the units.
 - In line with common practice in the IAM community, a dimensionless quantity is marked as unit `""` (empty string) rather than the word `dimensionless` (this works in a round-trip to-from-to pint).

Alternative over-ride:
 - A user can set `ignore_units=True`, in which case the operations will work just on the values and set `unit` to "unknown" (per suggestion by @znicholls in https://github.com/IAMconsortium/pyam/issues/536#issuecomment-847569955) or to the value of `ignore_units` (if it's a string).

Known issues:
- pandas fillna does not work with non-numeric values (in particular not with pint.Quantity).

Side note:
 - I simplified the test for the `apply()` function.

To-dos in later PRs (it's already quite a beast):
 - Make the fillna feature work with pint.
 - Add `ignore_units` to the `apply()` function.
 - Support custom pint registry instances.
 - Pass unit formatting options to pint as kwargs.

closes #537 
closes #535 


codecov[bot]: # [Codecov](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) Report
> Merging [#541](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) (7bb1333) into [main](https://codecov.io/gh/IAMconsortium/pyam/commit/a4c6c71a20666ac70a1812f77b61d80e18cfd96e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) (a4c6c71) will **decrease** coverage by `0.2%`.
> The diff coverage is `100.0%`.

[![Impacted file tree graph](https://codecov.io/gh/IAMconsortium/pyam/pull/541/graphs/tree.svg?width=650&height=150&src=pr&token=ThxWlNkonJ&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium)](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium)

```diff
@@           Coverage Diff           @@
##            main    #541     +/-   ##
=======================================
- Coverage   93.5%   93.2%   -0.3%     
=======================================
  Files         47      47             
  Lines       5167    5173      +6     
=======================================
- Hits        4832    4826      -6     
- Misses       335     347     +12     
```


| [Impacted Files](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium) | Coverage Δ | |
|---|---|---|
| [pyam/\_ops.py](https://codecov.io/gh/IAMconsortium/pyam/pull/541/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium#diff-cHlhbS9fb3BzLnB5) | `81.1% <100.0%> (-14.9%)` | :arrow_down: |
| [pyam/units.py](https://codecov.io/gh/IAMconsortium/pyam/pull/541/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium#diff-cHlhbS91bml0cy5weQ==) | `92.9% <100.0%> (ø)` | |
| [tests/test\_ops.py](https://codecov.io/gh/IAMconsortium/pyam/pull/541/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium#diff-dGVzdHMvdGVzdF9vcHMucHk=) | `96.7% <100.0%> (-3.3%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium). Last update [a4c6c71...7bb1333](https://codecov.io/gh/IAMconsortium/pyam/pull/541?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=IAMconsortium).

pjuergens: So it's just a comment and not a full review. I hope some of the others have time for an in-depth-review, cause I'd really like to use this feature :)
danielhuppmann: Thanks for the suggestion, @pjuergens - implemented it and updated the unit tests, because shorter units are generally preferable, I think. Also updated the description of this PR to reflect this change.

However, note that the underlying issue remains, because a user will generally not receive the same notation for a unit - in our case, `EJ / yr` is now returned as `EJ / a`. I don't see that as a problem, just something to be aware of...
danielhuppmann: Thanks for the comment @znicholls! I had a look at [pint-pandas](https://github.com/hgrecco/pint-pandas), but as you write, this assumes that an entire column in a dataframe has the same unit - this is not something that pyam satisfies, not even on a per-variable level. A user could currently have a scenario ensemble where "Primary Energy" is in EJ for some scenarios and "GWh" in others. Something to discuss for a pyam v2.0 release...

About speed and casting, I'd follow the strategy I was taught by @gidden - first make it work, then make it fast...
znicholls: > About speed and casting, I'd follow the strategy I was taught by @gidden - first make it work, then make it fast...

It's a great strategy!

Patch:
diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md
index 01bc110..0966bb1 100644
--- a/RELEASE_NOTES.md
+++ b/RELEASE_NOTES.md
@@ -1,5 +1,6 @@
 # Next Release
 
+- [#541](https://github.com/IAMconsortium/pyam/pull/541) Support units in binary operations
 - [#538](https://github.com/IAMconsortium/pyam/pull/538) Add option to set defaults in binary operations
 - [#537](https://github.com/IAMconsortium/pyam/pull/537) Enhance binary ops to support numerical arguments
 - [#533](https://github.com/IAMconsortium/pyam/pull/533) Add an `apply()` function for custom mathematical operations
diff --git a/pyam/_ops.py b/pyam/_ops.py
index cbeb2d0..4aff3b2 100644
--- a/pyam/_ops.py
+++ b/pyam/_ops.py
@@ -1,24 +1,39 @@
 import operator
 import pandas as pd
-from pyam.index import append_index_level, get_index_levels
+from pyam.index import append_index_level, get_index_levels, replace_index_values
 from pyam.utils import to_list
+from iam_units import registry
+from pint import Quantity
 
 
 # these functions have to be defined explicitly to allow calling them with keyword args
 def add(a, b):
-    return operator.add(a, b)
+    return operator.add(*_make_series(a, b))
 
 
 def subtract(a, b):
-    return operator.sub(a, b)
+    return operator.sub(*_make_series(a, b))
 
 
 def multiply(a, b):
-    return operator.mul(a, b)
+    return operator.mul(*_make_series(a, b))
 
 
 def divide(a, b):
-    return operator.truediv(a, b)
+    return operator.truediv(*_make_series(a, b))
+
+
+def _make_series(a, b):
+    """Cast instances of a pint.Quantity to a pd.Series
+
+    Calling an operation on a pd.Series and a pint.Quantity removes the units,
+    therefore a pint.Quantity is transformed to a pd.Series of quantities.
+    """
+    if isinstance(a, Quantity):
+        return pd.Series([a] * len(b), index=b.index), b
+    if isinstance(b, Quantity):
+        return a, pd.Series([b] * len(a), index=a.index)
+    return a, b
 
 
 KNOWN_OPS = {
@@ -29,7 +44,7 @@ KNOWN_OPS = {
 }
 
 
-def _op_data(df, name, method, axis, fillna=None, args=(), **kwds):
+def _op_data(df, name, method, axis, fillna=None, args=(), ignore_units=False, **kwds):
     """Internal implementation of numerical operations on timeseries"""
 
     if axis not in df._data.index.names:
@@ -46,13 +61,39 @@ def _op_data(df, name, method, axis, fillna=None, args=(), **kwds):
 
     # replace args and and kwds with values of `df._data` if applicable
     # _data_args and _data_kwds track if an argument was replaced by `df._data` values
-    _args, _data_args = [None] * len(args), [False] * len(args)
+    n = len(args)
+    _args, _data_args, _units_args = [None] * n, [False] * n, [None] * n
     for i, value in enumerate(args):
-        _args[i], _data_args[i] = _get_values(df, axis, value, cols, f"_arg{i}")
+        _args[i], _units_args[i], _data_args[i] = _get_values(
+            df, axis, value, cols, f"_arg{i}"
+        )
 
-    _data_kwds = {}
+    _data_kwds, _unit_kwds = {}, {}
     for i, (key, value) in enumerate(kwds.items()):
-        kwds[key], _data_kwds[key] = _get_values(df, axis, value, cols, key)
+        kwds[key], _unit_kwds[key], _data_kwds[key] = _get_values(
+            df, axis, value, cols, key
+        )
+
+    # fast-pass on units: override pint for some methods if all kwds have the same unit
+    if (
+        method in [add, subtract, divide]
+        and ignore_units is False
+        and fillna is None
+        and len(_unit_kwds["a"]) == 1
+        and len(_unit_kwds["b"]) == 1
+        and registry.Unit(_unit_kwds["a"][0]) == registry.Unit(_unit_kwds["b"][0])
+    ):
+        # activate ignore-units feature
+        ignore_units = _unit_kwds["a"][0] if method in [add, subtract] else ""
+        # downcast `pint.Quantity` to numerical value
+        kwds["a"], kwds["b"] = _to_value(kwds["a"]), _to_value(kwds["b"])
+
+    # cast args and kwds to pd.Series of pint.Quantity
+    if ignore_units is False:
+        for i, is_data in enumerate(_data_args):
+            _args[i] = _to_quantity(_args[i]) if is_data else _args[i]
+        for key, value in kwds.items():
+            kwds[key] = _to_quantity(value) if _data_kwds[key] else value
 
     # merge all args and kwds that are based on `df._data` to apply fillna
     if fillna:
@@ -70,14 +111,33 @@ def _op_data(df, name, method, axis, fillna=None, args=(), **kwds):
                 kwds[key] = _data[key]
 
     # apply method and check that returned object is valid
-    _value = method(*_args, **kwds)
-    if not isinstance(_value, pd.Series):
+    result = method(*_args, **kwds)
+    if not isinstance(result, pd.Series):
         msg = f"Value returned by `{method.__name__}` cannot be cast to an IamDataFrame"
-        raise ValueError(f"{msg}: {_value}")
-
-    # insert the index level and reset the series name
+        raise ValueError(f"{msg}: {result}")
+
+    # separate pint quantities into numerical value and unit (as index)
+    if ignore_units is False:
+        _value = pd.DataFrame(
+            [[i.magnitude, "{:~}".format(i.units)] for i in result.values],
+            columns=["value", "unit"],
+            index=result.index,
+        ).set_index("unit", append=True)
+        _value.index = replace_index_values(_value, "unit", {"dimensionless": ""})
+
+    # otherwise, set unit (as index) to "unknown" or the value given by "ignore_units"
+    else:
+        index = append_index_level(
+            result.reset_index("unit", drop=True).index,
+            codes=0,
+            level="unknown" if ignore_units is True else ignore_units,
+            name="unit",
+        )
+        _value = pd.Series(result.values, index=index, name="value")
+
+    # append the `name` to the index on the `axis`
     _value.index = append_index_level(_value.index, codes=0, level=name, name=axis)
-    return _value.rename(index=None)
+    return _value
 
 
 def _get_values(df, axis, value, cols, name):
@@ -98,11 +158,37 @@ def _get_values(df, axis, value, cols, name):
 
     Returns
     -------
-    Either filtered timeseries from `df` or `value`
+    Tuple of the following:
+     - Either `df.data` downselected by `{axis: value}` or `value`
+     - List of units of the timeseries data or `value`
+     - Bool whether first item was derived from `df.data`
 
     """
+    # check if `value` is a `pint.Quantity` and return unit specifically
+    if isinstance(value, Quantity):
+        return value, [value.units], False
+    # try selecting from `df.data`
     if any(v in get_index_levels(df._data, axis) for v in to_list(value)):
-        _data = df.filter(**{axis: value})._data.groupby(cols).sum().rename(index=name)
-        return _data, True
-    else:
-        return value, False
+        _df = df.filter(**{axis: value})
+        return _df._data.groupby(cols).sum().rename(index=name), _df.unit, True
+    # else, return value
+    return value, [], False
+
+
+def _to_quantity(data):
+    """Convert the values of an indexed pd.Series into pint.Quantity instances"""
+    return pd.Series(
+        [
+            registry.Quantity(v, u)
+            for v, u in zip(data.values, data.index.get_level_values("unit"))
+        ],
+        index=data.reset_index("unit", drop=True).index,
+        name=data.name,
+    )
+
+
+def _to_value(x):
+    """Return the value of a pint.Quantity"""
+    if isinstance(x, Quantity):
+        return x.magnitude
+    return x
diff --git a/pyam/core.py b/pyam/core.py
index c322cbe..e1743e4 100755
--- a/pyam/core.py
+++ b/pyam/core.py
@@ -1233,6 +1233,11 @@ class IamDataFrame(object):
         :class:`IamDataFrame` or **None**
             Aggregated timeseries data or None if `append=True`.
 
+        See Also
+        --------
+        add : Add timeseries data items along an `axis`.
+        aggregate_region : Aggregate timeseries data along the `region` dimension.
+
         Notes
         -----
         The aggregation function interprets any missing values (:any:`numpy.nan`)
@@ -1351,6 +1356,16 @@ class IamDataFrame(object):
         append : bool, default False
             append the aggregate timeseries to `self` and return None,
             else return aggregate timeseries as new :class:`IamDataFrame`
+
+        Returns
+        -------
+        :class:`IamDataFrame` or **None**
+            Aggregated timeseries data or None if `append=True`.
+
+        See Also
+        --------
+        add : Add timeseries data items `a` and `b` along an `axis`
+        aggregate : Aggregate timeseries data along the `variable` hierarchy.
         """
         _df = _aggregate_region(
             self,
@@ -1778,8 +1793,10 @@ class IamDataFrame(object):
         else:
             self.meta[col] = self.meta[col].apply(func, *args, **kwargs)
 
-    def add(self, a, b, name, axis="variable", fillna=None, append=False):
-        """Compute the addition of timeseries data between `a` and `b` along an `axis`
+    def add(
+        self, a, b, name, axis="variable", fillna=None, ignore_units=False, append=False
+    ):
+        """Add timeseries data items `a` and `b` along an `axis`
 
         This function computes `a + b`. If `a` or `b` are lists, the method applies
         :meth:`pandas.groupby().sum() <pandas.core.groupby.GroupBy.sum>` on each group.
@@ -1797,6 +1814,9 @@ class IamDataFrame(object):
         fillna : dict or scalar, optional
             Value to fill holes when rows are not defined for either `a` or `b`.
             Can be a scalar or a dictionary of the form :code:`{arg: default}`.
+        ignore_units : bool or str, optional
+            Perform operation on values without considering units. Set units of returned
+            data to `unknown` (if True) or the value of `ignore_units` (if str).
         append : bool, optional
             Whether to append aggregated timeseries data to this instance.
 
@@ -1804,15 +1824,37 @@ class IamDataFrame(object):
         -------
         :class:`IamDataFrame` or **None**
             Computed timeseries data or None if `append=True`.
+
+        See Also
+        --------
+        subtract, multiply, divide
+        apply : Apply a custom function on the timeseries data along any axis.
+        aggregate : Aggregate timeseries data along the `variable` hierarchy.
+        aggregate_region : Aggregate timeseries data along the `region` dimension.
+
+        Notes
+        -----
+        This function uses the :mod:`pint` package and the :mod:`iam-units` registry
+        (`read the docs <https://github.com/IAMconsortium/units>`_) to handle units.
+        :mod:`pyam` will keep notation consistent with the input format (if possible)
+        and otherwise uses abbreviated units :code:`'{:~}'.format(u)` (see
+        `here <https://pint.readthedocs.io/en/stable/tutorial.html#string-formatting>`_
+        for more information).
+
+        As a result, the notation of returned units may differ from the input format.
+        For example, the unit :code:`EJ/yr` may be reformatted to :code:`EJ / a`.
         """
-        _value = _op_data(self, name, "add", axis=axis, fillna=fillna, a=a, b=b)
+        kwds = dict(axis=axis, fillna=fillna, ignore_units=ignore_units)
+        _value = _op_data(self, name, "add", **kwds, a=a, b=b)
         if append:
             self.append(_value, inplace=True)
         else:
             return IamDataFrame(_value, meta=self.meta)
 
-    def subtract(self, a, b, name, axis="variable", fillna=None, append=False):
-        """Compute the difference of timeseries data between `a` and `b` along an `axis`
+    def subtract(
+        self, a, b, name, axis="variable", fillna=None, ignore_units=False, append=False
+    ):
+        """Compute the difference of timeseries data items `a` and `b` along an `axis`
 
         This function computes `a - b`. If `a` or `b` are lists, the method applies
         :meth:`pandas.groupby().sum() <pandas.core.groupby.GroupBy.sum>` on each group.
@@ -1830,6 +1872,9 @@ class IamDataFrame(object):
         fillna : dict or scalar, optional
             Value to fill holes when rows are not defined for either `a` or `b`.
             Can be a scalar or a dictionary of the form :code:`{arg: default}`.
+        ignore_units : bool or str, optional
+            Perform operation on values without considering units. Set units of returned
+            data to `unknown` (if True) or the value of `ignore_units` (if str).
         append : bool, optional
             Whether to append aggregated timeseries data to this instance.
 
@@ -1837,15 +1882,35 @@ class IamDataFrame(object):
         -------
         :class:`IamDataFrame` or **None**
             Computed timeseries data or None if `append=True`.
+
+        See Also
+        --------
+        add, multiply, divide
+        apply : Apply a custom function on the timeseries data along any axis.
+
+        Notes
+        -----
+        This function uses the :mod:`pint` package and the :mod:`iam-units` registry
+        (`read the docs <https://github.com/IAMconsortium/units>`_) to handle units.
+        :mod:`pyam` will keep notation consistent with the input format (if possible)
+        and otherwise uses abbreviated units :code:`'{:~}'.format(u)` (see
+        `here <https://pint.readthedocs.io/en/stable/tutorial.html#string-formatting>`_
+        for more information).
+
+        As a result, the notation of returned units may differ from the input format.
+        For example, the unit :code:`EJ/yr` may be reformatted to :code:`EJ / a`.
         """
-        _value = _op_data(self, name, "subtract", axis=axis, fillna=fillna, a=a, b=b)
+        kwds = dict(axis=axis, fillna=fillna, ignore_units=ignore_units)
+        _value = _op_data(self, name, "subtract", **kwds, a=a, b=b)
         if append:
             self.append(_value, inplace=True)
         else:
             return IamDataFrame(_value, meta=self.meta)
 
-    def multiply(self, a, b, name, axis="variable", fillna=None, append=False):
-        """Compute the product of timeseries data between `a` and `b` along an `axis`
+    def multiply(
+        self, a, b, name, axis="variable", fillna=None, ignore_units=False, append=False
+    ):
+        """Multiply timeseries data items `a` and `b` along an `axis`
 
         This function computes `a * b`. If `a` or `b` are lists, the method applies
         :meth:`pandas.groupby().sum() <pandas.core.groupby.GroupBy.sum>` on each group.
@@ -1863,6 +1928,9 @@ class IamDataFrame(object):
         fillna : dict or scalar, optional
             Value to fill holes when rows are not defined for either `a` or `b`.
             Can be a scalar or a dictionary of the form :code:`{arg: default}`.
+        ignore_units : bool or str, optional
+            Perform operation on values without considering units. Set units of returned
+            data to `unknown` (if True) or the value of `ignore_units` (if str).
         append : bool, optional
             Whether to append aggregated timeseries data to this instance.
 
@@ -1870,15 +1938,35 @@ class IamDataFrame(object):
         -------
         :class:`IamDataFrame` or **None**
             Computed timeseries data or None if `append=True`.
+
+        See Also
+        --------
+        add, subtract, divide
+        apply : Apply a custom function on the timeseries data along any axis.
+
+        Notes
+        -----
+        This function uses the :mod:`pint` package and the :mod:`iam-units` registry
+        (`read the docs <https://github.com/IAMconsortium/units>`_) to handle units.
+        :mod:`pyam` will keep notation consistent with the input format (if possible)
+        and otherwise uses abbreviated units :code:`'{:~}'.format(u)` (see
+        `here <https://pint.readthedocs.io/en/stable/tutorial.html#string-formatting>`_
+        for more information).
+
+        As a result, the notation of returned units may differ from the input format.
+        For example, the unit :code:`EJ/yr` may be reformatted to :code:`EJ / a`.
         """
-        _value = _op_data(self, name, "multiply", axis=axis, fillna=fillna, a=a, b=b)
+        kwds = dict(axis=axis, fillna=fillna, ignore_units=ignore_units)
+        _value = _op_data(self, name, "multiply", **kwds, a=a, b=b)
         if append:
             self.append(_value, inplace=True)
         else:
             return IamDataFrame(_value, meta=self.meta)
 
-    def divide(self, a, b, name, axis="variable", fillna=None, append=False):
-        """Compute the division of timeseries data between `a` and `b` along an `axis`
+    def divide(
+        self, a, b, name, axis="variable", fillna=None, ignore_units=False, append=False
+    ):
+        """Divide the timeseries data items `a` and `b` along an `axis`
 
         This function computes `a / b`. If `a` or `b` are lists, the method applies
         :meth:`pandas.groupby().sum() <pandas.core.groupby.GroupBy.sum>` on each group.
@@ -1896,6 +1984,9 @@ class IamDataFrame(object):
         fillna : dict or scalar, optional
             Value to fill holes when rows are not defined for either `a` or `b`.
             Can be a scalar or a dictionary of the form :code:`{arg: default}`.
+        ignore_units : bool or str, optional
+            Perform operation on values without considering units. Set units of returned
+            data to `unknown` (if True) or the value of `ignore_units` (if str).
         append : bool, optional
             Whether to append aggregated timeseries data to this instance.
 
@@ -1903,8 +1994,26 @@ class IamDataFrame(object):
         -------
         :class:`IamDataFrame` or **None**
             Computed timeseries data or None if `append=True`.
+
+        See Also
+        --------
+        add, subtract, multiply
+        apply : Apply a custom function on the timeseries data along any axis.
+
+        Notes
+        -----
+        This function uses the :mod:`pint` package and the :mod:`iam-units` registry
+        (`read the docs <https://github.com/IAMconsortium/units>`_) to handle units.
+        :mod:`pyam` will keep notation consistent with the input format (if possible)
+        and otherwise uses abbreviated units :code:`'{:~}'.format(u)` (see
+        `here <https://pint.readthedocs.io/en/stable/tutorial.html#string-formatting>`_
+        for more information).
+
+        As a result, the notation of returned units may differ from the input format.
+        For example, the unit :code:`EJ/yr` may be reformatted to :code:`EJ / a`.
         """
-        _value = _op_data(self, name, "divide", axis=axis, fillna=fillna, a=a, b=b)
+        kwds = dict(axis=axis, fillna=fillna, ignore_units=ignore_units)
+        _value = _op_data(self, name, "divide", **kwds, a=a, b=b)
         if append:
             self.append(_value, inplace=True)
         else:
@@ -1944,6 +2053,17 @@ class IamDataFrame(object):
         -------
         :class:`IamDataFrame` or **None**
             Computed timeseries data or None if `append=True`.
+
+        Notes
+        -----
+        This function uses the :mod:`pint` package and the :mod:`iam-units` registry
+        (`read the docs <https://github.com/IAMconsortium/units>`_) to handle units.
+        :mod:`pyam` uses abbreviated units :code:`'{:~}'.format(u)` (see
+        `here <https://pint.readthedocs.io/en/stable/tutorial.html#string-formatting>`_
+        for more information).
+
+        As a result, the notation of returned units may differ from the input format.
+        For example, the unit :code:`EJ/yr` may be reformatted to :code:`EJ / a`.
         """
         _value = _op_data(self, name, func, axis=axis, fillna=fillna, args=args, **kwds)
         if append:
diff --git a/pyam/units.py b/pyam/units.py
index f6d5b15..1dbc072 100644
--- a/pyam/units.py
+++ b/pyam/units.py
@@ -36,7 +36,7 @@ def convert_unit(
     # Make versions without -equiv
     _current, _to = [i.replace("-equiv", "") for i in [current, to]]
     # Pair of (magnitude, unit)
-    qty = [ret.data.loc[where, "value"].values, _current]
+    qty = [ret._data.loc[where].values, _current]
 
     try:
         # Create a vector pint.Quantity and convert it ordinarily


Test patch:
diff --git a/tests/test_ops.py b/tests/test_ops.py
index 10dd4a9..1c0be34 100644
--- a/tests/test_ops.py
+++ b/tests/test_ops.py
@@ -1,6 +1,8 @@
 import pandas as pd
 import pytest
 import operator
+import pint
+from iam_units import registry
 
 from pyam import IamDataFrame
 from pyam.testing import assert_iamframe_equal
@@ -44,27 +46,35 @@ def test_add_raises(test_df_year):
 
 
 @pytest.mark.parametrize(
-    "arg, df_func, fillna",
+    "arg, df_func, fillna, ignore_units",
     (
-        ("Primary Energy|Coal", df_ops_variable, None),
-        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}),
-        ("Primary Energy|Coal", df_ops_variable_default, 5),
-        (2, df_ops_variable_number, None),
+        ("Primary Energy|Coal", df_ops_variable, None, False),
+        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}, "foo"),
+        ("Primary Energy|Coal", df_ops_variable_default, 5, "foo"),
+        (registry.Quantity(2, "EJ/yr"), df_ops_variable_number, None, False),
+        (2, df_ops_variable_number, None, "foo"),
     ),
 )
 @pytest.mark.parametrize("append", (False, True))
-def test_add_variable(test_df_year, arg, df_func, fillna, append):
+def test_add_variable(test_df_year, arg, df_func, fillna, ignore_units, append):
     """Verify that in-dataframe addition works on the default `variable` axis"""
 
-    exp = df_func(operator.add, "Sum", unit="EJ/yr", meta=test_df_year.meta)
+    unit = "EJ/yr" if ignore_units is False else ignore_units
+    exp = df_func(operator.add, "Sum", unit=unit, meta=test_df_year.meta)
 
     args = ("Primary Energy", arg, "Sum")
+    kwds = dict(ignore_units=ignore_units, fillna=fillna)
     if append:
         obs = test_df_year.copy()
-        obs.add(*args, fillna=fillna, append=True)
+        obs.add(*args, **kwds, append=True)
         assert_iamframe_equal(test_df_year.append(exp), obs)
     else:
-        assert_iamframe_equal(exp, test_df_year.add(*args, fillna=fillna))
+        # check that incompatible units raise the expected error
+        if ignore_units:
+            with pytest.raises(pint.DimensionalityError):
+                test_df_year.add(*args, fillna=fillna, ignore_units=False)
+
+        assert_iamframe_equal(exp, test_df_year.add(*args, **kwds))
 
 
 @pytest.mark.parametrize("append", (False, True))
@@ -91,27 +101,35 @@ def test_add_scenario(test_df_year, append):
 
 
 @pytest.mark.parametrize(
-    "arg, df_func, fillna",
+    "arg, df_func, fillna, ignore_units",
     (
-        ("Primary Energy|Coal", df_ops_variable, None),
-        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}),
-        ("Primary Energy|Coal", df_ops_variable_default, 5),
-        (2, df_ops_variable_number, None),
+        ("Primary Energy|Coal", df_ops_variable, None, False),
+        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}, "foo"),
+        ("Primary Energy|Coal", df_ops_variable_default, 5, "foo"),
+        (registry.Quantity(2, "EJ/yr"), df_ops_variable_number, None, False),
+        (2, df_ops_variable_number, None, "foo"),
     ),
 )
 @pytest.mark.parametrize("append", (False, True))
-def test_subtract_variable(test_df_year, arg, df_func, fillna, append):
+def test_subtract_variable(test_df_year, arg, df_func, fillna, append, ignore_units):
     """Verify that in-dataframe subtraction works on the default `variable` axis"""
 
-    exp = df_func(operator.sub, "Diff", unit="EJ/yr", meta=test_df_year.meta)
+    unit = "EJ/yr" if ignore_units is False else ignore_units
+    exp = df_func(operator.sub, "Diff", unit=unit, meta=test_df_year.meta)
 
+    args = ("Primary Energy", arg, "Diff")
+    kwds = dict(ignore_units=ignore_units, fillna=fillna)
     if append:
         obs = test_df_year.copy()
-        obs.subtract("Primary Energy", arg, "Diff", fillna=fillna, append=True)
+        obs.subtract(*args, **kwds, append=True)
         assert_iamframe_equal(test_df_year.append(exp), obs)
     else:
-        obs = test_df_year.subtract("Primary Energy", arg, "Diff", fillna=fillna)
-        assert_iamframe_equal(exp, obs)
+        # check that incompatible units raise the expected error
+        if ignore_units:
+            with pytest.raises(pint.DimensionalityError):
+                test_df_year.add(*args, fillna=fillna, ignore_units=False)
+
+        assert_iamframe_equal(exp, test_df_year.subtract(*args, **kwds))
 
 
 @pytest.mark.parametrize("append", (False, True))
@@ -138,27 +156,38 @@ def test_subtract_scenario(test_df_year, append):
 
 
 @pytest.mark.parametrize(
-    "arg, df_func, fillna",
+    "arg, df_func, fillna, ignore_units",
     (
-        ("Primary Energy|Coal", df_ops_variable, None),
-        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}),
-        ("Primary Energy|Coal", df_ops_variable_default, 5),
-        (2, df_ops_variable_number, None),
+        ("Primary Energy|Coal", df_ops_variable, None, False),
+        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}, "foo"),
+        ("Primary Energy|Coal", df_ops_variable_default, 5, "foo"),
+        # note that multiplying with pint reformats the unit
+        (2, df_ops_variable_number, None, False),
     ),
 )
 @pytest.mark.parametrize("append", (False, True))
-def test_multiply_variable(test_df_year, arg, df_func, fillna, append):
+def test_multiply_variable(test_df_year, arg, df_func, fillna, ignore_units, append):
     """Verify that in-dataframe addition works on the default `variable` axis"""
 
-    exp = df_func(operator.mul, "Prod", unit="EJ/yr", meta=test_df_year.meta)
+    if ignore_units:
+        unit = ignore_units
+    else:
+        unit = "EJ / a" if isinstance(arg, int) else "EJ ** 2 / a ** 2"
+    exp = df_func(operator.mul, "Prod", unit=unit, meta=test_df_year.meta)
 
     args = ("Primary Energy", arg, "Prod")
+    kwds = dict(ignore_units=ignore_units, fillna=fillna)
     if append:
         obs = test_df_year.copy()
-        obs.multiply(*args, fillna=fillna, append=True)
+        obs.multiply(*args, **kwds, append=True)
         assert_iamframe_equal(test_df_year.append(exp), obs)
     else:
-        assert_iamframe_equal(exp, test_df_year.multiply(*args, fillna=fillna))
+        # check that incompatible units raise the expected error
+        if ignore_units:
+            with pytest.raises(pint.DimensionalityError):
+                test_df_year.add(*args, fillna=fillna, ignore_units=False)
+
+        assert_iamframe_equal(exp, test_df_year.multiply(*args, **kwds))
 
 
 @pytest.mark.parametrize("append", (False, True))
@@ -172,7 +201,7 @@ def test_multiply_scenario(test_df_year, append):
         scenario=v[2],
         region="World",
         variable="Primary Energy",
-        unit="EJ/yr",
+        unit="EJ ** 2 / a ** 2",
     )
 
     if append:
@@ -185,27 +214,39 @@ def test_multiply_scenario(test_df_year, append):
 
 
 @pytest.mark.parametrize(
-    "arg, df_func, fillna",
+    "arg, df_func, fillna, ignore_units",
     (
-        ("Primary Energy|Coal", df_ops_variable, None),
-        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}),
-        ("Primary Energy|Coal", df_ops_variable_default, 5),
-        (2, df_ops_variable_number, None),
+        ("Primary Energy|Coal", df_ops_variable, None, False),
+        ("Primary Energy|Coal", df_ops_variable_default, {"c": 7, "b": 5}, "foo"),
+        ("Primary Energy|Coal", df_ops_variable_default, 5, "foo"),
+        (registry.Quantity(2, "EJ/yr"), df_ops_variable_number, None, False),
+        (2, df_ops_variable_number, None, False),
     ),
 )
 @pytest.mark.parametrize("append", (False, True))
-def test_divide_variable(test_df_year, arg, df_func, fillna, append):
+def test_divide_variable(test_df_year, arg, df_func, fillna, append, ignore_units):
     """Verify that in-dataframe addition works on the default `variable` axis"""
 
-    exp = df_func(operator.truediv, "Ratio", unit="EJ/yr", meta=test_df_year.meta)
+    # note that dividing with pint reformats the unit
+    if ignore_units:
+        unit = ignore_units
+    else:
+        unit = "EJ / a" if isinstance(arg, int) else ""
+    exp = df_func(operator.truediv, "Ratio", unit=unit, meta=test_df_year.meta)
 
     args = ("Primary Energy", arg, "Ratio")
+    kwds = dict(ignore_units=ignore_units, fillna=fillna)
     if append:
         obs = test_df_year.copy()
-        obs.divide(*args, fillna=fillna, append=True)
+        obs.divide(*args, **kwds, append=True)
         assert_iamframe_equal(test_df_year.append(exp), obs)
     else:
-        assert_iamframe_equal(exp, test_df_year.divide(*args, fillna=fillna))
+        # check that incompatible units raise the expected error
+        if ignore_units:
+            with pytest.raises(pint.DimensionalityError):
+                test_df_year.add(*args, fillna=fillna, ignore_units=False)
+
+        assert_iamframe_equal(exp, test_df_year.divide(*args, **kwds))
 
 
 @pytest.mark.parametrize("append", (False, True))
@@ -219,7 +260,7 @@ def test_divide_scenario(test_df_year, append):
         scenario=v[2],
         region="World",
         variable="Primary Energy",
-        unit="EJ/yr",
+        unit="",
     )
 
     if append:
@@ -232,39 +273,34 @@ def test_divide_scenario(test_df_year, append):
 
 
 @pytest.mark.parametrize("append", (False, True))
-def test_apply_variable(plot_stackplot_df, append):
+def test_apply_variable(test_df_year, append):
     """Verify that in-dataframe apply works on the default `variable` axis"""
 
-    def custom_func(a, b, c, d, e):
-        return a / c + b / d + e
+    def custom_func(a, b, c, d):
+        return a * b + c * d
+
+    v = "new variable"
 
-    args = ["Emissions|CO2|Tar", "Emissions|CO2|Cars", "Emissions|CO2|LUC"]
-    kwds = {"d": "Emissions|CO2|Agg", "e": 5}
     exp = IamDataFrame(
         pd.DataFrame(
-            [
-                0.3 / (-0.3) + 1.6 / 0.5 + 5,
-                0.35 / (-0.6) + 3.8 / (-0.1) + 5,
-                0.35 / (-1.2) + 3.0 / (-0.5) + 5,
-                0.33 / (-1.0) + 2.5 / (-0.7) + 5,
-            ],
-            index=[2005, 2010, 2015, 2020],
+            [custom_func(1, 2, 0.5, 3), custom_func(6, 2, 3, 3)], index=[2005, 2010]
         ).T,
-        model="IMG",
-        scenario="a_scen",
-        region="World",
-        variable="new variable",
-        unit="Mt CO2/yr",
+        **DF_ARGS,
+        scenario="scen_a",
+        variable=v,
+        unit="EJ / a",  # applying operations with pint reformats the unit
+        meta=test_df_year.meta,
     )
 
+    args = ["Primary Energy", 2]
+    kwds = dict(d=3, c="Primary Energy|Coal")
+
     if append:
-        obs = plot_stackplot_df.copy()
+        obs = test_df_year.copy()
         obs.apply(custom_func, name="new variable", append=True, args=args, **kwds)
-        assert_iamframe_equal(plot_stackplot_df.append(exp), obs)
+        assert_iamframe_equal(test_df_year.append(exp), obs)
     else:
-        obs = plot_stackplot_df.apply(
-            custom_func, name="new variable", args=args, **kwds
-        )
+        obs = test_df_year.apply(custom_func, name=v, args=args, **kwds)
         assert_iamframe_equal(exp, obs)