diff --git a/doc/sphinx/src/databox.rst b/doc/sphinx/src/databox.rst
index 624e187cf..fd4114161 100644
--- a/doc/sphinx/src/databox.rst
+++ b/doc/sphinx/src/databox.rst
@@ -14,7 +14,13 @@ To use databox, simply include the relevant header:
 
   #include <databox.hpp>
 
-``DatBox`` is templated on underyling data type, which defaults to the
+``DataBox`` Templates
+^^^^^^^^^^^^^^^^^^^^^
+
+Underlying Data Type
+--------------------
+
+``DataBox`` is templated on the underyling data type, which defaults to the
 ``Real`` type provided by ``ports-of-call``. (This is usually a
 ``double``.)
 
@@ -30,6 +36,9 @@ single type, you may wish to declare a type alias such as:
 
    using DataBox = Spiner::DataBox<double>
 
+Interpolation Gridding
+----------------------
+
 Spiner is also templated on how the interpolation gridding works. This
 template parameter is called ``Grid_t``. The available options at this time are:
 
@@ -47,6 +56,46 @@ as, for example:
 More detail on the interpolation gridding is available below and in
 the interpolation section.
 
+Transformations
+---------------
+
+Spiner performs linear interpolation of the data, but can apply transformations
+to the data in order to enable interpolation methods like log-log, log-lin,
+etc.  Spiner will perform linear interpolation of the transformed variables.
+
+When constructing the gridding, ``Spiner::RegularGrid1D<T>`` and
+``Spiner::PiecewiseGrid1D<T>`` are templated on the independent variable
+transformation.  Note that all independent variables will use the same
+transformation.  ``DataBox`` is templated on the dependent variable
+transformation.  In all cases, the default is a linear "transformation" (no
+transformation to the data), so if you do not specify the transformation then
+you get linear interpolation.  The available transformations are in
+``spiner/transformations.hpp`` and users can write custom transformations
+following the examples there.
+
+For example, if the user wants to transform the independent variables with a
+logarithm and the dependent variable with a custom arctangent transformation,
+the declaration of the ``DataBox`` might look like
+
+.. code-block:: cpp
+
+   using DataBox = Spiner::DataBox<double, Spiner::RegularGrid1D<double,
+   Spiner::TransformLogarithmic>, CustomArctanTransform>;
+
+This would result in interpolation according to the equation
+
+.. math::
+
+   v &= \arctan(y)
+
+   u_i &= \log(x_i)
+
+   v &= b + \sum m_i u_i
+
+where :math:`x_i` are the independent variables.  By convention, the Spiner
+documentation uses :math:`u_i` as the transformed independent variables, and
+:math:`v` as the transformed dependent variable.
+
 .. note::
    In C++17 and later, you can also get the default type specialization
    by simply omitting the template arguments.
diff --git a/spiner/databox.hpp b/spiner/databox.hpp
index 793664364..caa752342 100644
--- a/spiner/databox.hpp
+++ b/spiner/databox.hpp
@@ -51,6 +51,7 @@ enum class DataStatus {
 enum class AllocationTarget { Host, Device };
 
 template <typename T = Real, typename Grid_t = RegularGrid1D<T>,
+          typename Transform = TransformLinear,
           typename Concept =
               typename std::enable_if<std::is_arithmetic<T>::value, bool>::type>
 class DataBox {
@@ -101,7 +102,7 @@ class DataBox {
     setAllIndexed_();
   }
   PORTABLE_INLINE_FUNCTION
-  DataBox(const DataBox<T, Grid_t, Concept> &src) noexcept
+  DataBox(const DataBox<T, Grid_t, Transform, Concept> &src) noexcept
       : rank_(src.rank_), status_(src.status_), data_(src.data_) {
     setAllIndexed_();
     dataView_.InitWithShallowSlice(src.dataView_, 6, 0, src.dim(6));
@@ -113,8 +114,8 @@ class DataBox {
 
   // Slice constructor
   PORTABLE_INLINE_FUNCTION
-  DataBox(const DataBox<T, Grid_t, Concept> &b, const int dim, const int indx,
-          const int nvar) noexcept
+  DataBox(const DataBox<T, Grid_t, Transform, Concept> &b, const int dim,
+          const int indx, const int nvar) noexcept
       : status_(DataStatus::Unmanaged), data_(b.data_) {
     dataView_.InitWithShallowSlice(b.dataView_, dim, indx, nvar);
     rank_ = dataView_.GetRank();
@@ -150,29 +151,44 @@ class DataBox {
     resize(AllocationTarget::Host, std::forward<Args>(args)...);
   }
 
+  // Data array accessors
+  template <typename... Args>
+  PORTABLE_INLINE_FUNCTION T get_data_value(Args... args) const {
+    return Transform::reverse(dataView_(std::forward<Args>(args)...));
+  }
+  template <typename... Args>
+  PORTABLE_INLINE_FUNCTION void set_data_value(const T new_value,
+                                               Args... args) const {
+    dataView_(std::forward<Args>(args)...) = Transform::forward(new_value);
+  }
+
   // Index operators
   // examle calls:
   // T x = db(n4, n3, n2, n1);
   template <typename... Args>
-  PORTABLE_INLINE_FUNCTION T &operator()(Args... args) {
-    return dataView_(std::forward<Args>(args)...);
-  }
-  template <typename... Args>
   PORTABLE_INLINE_FUNCTION T &operator()(Args... args) const {
+    static_assert(
+        std::is_same<Transform, TransformLinear>::value,
+        "DataBox::operator() is only available if the independent variable "
+        "data transform is Spiner::TransformLinear, because otherwise "
+        "accessing the underlying data is not doing what you expect.  Use "
+        "DataBox::get_data_value and DataBox::set_data_value instead, because "
+        "those correctly account for the independent value data "
+        "transformation.");
     return dataView_(std::forward<Args>(args)...);
   }
 
   // Slice operation
   PORTABLE_INLINE_FUNCTION
-  DataBox<T, Grid_t, Concept> slice(const int dim, const int indx,
-                                    const int nvar) const {
+  DataBox<T, Grid_t, Transform, Concept> slice(const int dim, const int indx,
+                                               const int nvar) const {
     return DataBox(*this, dim, indx, nvar);
   }
-  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Concept>
+  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Transform, Concept>
   slice(const int indx) const {
     return slice(rank_, indx, 1);
   }
-  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Concept>
+  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Transform, Concept>
   slice(const int ix2, const int ix1) const {
     // DataBox a(*this, rank_, ix2, 1);
     // return DataBox(a, a.rank_, ix1, 1);
@@ -213,13 +229,14 @@ class DataBox {
   // WARNING: requires memory to be pre-allocated.
   // TODO: add 3d and higher interpFromDB if necessary
   PORTABLE_INLINE_FUNCTION void
-  interpFromDB(const DataBox<T, Grid_t, Concept> &db, const T x);
+  interpFromDB(const DataBox<T, Grid_t, Transform, Concept> &db, const T x);
   PORTABLE_INLINE_FUNCTION void
-  interpFromDB(const DataBox<T, Grid_t, Concept> &db, const T x2, const T x1);
+  interpFromDB(const DataBox<T, Grid_t, Transform, Concept> &db, const T x2,
+               const T x1);
   template <typename... Args>
-  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Concept>
+  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Transform, Concept>
   interpToDB(Args... args) {
-    DataBox<T, Grid_t, Concept> db;
+    DataBox<T, Grid_t, Transform, Concept> db;
     db.interpFromDB(*this, std::forward<Args>(args)...);
     return db;
   }
@@ -245,10 +262,11 @@ class DataBox {
   // Does no checks that memory is available.
   // Optionally copies shape of source with ndims fewer slowest-moving
   // dimensions
-  PORTABLE_INLINE_FUNCTION void copyShape(const DataBox<T, Grid_t, Concept> &db,
-                                          const int ndims = 0);
+  PORTABLE_INLINE_FUNCTION void
+  copyShape(const DataBox<T, Grid_t, Transform, Concept> &db,
+            const int ndims = 0);
   // Returns new databox with same memory and metadata
-  inline void copyMetadata(const DataBox<T, Grid_t, Concept> &src);
+  inline void copyMetadata(const DataBox<T, Grid_t, Transform, Concept> &src);
 
 #ifdef SPINER_USE_HDF
   inline herr_t saveHDF() const { return saveHDF(SP5::DB::FILENAME); }
@@ -264,9 +282,9 @@ class DataBox {
   inline Grid_t &range(const int i) { return grids_[i]; }
 
   // Assignment and move, both perform shallow copy
-  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Concept> &
-  operator=(const DataBox<T, Grid_t, Concept> &other);
-  inline void copy(const DataBox<T, Grid_t, Concept> &src);
+  PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Transform, Concept> &
+  operator=(const DataBox<T, Grid_t, Transform, Concept> &other);
+  inline void copy(const DataBox<T, Grid_t, Transform, Concept> &src);
 
   // utility info
   inline DataStatus dataStatus() const { return status_; }
@@ -274,8 +292,10 @@ class DataBox {
   inline bool ownsAllocatedMemory() {
     return (status_ != DataStatus::Unmanaged);
   }
-  inline bool operator==(const DataBox<T, Grid_t, Concept> &other) const;
-  inline bool operator!=(const DataBox<T, Grid_t, Concept> &other) const {
+  inline bool
+  operator==(const DataBox<T, Grid_t, Transform, Concept> &other) const;
+  inline bool
+  operator!=(const DataBox<T, Grid_t, Transform, Concept> &other) const {
     return !(*this == other);
   }
 
@@ -370,11 +390,11 @@ class DataBox {
   }
   // ------------------------------------
 
-  DataBox<T, Grid_t, Concept>
+  DataBox<T, Grid_t, Transform, Concept>
   getOnDevice() const { // getOnDevice is always a deep copy
     if (size() == 0 ||
         status_ == DataStatus::Empty) { // edge case for unallocated
-      DataBox<T, Grid_t, Concept> a;
+      DataBox<T, Grid_t, Transform, Concept> a;
       return a;
     }
     // create device memory (host memory if no device)
@@ -382,8 +402,8 @@ class DataBox {
     // copy to device
     portableCopyToDevice(device_data, data_, sizeBytes());
     // create new databox of size size
-    DataBox<T, Grid_t, Concept> a{device_data, dim(6), dim(5), dim(4),
-                                  dim(3),      dim(2), dim(1)};
+    DataBox<T, Grid_t, Transform, Concept> a{
+        device_data, dim(6), dim(5), dim(4), dim(3), dim(2), dim(1)};
     a.copyShape(*this);
     // set correct allocation status of the new databox
     // note this is ALWAYS device, even if host==device.
@@ -439,25 +459,28 @@ class DataBox {
 };
 
 // Read an array, shallow
-template <typename T, typename Grid_t, typename Concept>
-inline void DataBox<T, Grid_t, Concept>::setArray(PortableMDArray<T> &A) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline void
+DataBox<T, Grid_t, Transform, Concept>::setArray(PortableMDArray<T> &A) {
   dataView_ = A;
   rank_ = A.GetRank();
   setAllIndexed_();
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 PORTABLE_INLINE_FUNCTION T
-DataBox<T, Grid_t, Concept>::interpToReal(const T x) const noexcept {
+DataBox<T, Grid_t, Transform, Concept>::interpToReal(const T x) const noexcept {
   assert(canInterpToReal_(1));
   int ix;
   weights_t<T> w;
   grids_[0].weights(x, ix, w);
-  return w[0] * dataView_(ix) + w[1] * dataView_(ix + 1);
+  const auto v = w[0] * dataView_(ix) + w[1] * dataView_(ix + 1);
+  return Transform::reverse(v);
 }
 
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_FORCEINLINE_FUNCTION T
+DataBox<T, Grid_t, Transform, Concept>::interpToReal(
     const T x2, const T x1) const noexcept {
   assert(canInterpToReal_(2));
   int ix1, ix2;
@@ -466,14 +489,16 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
   grids_[1].weights(x2, ix2, w2);
   // TODO: prefectch corners for speed?
   // TODO: re-order access pattern?
-  return (w2[0] *
-              (w1[0] * dataView_(ix2, ix1) + w1[1] * dataView_(ix2, ix1 + 1)) +
-          w2[1] * (w1[0] * dataView_(ix2 + 1, ix1) +
-                   w1[1] * dataView_(ix2 + 1, ix1 + 1)));
+  const auto v =
+      (w2[0] * (w1[0] * dataView_(ix2, ix1) + w1[1] * dataView_(ix2, ix1 + 1)) +
+       w2[1] * (w1[0] * dataView_(ix2 + 1, ix1) +
+                w1[1] * dataView_(ix2 + 1, ix1 + 1)));
+  return Transform::reverse(v);
 }
 
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_FORCEINLINE_FUNCTION T
+DataBox<T, Grid_t, Transform, Concept>::interpToReal(
     const T x3, const T x2, const T x1) const noexcept {
   assert(canInterpToReal_(3));
   int ix[3];
@@ -483,7 +508,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
   grids_[2].weights(x3, ix[2], w[2]);
   // TODO: prefect corners for speed?
   // TODO: re-order access pattern?
-  return (
+  const auto v = (
       w[2][0] * (w[1][0] * (w[0][0] * dataView_(ix[2], ix[1], ix[0]) +
                             w[0][1] * dataView_(ix[2], ix[1], ix[0] + 1)) +
                  w[1][1] * (w[0][0] * dataView_(ix[2], ix[1] + 1, ix[0]) +
@@ -493,10 +518,12 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
                       w[0][1] * dataView_(ix[2] + 1, ix[1], ix[0] + 1)) +
            w[1][1] * (w[0][0] * dataView_(ix[2] + 1, ix[1] + 1, ix[0]) +
                       w[0][1] * dataView_(ix[2] + 1, ix[1] + 1, ix[0] + 1))));
+  return Transform::reverse(v);
 }
 
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_FORCEINLINE_FUNCTION T
+DataBox<T, Grid_t, Transform, Concept>::interpToReal(
     const T x3, const T x2, const T x1, const int idx) const noexcept {
   assert(rank_ == 4);
   for (int r = 1; r < rank_; ++r) {
@@ -510,7 +537,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
   grids_[3].weights(x3, ix[2], w[2]);
   // TODO: prefect corners for speed?
   // TODO: re-order access pattern?
-  return (
+  const auto v = (
       w[2][0] *
           (w[1][0] * (w[0][0] * dataView_(ix[2], ix[1], ix[0], idx) +
                       w[0][1] * dataView_(ix[2], ix[1], ix[0] + 1, idx)) +
@@ -522,12 +549,14 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
            w[1][1] *
                (w[0][0] * dataView_(ix[2] + 1, ix[1] + 1, ix[0], idx) +
                 w[0][1] * dataView_(ix[2] + 1, ix[1] + 1, ix[0] + 1, idx))));
+  return Transform::reverse(v);
 }
 
 // DH: this is a large function to force an inline, perhaps just make it a
 // suggestion to the compiler?
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_FORCEINLINE_FUNCTION T
+DataBox<T, Grid_t, Transform, Concept>::interpToReal(
     const T x4, const T x3, const T x2, const T x1) const noexcept {
   assert(canInterpToReal_(4));
   T x[] = {x1, x2, x3, x4};
@@ -539,7 +568,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
   // TODO(JMM): This is getty pretty gross. Should we automate?
   // Hand-written is probably faster, though.
   // Breaking line-limit to make this easier to read
-  return (
+  const auto v = (
       w[3][0] *
           (w[2][0] *
                (w[1][0] *
@@ -575,10 +604,12 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
                                                      ix[1] + 1, ix[0] + 1))))
 
   );
+  return Transform::reverse(v);
 }
 
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_FORCEINLINE_FUNCTION T
+DataBox<T, Grid_t, Transform, Concept>::interpToReal(
     const T x4, const T x3, const T x2, const int idx,
     const T x1) const noexcept {
   assert(rank_ == 5);
@@ -598,7 +629,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
   // TODO(JMM): This is getty pretty gross. Should we automate?
   // Hand-written is probably faster, though.
   // Breaking line-limit to make this easier to read
-  return (
+  const auto v = (
       w[3][0] *
           (w[2][0] *
                (w[1][0] *
@@ -638,12 +669,13 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox<T, Grid_t, Concept>::interpToReal(
                                                idx, ix[0] + 1))))
 
   );
+  return Transform::reverse(v);
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 PORTABLE_INLINE_FUNCTION void
-DataBox<T, Grid_t, Concept>::interpFromDB(const DataBox<T, Grid_t, Concept> &db,
-                                          const T x) {
+DataBox<T, Grid_t, Transform, Concept>::interpFromDB(
+    const DataBox<T, Grid_t, Transform, Concept> &db, const T x) {
   assert(db.indices_[db.rank_ - 1] == IndexType::Interpolated);
   assert(db.grids_[db.rank_ - 1].isWellFormed());
   assert(size() == (db.size() / db.dim(db.rank_)));
@@ -653,7 +685,8 @@ DataBox<T, Grid_t, Concept>::interpFromDB(const DataBox<T, Grid_t, Concept> &db,
   copyShape(db, 1);
 
   db.grids_[db.rank_ - 1].weights(x, ix, w);
-  DataBox<T, Grid_t, Concept> lower(db.slice(ix)), upper(db.slice(ix + 1));
+  DataBox<T, Grid_t, Transform, Concept> lower(db.slice(ix)),
+      upper(db.slice(ix + 1));
   // lower = db.slice(ix);
   // upper = db.slice(ix+1);
   for (int i = 0; i < size(); i++) {
@@ -661,10 +694,10 @@ DataBox<T, Grid_t, Concept>::interpFromDB(const DataBox<T, Grid_t, Concept> &db,
   }
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 PORTABLE_INLINE_FUNCTION void
-DataBox<T, Grid_t, Concept>::interpFromDB(const DataBox<T, Grid_t, Concept> &db,
-                                          const T x2, const T x1) {
+DataBox<T, Grid_t, Transform, Concept>::interpFromDB(
+    const DataBox<T, Grid_t, Transform, Concept> &db, const T x2, const T x1) {
   assert(db.rank_ >= 2);
   assert(db.indices_[db.rank_ - 1] == IndexType::Interpolated);
   assert(db.grids_[db.rank_ - 1].isWellFormed());
@@ -678,7 +711,7 @@ DataBox<T, Grid_t, Concept>::interpFromDB(const DataBox<T, Grid_t, Concept> &db,
 
   db.grids_[db.rank_ - 2].weights(x1, ix1, w1);
   db.grids_[db.rank_ - 1].weights(x2, ix2, w2);
-  DataBox<T, Grid_t, Concept> corners[2][2]{
+  DataBox<T, Grid_t, Transform, Concept> corners[2][2]{
       {db.slice(ix2, ix1), db.slice(ix2 + 1, ix1)},
       {db.slice(ix2, ix1 + 1), db.slice(ix2 + 1, ix1 + 1)}};
   //    copyShape(db,2);
@@ -707,10 +740,9 @@ DataBox<T, Grid_t, Concept>::interpFromDB(const DataBox<T, Grid_t, Concept> &db,
 // Reshapes from other databox, but does not allocate memory.
 // Does no checks that memory is available.
 // Optionally copies shape of source with ndims fewer slowest-moving dimensions
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_INLINE_FUNCTION void
-DataBox<T, Grid_t, Concept>::copyShape(const DataBox<T, Grid_t, Concept> &db,
-                                       const int ndims) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_INLINE_FUNCTION void DataBox<T, Grid_t, Transform, Concept>::copyShape(
+    const DataBox<T, Grid_t, Transform, Concept> &db, const int ndims) {
   rank_ = db.rank_ - ndims;
   int dims[MAXRANK];
   for (int i = 0; i < MAXRANK; i++)
@@ -727,9 +759,9 @@ DataBox<T, Grid_t, Concept>::copyShape(const DataBox<T, Grid_t, Concept> &db,
 }
 // reallocates and then copies shape from other databox
 // everything but the actual copy in a deep copy
-template <typename T, typename Grid_t, typename Concept>
-inline void DataBox<T, Grid_t, Concept>::copyMetadata(
-    const DataBox<T, Grid_t, Concept> &src) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline void DataBox<T, Grid_t, Transform, Concept>::copyMetadata(
+    const DataBox<T, Grid_t, Transform, Concept> &src) {
   AllocationTarget t =
       (src.status_ == DataStatus::AllocatedDevice ? AllocationTarget::Device
                                                   : AllocationTarget::Host);
@@ -743,9 +775,9 @@ inline void DataBox<T, Grid_t, Concept>::copyMetadata(
 }
 
 #ifdef SPINER_USE_HDF
-template <typename T, typename Grid_t, typename Concept>
-inline herr_t
-DataBox<T, Grid_t, Concept>::saveHDF(const std::string &filename) const {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline herr_t DataBox<T, Grid_t, Transform, Concept>::saveHDF(
+    const std::string &filename) const {
   herr_t status;
   hid_t file;
 
@@ -755,10 +787,9 @@ DataBox<T, Grid_t, Concept>::saveHDF(const std::string &filename) const {
   return status;
 }
 
-template <typename T, typename Grid_t, typename Concept>
-inline herr_t
-DataBox<T, Grid_t, Concept>::saveHDF(hid_t loc,
-                                     const std::string &groupname) const {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline herr_t DataBox<T, Grid_t, Transform, Concept>::saveHDF(
+    hid_t loc, const std::string &groupname) const {
   hid_t group, grids;
   herr_t status = 0;
   static_assert(std::is_same<T, double>::value || std::is_same<T, float>::value,
@@ -817,18 +848,19 @@ DataBox<T, Grid_t, Concept>::saveHDF(hid_t loc,
   return status;
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 inline herr_t
-DataBox<T, Grid_t, Concept>::loadHDF(const std::string &filename) {
+DataBox<T, Grid_t, Transform, Concept>::loadHDF(const std::string &filename) {
   herr_t status;
   hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
   status = loadHDF(file, SP5::DB::GRPNAME);
   return status;
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 inline herr_t
-DataBox<T, Grid_t, Concept>::loadHDF(hid_t loc, const std::string &groupname) {
+DataBox<T, Grid_t, Transform, Concept>::loadHDF(hid_t loc,
+                                                const std::string &groupname) {
   hid_t group, grids;
   herr_t status = 0;
   std::vector<int> index_types;
@@ -883,9 +915,10 @@ DataBox<T, Grid_t, Concept>::loadHDF(hid_t loc, const std::string &groupname) {
 #endif // SPINER_USE_HDF
 
 // Performs shallow copy by default
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Concept> &
-DataBox<T, Grid_t, Concept>::operator=(const DataBox<T, Grid_t, Concept> &src) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_INLINE_FUNCTION DataBox<T, Grid_t, Transform, Concept> &
+DataBox<T, Grid_t, Transform, Concept>::operator=(
+    const DataBox<T, Grid_t, Transform, Concept> &src) {
   if (this != &src) {
     rank_ = src.rank_;
     status_ = src.status_;
@@ -900,17 +933,17 @@ DataBox<T, Grid_t, Concept>::operator=(const DataBox<T, Grid_t, Concept> &src) {
 }
 
 // Performs a deep copy
-template <typename T, typename Grid_t, typename Concept>
-inline void
-DataBox<T, Grid_t, Concept>::copy(const DataBox<T, Grid_t, Concept> &src) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline void DataBox<T, Grid_t, Transform, Concept>::copy(
+    const DataBox<T, Grid_t, Transform, Concept> &src) {
   copyMetadata(src);
   for (int i = 0; i < src.size(); i++)
     dataView_(i) = src(i);
 }
 
-template <typename T, typename Grid_t, typename Concept>
-inline bool DataBox<T, Grid_t, Concept>::operator==(
-    const DataBox<T, Grid_t, Concept> &other) const {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline bool DataBox<T, Grid_t, Transform, Concept>::operator==(
+    const DataBox<T, Grid_t, Transform, Concept> &other) const {
   if (rank_ != other.rank_) return false;
   for (int i = 0; i < rank_; i++) {
     if (indices_[i] != other.indices_[i]) return false;
@@ -923,42 +956,44 @@ inline bool DataBox<T, Grid_t, Concept>::operator==(
 }
 
 // TODO: should this be std::reduce?
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_INLINE_FUNCTION T DataBox<T, Grid_t, Concept>::min() const {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_INLINE_FUNCTION T DataBox<T, Grid_t, Transform, Concept>::min() const {
   T min = std::numeric_limits<T>::infinity();
   for (int i = 0; i < size(); i++) {
-    min = std::min(min, dataView_(i));
+    min = std::min(min, Transform::reverse(dataView_(i)));
   }
   return min;
 }
 
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_INLINE_FUNCTION T DataBox<T, Grid_t, Concept>::max() const {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_INLINE_FUNCTION T DataBox<T, Grid_t, Transform, Concept>::max() const {
   T max = -std::numeric_limits<T>::infinity();
   for (int i = 0; i < size(); i++) {
-    max = std::max(max, dataView_(i));
+    max = std::max(max, Transform::reverse(dataView_(i)));
   }
   return max;
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 PORTABLE_INLINE_FUNCTION Grid_t
-DataBox<T, Grid_t, Concept>::range(int i) const {
+DataBox<T, Grid_t, Transform, Concept>::range(int i) const {
   assert(0 <= i && i < rank_);
   assert(indices_[i] == IndexType::Interpolated);
   return grids_[i];
 }
 
-template <typename T, typename Grid_t, typename Concept>
-PORTABLE_INLINE_FUNCTION void DataBox<T, Grid_t, Concept>::setAllIndexed_() {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+PORTABLE_INLINE_FUNCTION void
+DataBox<T, Grid_t, Transform, Concept>::setAllIndexed_() {
   for (int i = 0; i < rank_; i++) {
     indices_[i] = IndexType::Indexed;
   }
 }
 
-template <typename T, typename Grid_t, typename Concept>
+template <typename T, typename Grid_t, typename Transform, typename Concept>
 PORTABLE_INLINE_FUNCTION bool
-DataBox<T, Grid_t, Concept>::canInterpToReal_(const int interpOrder) const {
+DataBox<T, Grid_t, Transform, Concept>::canInterpToReal_(
+    const int interpOrder) const {
   if (rank_ != interpOrder) return false;
   for (int i = 0; i < rank_; i++) {
     if (indices_[i] != IndexType::Interpolated) return false;
@@ -967,13 +1002,13 @@ DataBox<T, Grid_t, Concept>::canInterpToReal_(const int interpOrder) const {
   return true;
 }
 
-template <typename T, typename Grid_t, typename Concept>
-inline DataBox<T, Grid_t, Concept>
-getOnDeviceDataBox(const DataBox<T, Grid_t, Concept> &a_host) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline DataBox<T, Grid_t, Transform, Concept>
+getOnDeviceDataBox(const DataBox<T, Grid_t, Transform, Concept> &a_host) {
   return a_host.getOnDevice();
 }
-template <typename T, typename Grid_t, typename Concept>
-inline void free(DataBox<T, Grid_t, Concept> &db) {
+template <typename T, typename Grid_t, typename Transform, typename Concept>
+inline void free(DataBox<T, Grid_t, Transform, Concept> &db) {
   db.finalize();
 }
 
diff --git a/spiner/piecewise_grid_1d.hpp b/spiner/piecewise_grid_1d.hpp
index 49501c964..206a87ab1 100644
--- a/spiner/piecewise_grid_1d.hpp
+++ b/spiner/piecewise_grid_1d.hpp
@@ -35,6 +35,7 @@
 namespace Spiner {
 
 template <typename T = Real, int NGRIDSMAX = 5,
+          typename Transform = TransformLinear,
           typename =
               typename std::enable_if<std::is_arithmetic<T>::value, bool>::type>
 class PiecewiseGrid1D {
@@ -46,7 +47,7 @@ class PiecewiseGrid1D {
   // This is functionally equivalent because grids_ will
   // be initialized to default values
   PORTABLE_INLINE_FUNCTION PiecewiseGrid1D() {}
-  PiecewiseGrid1D(const std::vector<RegularGrid1D<T>> grids) {
+  PiecewiseGrid1D(const std::vector<RegularGrid1D<T, Transform>> grids) {
     NGRIDS_ = grids.size();
     PORTABLE_ALWAYS_REQUIRE(
         NGRIDS_ <= NGRIDSMAX,
@@ -65,8 +66,8 @@ class PiecewiseGrid1D {
       }
     }
   }
-  PiecewiseGrid1D(std::initializer_list<RegularGrid1D<T>> grids)
-      : PiecewiseGrid1D(std::vector<RegularGrid1D<T>>(grids)) {}
+  PiecewiseGrid1D(std::initializer_list<RegularGrid1D<T, Transform>> grids)
+      : PiecewiseGrid1D(std::vector<RegularGrid1D<T, Transform>>(grids)) {}
 
   template <typename F>
   PORTABLE_INLINE_FUNCTION int findGrid(const F &direction) const {
@@ -124,14 +125,14 @@ class PiecewiseGrid1D {
   }
 
   PORTABLE_INLINE_FUNCTION
-  bool operator==(const PiecewiseGrid1D<T, NGRIDSMAX> &other) const {
+  bool operator==(const PiecewiseGrid1D<T, NGRIDSMAX, Transform> &other) const {
     for (int ig = 0; ig < NGRIDS_; ++ig) {
       if (grids_[ig] != other.grids_[ig]) return false;
     }
     return true;
   }
   PORTABLE_INLINE_FUNCTION
-  bool operator!=(const PiecewiseGrid1D<T, NGRIDSMAX> &other) const {
+  bool operator!=(const PiecewiseGrid1D<T, NGRIDSMAX, Transform> &other) const {
     return !(*this == other);
   }
   PORTABLE_INLINE_FUNCTION T min() const { return grids_[0].min(); }
@@ -225,7 +226,7 @@ class PiecewiseGrid1D {
            SP5::H1D::GRID_FORMAT[1];
   }
 
-  RegularGrid1D<T> grids_[NGRIDSMAX];
+  RegularGrid1D<T, Transform> grids_[NGRIDSMAX];
   int pointTotals_[NGRIDSMAX];
   int NGRIDS_;
 };
diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp
index d69076d65..fa3772d41 100644
--- a/spiner/regular_grid_1d.hpp
+++ b/spiner/regular_grid_1d.hpp
@@ -31,6 +31,7 @@
 #include "ports-of-call/portable_errors.hpp"
 #include "sp5.hpp"
 #include "spiner_types.hpp"
+#include "transformations.hpp"
 
 namespace Spiner {
 
@@ -44,7 +45,7 @@ struct weights_t {
   }
 };
 
-template <typename T = Real,
+template <typename T = Real, typename Transform = TransformLinear,
           typename std::enable_if<std::is_arithmetic<T>::value, bool>::type =
               true>
 class RegularGrid1D {
@@ -55,34 +56,23 @@ class RegularGrid1D {
 
   // Constructors
   PORTABLE_INLINE_FUNCTION RegularGrid1D()
-      : min_(rNaN), max_(rNaN), dx_(rNaN), idx_(rNaN), N_(iNaN) {}
-  PORTABLE_INLINE_FUNCTION RegularGrid1D(T min, T max, size_t N)
-      : min_(min), max_(max), dx_((max - min) / ((T)(N - 1))), idx_(1 / dx_),
-        N_(N) {
-    PORTABLE_ALWAYS_REQUIRE(min_ < max_ && N_ > 0, "Valid grid");
-  }
-
-  // Forces x in the interval
-  PORTABLE_INLINE_FUNCTION int bound(int ix) const {
-#ifndef SPINER_DISABLE_BOUNDS_CHECKS
-    if (ix < 0) ix = 0;
-    if (ix >= (int)N_ - 1) ix = (int)N_ - 2; // Ensures ix+1 exists
-#endif
-    return ix;
-  }
-
-  // Gets real value at index
-  PORTABLE_INLINE_FUNCTION T x(const int i) const { return i * dx_ + min_; }
-  PORTABLE_INLINE_FUNCTION int index(const T x) const {
-    return bound(idx_ * (x - min_));
+      : umin_(rNaN), umax_(rNaN), du_(rNaN), inv_du_(rNaN), N_(iNaN) {}
+  PORTABLE_INLINE_FUNCTION RegularGrid1D(T xmin, T xmax, size_t N)
+      : xmin_(xmin), xmax_(xmax), umin_(Transform::forward(xmin)),
+        umax_(Transform::forward(xmax)),
+        du_((umax_ - umin_) / static_cast<T>(N - 1)), inv_du_(1 / du_), N_(N) {
+    // A transform could be monotonically decreasing, so there's no guarantee
+    // that umin_ < umax_
+    PORTABLE_ALWAYS_REQUIRE(xmin_ < xmax_ && N_ > 0, "Valid grid");
   }
 
   // Returns closest index and weights for interpolation
   PORTABLE_INLINE_FUNCTION void weights(const T &x, int &ix,
                                         weights_t<T> &w) const {
-    ix = index(x);
-    const auto floor = static_cast<T>(ix) * dx_ + min_;
-    w[1] = idx_ * (x - floor);
+    const T u = Transform::forward(x);
+    ix = index_u(u);
+    const auto floor = static_cast<T>(ix) * du_ + umin_;
+    w[1] = inv_du_ * (u - floor);
     w[0] = (1. - w[1]);
   }
 
@@ -95,26 +85,50 @@ class RegularGrid1D {
     return w[0] * A(ix) + w[1] * A(ix + 1);
   }
 
-  // utitilies
+  // (in)equality comparison
   PORTABLE_INLINE_FUNCTION bool
-  operator==(const RegularGrid1D<T> &other) const {
-    return (min_ == other.min_ && max_ == other.max_ && dx_ == other.dx_ &&
-            idx_ == other.idx_ && N_ == other.N_);
+  operator==(const RegularGrid1D<T, Transform> &other) const {
+    return (umin_ == other.umin_ && umax_ == other.umax_ && du_ == other.du_ &&
+            N_ == other.N_);
   }
   PORTABLE_INLINE_FUNCTION bool
-  operator!=(const RegularGrid1D<T> &other) const {
+  operator!=(const RegularGrid1D<T, Transform> &other) const {
     return !(*this == other);
   }
-  PORTABLE_INLINE_FUNCTION T min() const { return min_; }
-  PORTABLE_INLINE_FUNCTION T max() const { return max_; }
-  PORTABLE_INLINE_FUNCTION T dx() const { return dx_; }
+
+  // queries
+  PORTABLE_INLINE_FUNCTION T min() const { return xmin_; }
+  PORTABLE_INLINE_FUNCTION T max() const { return xmax_; }
   PORTABLE_INLINE_FUNCTION size_t nPoints() const { return N_; }
   PORTABLE_INLINE_FUNCTION bool isnan() const {
-    return (std::isnan(min_) || std::isnan(max_) || std::isnan(dx_) ||
-            std::isnan(idx_) || std::isnan((T)N_));
+    return (std::isnan(xmin_) || std::isnan(xmax_) || std::isnan(umin_) ||
+            std::isnan(umax_) || std::isnan(du_) || std::isnan(inv_du_) ||
+            std::isnan((T)N_));
   }
   PORTABLE_INLINE_FUNCTION bool isWellFormed() const { return !isnan(); }
+  // TODO: min() and x(0) won't necessarily match.
+  //       max() and x(nPoints-1) won't necessarily match.
+  //       Possible fixes:
+  //       -- If (i == 0) return xmin_.  Introduces two "if" statements every
+  //          time you look up at x value.
+  //       -- Allow min() and x(0) to differ.  This could be confusing and
+  //          cause issues in corner cases.
+  //       -- Allow min() to be different from xmin passed in by the user,
+  //          which may be unexpected (and possibly undesirable) behavior.
+  //       -- Apply an additional linear transformation calibrated so that the
+  //          bounds are exact.  This may muck with any calibrations of the
+  //          transform itself.
+  //       -- Carry an additional array of x values to be used here.  This
+  //          introduces a lot of memory, when no memory was used before.
+  // Translate between x coordinate and index
+  PORTABLE_INLINE_FUNCTION T x(const int i) const {
+    return Transform::reverse(u(i));
+  }
+  PORTABLE_INLINE_FUNCTION int index(const T x) const {
+    return index_u(Transform::forward(x));
+  }
 
+  // HDF
 #ifdef SPINER_USE_HDF
   inline herr_t saveHDF(hid_t loc, const std::string &name) const {
     static_assert(
@@ -123,7 +137,7 @@ class RegularGrid1D {
     auto H5T_T =
         std::is_same<T, double>::value ? H5T_NATIVE_DOUBLE : H5T_NATIVE_FLOAT;
     herr_t status;
-    T range[] = {min_, max_, dx_};
+    T range[] = {umin_, umax_, du_};
     hsize_t range_dims[] = {3};
     int n = static_cast<int>(N_);
     status = H5LTmake_dataset(loc, name.c_str(), SP5::RG1D::RANGE_RANK,
@@ -144,10 +158,10 @@ class RegularGrid1D {
     T range[3];
     int n;
     status = H5LTread_dataset(loc, name.c_str(), H5T_T, range);
-    min_ = range[0];
-    max_ = range[1];
-    dx_ = range[2];
-    idx_ = 1. / dx_;
+    umin_ = range[0];
+    umax_ = range[1];
+    du_ = range[2];
+    inv_du_ = 1. / du_;
     status += H5LTget_attribute_int(loc, name.c_str(), SP5::RG1D::N, &n);
     N_ = n;
     return status;
@@ -155,8 +169,24 @@ class RegularGrid1D {
 #endif
 
  private:
-  T min_, max_;
-  T dx_, idx_;
+  // Forces x in the interval
+  PORTABLE_INLINE_FUNCTION int bound(int ix) const {
+#ifndef SPINER_DISABLE_BOUNDS_CHECKS
+    if (ix < 0) ix = 0;
+    if (ix >= (int)N_ - 1) ix = (int)N_ - 2; // Ensures ix+1 exists
+#endif
+    return ix;
+  }
+
+  // Translate between u (transformed variable) coordinate and index
+  PORTABLE_INLINE_FUNCTION T u(const int i) const { return i * du_ + umin_; }
+  PORTABLE_INLINE_FUNCTION int index_u(const T u) const {
+    return bound(inv_du_ * (u - umin_));
+  }
+
+  T xmin_, xmax_;
+  T umin_, umax_;
+  T du_, inv_du_;
   size_t N_;
 };
 
diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp
new file mode 100644
index 000000000..f491c3ba0
--- /dev/null
+++ b/spiner/transformations.hpp
@@ -0,0 +1,90 @@
+#ifndef _SPINER_TRANSFORM_HPP_
+#define _SPINER_TRANSFORM_HPP_
+//======================================================================
+// © (or copyright) 2019-2021. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+//======================================================================
+
+#include "ports-of-call/portability.hpp"
+
+#include <cmath>
+#include <limits>
+
+namespace Spiner {
+
+// Note on notation:
+// -- "real" space is called x
+// -- "transformed" space is called u
+// -- u = forward(x)
+// -- x = reverse(u)
+
+// linear transformation (aka no-op): y = x
+struct TransformLinear {
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) {
+    return x;
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) {
+    return u;
+  }
+};
+
+// logarithmic transformation: y = log(x + small)
+struct TransformLogarithmic {
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) {
+    constexpr T eps = eps_f<T>();
+    return std::log(std::abs(x) + eps);
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) {
+    constexpr T eps = eps_r<T>();
+    return std::exp(u) - eps;
+  }
+
+ private:
+  // When possible, we use asymetric epsilon values to ensure that
+  // reverse(forward(0)) is exact.  In general, a performant calculation is
+  // more important than getting this value exactly correct, so we require that
+  // our epsilon values be constexpr.
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T eps_f() {
+    return std::numeric_limits<T>::min();
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T eps_r() {
+    // If std::exp and std::log were constexpr, we could explicitly calculate
+    // the right epsilon value as a constexpr value:
+    //     return std::exp(forward<T>(static_cast<T>(0)));
+    // Unfortunately, we have to wait for C++26 for constexpr math.  We
+    // hard-code certain known values, but default to a symmetric epsilon if
+    // that's the best that we can do.
+    if constexpr (std::is_same_v<std::decay_t<T>, float> &&
+                  std::numeric_limits<T>::is_iec559) {
+      return 1.175490707446280263444352135525329e-38f;
+    } else if constexpr (std::is_same_v<std::decay_t<T>, double> &&
+                         std::numeric_limits<T>::is_iec559) {
+      return 2.225073858507262647230317031903882e-308;
+    } else {
+      return eps_f<T>();
+    }
+  }
+};
+
+// TODO: log_NQT and arcsinh_NQT, but these require adding a dependency on
+//       https://github.com/lanl/not-quite-transcendental.  I may leave this for
+//       Jonah ;-)
+
+} // namespace Spiner
+#endif // _SPINER_TRANSFORM_HPP_
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 8e019ce3e..d843689d3 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -57,7 +57,7 @@ target_link_libraries(spiner_test
     ${spinerTest_POPULATED_TARGETS}
 )
 # add unit tests
-add_executable(test.bin test.cpp)
+add_executable(test.bin test.cpp regular_grid_1d.cpp transformations.cpp)
 target_link_libraries(test.bin
   PRIVATE
     spiner::spiner
diff --git a/test/benchmark.cpp b/test/benchmark.cpp
index 07a7195dd..e59bcbfb8 100644
--- a/test/benchmark.cpp
+++ b/test/benchmark.cpp
@@ -74,8 +74,10 @@ int main(int argc, char *argv[]) {
 
     RegularGrid1D gcoarse(xmin, xmax, ncoarse);
     std::vector<RegularGrid1D> gfine;
+    std::vector<Real> dxfine;
     for (const auto &n : nfine) {
       gfine.push_back(RegularGrid1D(xmin, xmax, n));
+      dxfine.push_back((xmax - xmin) / n);
     }
 
     std::cout << "# ncoarse = " << ncoarse << std::endl;
@@ -98,7 +100,8 @@ int main(int argc, char *argv[]) {
     for (int ifine = 0; ifine < nfine.size(); ++ifine) {
       auto n = nfine[ifine];
       auto g = gfine[ifine];
-      Real d3x = g.dx() * g.dx() * g.dx();
+      auto dx = dxfine[ifine];
+      Real d3x = dx * dx * dx;
       Real L2_error;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::fence();
diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp
new file mode 100644
index 000000000..240519160
--- /dev/null
+++ b/test/regular_grid_1d.cpp
@@ -0,0 +1,217 @@
+// © (or copyright) 2019-2021. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
+#include <ports-of-call/portability.hpp>
+#include <spiner/regular_grid_1d.hpp>
+
+#include <catch2/catch_session.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include <cmath>
+
+using Catch::Matchers::WithinRel;
+
+// test transform: reverse(forward(xlo)) and reverse(forward(xhi)) don't map
+// 100% correctly and end up inside the bounds of [xlo, xhi], so long as xlo
+// and xhi bracket xref = 0.5
+struct TxNarrow {
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) {
+    constexpr T xref = 0.5;
+    constexpr T eps = 10 * std::numeric_limits<T>::epsilon();
+    return (static_cast<T>(1) - eps) * (x - xref) + xref;
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) {
+    constexpr T uref = 0.5;
+    constexpr T eps = 10 * std::numeric_limits<T>::epsilon();
+    return (static_cast<T>(1) - eps) * (u - uref) + uref;
+  }
+};
+
+// test transform: reverse(forward(xlo)) and reverse(forward(xhi)) don't map
+// 100% correctly and end up outside the bounds of [xlo, xhi], so long as xlo
+// and xhi bracket xref = 0.5
+struct TxExpand {
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) {
+    constexpr T xref = 0.5;
+    constexpr T eps = 10 * std::numeric_limits<T>::epsilon();
+    return (static_cast<T>(1) + eps) * (x - xref) + xref;
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) {
+    constexpr T uref = 0.5;
+    constexpr T eps = 10 * std::numeric_limits<T>::epsilon();
+    return (static_cast<T>(1) + eps) * (u - uref) + uref;
+  }
+};
+
+// test transform: monotonically decreasing instead of increasing
+struct TxDecrease {
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) {
+    return 1 - x;
+  }
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) {
+    return 1 - u;
+  }
+};
+
+TEST_CASE("RegularGrid1D", "[RegularGrid1D]") {
+  SECTION("A regular grid 1d emits appropriate metadata") {
+    constexpr Real min = -1;
+    constexpr Real max = 1;
+    constexpr size_t N = 10;
+    Spiner::RegularGrid1D<Real> g(min, max, N);
+    REQUIRE(g.min() == min);
+    REQUIRE(g.max() == max);
+    REQUIRE(g.nPoints() == N);
+  }
+}
+
+TEST_CASE("RegularGrid1D with production transformations",
+          "[RegularGrid1D][transformations]") {
+  constexpr double min = 1;
+  constexpr double max = 1024;
+  constexpr size_t N = 11;
+  Spiner::RegularGrid1D<double, Spiner::TransformLinear> glin(min, max, N);
+  Spiner::RegularGrid1D<double, Spiner::TransformLogarithmic> glog(min, max, N);
+  REQUIRE(glin.min() == min);
+  REQUIRE(glog.min() == min);
+  REQUIRE(glin.max() == max);
+  REQUIRE(glog.max() == max);
+  REQUIRE(glin.nPoints() == N);
+  REQUIRE(glog.nPoints() == N);
+
+  // Check all fixed points (lin)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = 102.3 * double(n) + 1;
+    CHECK_THAT(glin.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check all fixed points (log)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = std::pow(double(2), n);
+    CHECK_THAT(glog.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check weights (lin)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double xlin = 0.5 * (glin.x(n - 1) + glin.x(n));
+    int ilin;
+    Spiner::weights_t<double> wlin;
+    glin.weights(xlin, ilin, wlin);
+    CHECK(ilin == n - 1);
+    CHECK_THAT(wlin.first, WithinRel(0.5));
+    CHECK_THAT(wlin.second, WithinRel(0.5));
+  }
+
+  // Check weights (log)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double xlog = std::sqrt(glog.x(n - 1) * glog.x(n));
+    int ilog;
+    Spiner::weights_t<double> wlog;
+    glog.weights(xlog, ilog, wlog);
+    CHECK(ilog == n - 1);
+    CHECK_THAT(wlog.first, WithinRel(0.5));
+    CHECK_THAT(wlog.second, WithinRel(0.5));
+  }
+}
+
+TEST_CASE("RegularGrid1D with test transformations",
+          "[RegularGrid1D][transformations]") {
+  constexpr double min = 0.0;
+  constexpr double max = 1.0;
+  constexpr size_t N = 101;
+  Spiner::RegularGrid1D<double, TxNarrow> gn(min, max, N);
+  Spiner::RegularGrid1D<double, TxExpand> ge(min, max, N);
+  Spiner::RegularGrid1D<double, TxDecrease> gd(min, max, N);
+
+  // Basic tests (narrow)
+  REQUIRE(gn.min() == min);
+  REQUIRE(gn.max() == max);
+  REQUIRE(gn.nPoints() == N);
+  // TODO: Do we want the bounds to be exact?
+  CHECK(gn.x(0) == min);
+  CHECK(gn.x(N - 1) == max);
+
+  // Basic tests (expand)
+  REQUIRE(ge.min() == min);
+  REQUIRE(ge.max() == max);
+  REQUIRE(ge.nPoints() == N);
+  // TODO: Do we want the bounds to be exact?
+  CHECK(ge.x(0) == min);
+  CHECK(ge.x(N - 1) == max);
+
+  // Basic tests (decrease)
+  REQUIRE(gd.min() == min);
+  REQUIRE(gd.max() == max);
+  REQUIRE(gd.nPoints() == N);
+  // TODO: Do we want the bounds to be exact?
+  CHECK(gd.x(0) == min);
+  CHECK(gd.x(N - 1) == max);
+
+  // Check all fixed points (narrow)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = 0.01 * double(n);
+    CHECK_THAT(gn.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check all fixed points (expand)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = 0.01 * double(n);
+    CHECK_THAT(ge.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check all fixed points (decrease)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = 0.01 * double(n);
+    CHECK_THAT(gd.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check weights (narrow)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double x = 0.5 * (gn.x(n - 1) + gn.x(n));
+    int index;
+    Spiner::weights_t<double> weights;
+    gn.weights(x, index, weights);
+    CHECK(index == n - 1);
+    CHECK_THAT(weights.first, WithinRel(0.5, 1.0e-12));
+    CHECK_THAT(weights.second, WithinRel(0.5, 1.0e-12));
+  }
+
+  // Check weights (expand)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double x = 0.5 * (ge.x(n - 1) + ge.x(n));
+    int index;
+    Spiner::weights_t<double> weights;
+    ge.weights(x, index, weights);
+    CHECK(index == n - 1);
+    CHECK_THAT(weights.first, WithinRel(0.5, 1.0e-12));
+    CHECK_THAT(weights.second, WithinRel(0.5, 1.0e-12));
+  }
+
+  // Check weights (decrease)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double x = 0.5 * (gd.x(n - 1) + gd.x(n));
+    int index;
+    Spiner::weights_t<double> weights;
+    gd.weights(x, index, weights);
+    CHECK(index == n - 1);
+    CHECK_THAT(weights.first, WithinRel(0.5, 1.0e-12));
+    CHECK_THAT(weights.second, WithinRel(0.5, 1.0e-12));
+  }
+}
diff --git a/test/test.cpp b/test/test.cpp
index 30671ff93..a4ddc6eb0 100644
--- a/test/test.cpp
+++ b/test/test.cpp
@@ -52,6 +52,16 @@ PORTABLE_INLINE_FUNCTION Real linearFunction(Real b, Real a, Real z, Real y,
                                              Real x) {
   return x + y + z + a + b;
 }
+PORTABLE_INLINE_FUNCTION double log_lin_func(const double x, const double y) {
+  // log-log: f is linear function in log(x) and log(y)
+  // --> f = a + b * log(x) + c * log(y)
+  return 0.7 * std::log(x) + 0.3 * std::log(y) - 0.5;
+}
+PORTABLE_INLINE_FUNCTION double log_log_func(const double x, const double y) {
+  // log-log: log(f) is linear function in log(x) and log(y)
+  // --> log(f) = a + b * log(x) + c * log(y)
+  return std::exp(0.7 * std::log(x) + 0.3 * std::log(y) - 0.5);
+}
 
 SCENARIO("PortableMDArrays can be allocated from a pointer",
          "[PortableMDArray]") {
@@ -89,19 +99,6 @@ SCENARIO("PortableMDArrays can be allocated from a pointer",
   }
 }
 
-TEST_CASE("RegularGrid1D", "[RegularGrid1D]") {
-  SECTION("A regular grid 1d emits appropriate metadata") {
-    constexpr Real min = -1;
-    constexpr Real max = 1;
-    constexpr size_t N = 10;
-    RegularGrid1D g(min, max, N);
-    REQUIRE(g.min() == min);
-    REQUIRE(g.max() == max);
-    REQUIRE(g.nPoints() == N);
-    REQUIRE(g.dx() == (max - min) / ((Real)(N - 1)));
-  }
-}
-
 TEST_CASE("PiecewiseGrid1D", "[PiecewiseGrid1D]") {
   GIVEN("Some regular grid 1Ds") {
     RegularGrid1D g1(0, 0.25, 3);
@@ -509,6 +506,88 @@ TEST_CASE("DataBox interpolation", "[DataBox]") {
   free(db); // free databox
 }
 
+TEST_CASE("DataBox Interpolation with log-lin transformations", "[DataBox]") {
+  using Transform = Spiner::TransformLogarithmic;
+  using RG1D = Spiner::RegularGrid1D<double, Transform>;
+  using DB = Spiner::DataBox<double, RG1D, Spiner::TransformLinear>;
+
+  constexpr int RANK = 2;
+  const int Npt = 32;
+  const int Nsample = Npt * 16;
+
+  const double xmin = 1;
+  const double xmax = 10;
+
+  DB db(Spiner::AllocationTarget::Device, Npt, Npt);
+  for (int i = 0; i < RANK; i++) {
+    db.setRange(i, xmin, xmax, Npt);
+  }
+
+  portableFor(
+      "fill databox", 0, Npt, 0, Npt,
+      PORTABLE_LAMBDA(const int ix, const int iy) {
+        RG1D grid(xmin, xmax, Npt);
+        const double x = grid.x(ix);
+        const double y = grid.x(iy);
+        db(ix, iy) = log_lin_func(x, y);
+      });
+
+  double error = 0;
+  portableReduce(
+      "interpolate databox", 0, Nsample, 0, Nsample,
+      PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) {
+        RG1D grid(xmin, xmax, Nsample);
+        const double x = grid.x(ix);
+        const double y = grid.x(iy);
+        const double f_true = log_lin_func(x, y);
+        const double difference = db.interpToReal(x, y) - f_true;
+        accumulate += (difference * difference);
+      },
+      error);
+  REQUIRE(error <= EPSTEST);
+}
+
+TEST_CASE("DataBox Interpolation with log-log transformations", "[DataBox]") {
+  using Transform = Spiner::TransformLogarithmic;
+  using RG1D = Spiner::RegularGrid1D<double, Transform>;
+  using DB = Spiner::DataBox<double, RG1D, Transform>;
+
+  constexpr int RANK = 2;
+  const int Npt = 32;
+  const int Nsample = Npt * 16;
+
+  const double xmin = 1;
+  const double xmax = 10;
+
+  DB db(Spiner::AllocationTarget::Device, Npt, Npt);
+  for (int i = 0; i < RANK; i++) {
+    db.setRange(i, xmin, xmax, Npt);
+  }
+
+  portableFor(
+      "fill databox", 0, Npt, 0, Npt,
+      PORTABLE_LAMBDA(const int ix, const int iy) {
+        RG1D grid(xmin, xmax, Npt);
+        const double x = grid.x(ix);
+        const double y = grid.x(iy);
+        db.set_data_value(log_log_func(x, y), ix, iy);
+      });
+
+  double error = 0;
+  portableReduce(
+      "interpolate databox", 0, Nsample, 0, Nsample,
+      PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) {
+        RG1D grid(xmin, xmax, Nsample);
+        const double x = grid.x(ix);
+        const double y = grid.x(iy);
+        const double f_true = log_log_func(x, y);
+        const double difference = db.interpToReal(x, y) - f_true;
+        accumulate += (difference * difference);
+      },
+      error);
+  REQUIRE(error <= EPSTEST);
+}
+
 TEST_CASE("DataBox Interpolation with piecewise grids",
           "[DataBox][PiecewiseGrid1D]") {
   GIVEN("A piecewise grid") {
diff --git a/test/transformations.cpp b/test/transformations.cpp
new file mode 100644
index 000000000..c95259d18
--- /dev/null
+++ b/test/transformations.cpp
@@ -0,0 +1,109 @@
+// © (or copyright) 2019-2021. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
+#include <spiner/transformations.hpp>
+
+#include <catch2/catch_session.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include <cmath>
+
+using Catch::Matchers::WithinAbs;
+using Catch::Matchers::WithinRel;
+using Catch::Matchers::WithinULP;
+
+namespace {
+PORTABLE_FUNCTION double relerr(const double z, const double zz) {
+  if (z == 0) {
+    return std::abs(zz - z);
+  } else {
+    return std::abs(double(1) - zz / z);
+  }
+}
+}; // namespace
+
+TEST_CASE("transformation: linear", "[transformations]") {
+  // This one is almost too simple to meaningfully test, but we can at least
+  // ensure that it compiles and does the trivially-obvious things.
+  using Transform = Spiner::TransformLinear;
+
+  // Test on CPU
+  const double x0 = 3.14159;
+  CHECK(Transform::forward(x0) == x0);
+  CHECK(Transform::reverse(x0) == x0);
+  const float x1 = 2.71828;
+  CHECK(Transform::forward(x1) == x1);
+  CHECK(Transform::reverse(x1) == x1);
+
+  // Test on GPU (or CPU if no GPU available)
+  const int num_threads = 10;
+  double accum = 0;
+  portableReduce(
+      "parallel_reduce", 0, num_threads,
+      PORTABLE_LAMBDA(const int n, double &reduce) {
+        const double x = static_cast<double>(n);
+        const double y = x;
+        const double yy = Transform::forward(x);
+        const double xx = Transform::reverse(y);
+        reduce += 0.5 * (relerr(x, xx) + relerr(y, yy));
+      },
+      accum);
+  CHECK(accum / num_threads == 0);
+}
+
+TEST_CASE("transformation: logarithmic", "[transformations]") {
+  using Transform = Spiner::TransformLogarithmic;
+
+  // Test on CPU
+  auto matcher = [](const double d) {
+    // return WithinULP(d, 500);
+    return WithinRel(d, 1.0e-12);
+  };
+  // Scan across most of the range
+  for (int p = -307; p <= 307; ++p) {
+    const double x = std::pow(double(10), p);
+    const double y = std::log(x);
+    // Basic checks
+    // -- The "exact" calculation (y = log(x)) won't match the transformation
+    //    (y = log(x + epsilon) for very small values of x, so skip checks.
+    if (p >= -295) {
+      CHECK_THAT(Transform::forward(x), matcher(y));
+      CHECK_THAT(Transform::reverse(y), matcher(x));
+    }
+    // Round trip
+    CHECK_THAT(Transform::reverse(Transform::forward(x)), matcher(x));
+  }
+  // Special value
+  // -- Transform::forward(0) will infer the type to be an integer, and you
+  //    will get a WILDLY incorrect answer for the round trip.
+  CHECK(std::isfinite(Transform::forward(0.0)));
+  CHECK(Transform::reverse(Transform::forward(static_cast<double>(0))) == 0);
+  CHECK(Transform::reverse(Transform::forward(static_cast<float>(0))) == 0);
+
+  // Test on GPU (or CPU if no GPU available)
+  const int num_threads = 101;
+  double accum = 0;
+  portableReduce(
+      "parallel_reduce", 0, num_threads,
+      PORTABLE_LAMBDA(const int n, double &reduce) {
+        const int p = n - num_threads / 2;
+        const double x = std::pow(double(10.0), p);
+        const double y = std::log(x);
+        const double yy = Transform::forward(x);
+        const double xx = Transform::reverse(y);
+        reduce += 0.5 * (relerr(x, xx) + relerr(y, yy));
+      },
+      accum);
+  CHECK(accum / num_threads <= 6.0e-12);
+}