diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index ff5380c537..e82c4051f9 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -277,7 +277,7 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
         direction = opposite_direction(direction);
 
     if (P->iso_obj != nullptr &&
-        dynamic_cast<NS_PROJ::operation::CoordinateOperation*>(P->iso_obj.get()) == nullptr ) {
+        !P->iso_obj_is_coordinate_operation ) {
         pj_log(P->ctx, PJ_LOG_ERROR, "Object is not a coordinate operation");
         proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
         return proj_coord_error ();
@@ -326,8 +326,11 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
                 }
                 P->iCurCoordOp = iBest;
             }
-            PJ_COORD res = direction == PJ_FWD ?
-                        pj_fwd4d( coord, alt.pj ) : pj_inv4d( coord, alt.pj );
+            PJ_COORD res = coord;
+            if( direction == PJ_FWD )
+                pj_fwd4d( res, alt.pj );
+            else
+                pj_inv4d( res, alt.pj );
             if( proj_errno(alt.pj) == PROJ_ERR_OTHER_NETWORK_ERROR ) {
                 return proj_coord_error ();
             }
@@ -368,11 +371,12 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
                         P->iCurCoordOp = i;
                     }
                     if( direction == PJ_FWD ) {
-                        return pj_fwd4d( coord, alt.pj );
+                        pj_fwd4d( coord, alt.pj );
                     }
                     else {
-                        return pj_inv4d( coord, alt.pj );
+                        pj_inv4d( coord, alt.pj );
                     }
+                    return coord;
                 }
             }
         }
@@ -383,9 +387,10 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
 
     P->iCurCoordOp = 0; // dummy value, to be used by proj_trans_get_last_used_operation()
     if (direction == PJ_FWD)
-        return pj_fwd4d (coord, P);
+        pj_fwd4d (coord, P);
     else
-        return pj_inv4d (coord, P);
+        pj_inv4d (coord, P);
+    return coord;
 }
 
 /*****************************************************************************/
diff --git a/src/conversions/axisswap.cpp b/src/conversions/axisswap.cpp
index e944b7dab9..c19e2233a6 100644
--- a/src/conversions/axisswap.cpp
+++ b/src/conversions/axisswap.cpp
@@ -57,6 +57,8 @@ It is only necessary to specify the axes that are affected by the swap
 #include <stdlib.h>
 #include <string.h>
 
+#include <algorithm>
+
 #include "proj.h"
 #include "proj_internal.h"
 
@@ -75,20 +77,14 @@ static int sign(int x) {
 
 static PJ_XY forward_2d(PJ_LP lp, PJ *P) {
     struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
-    unsigned int i;
-    PJ_COORD out, in;
-
-    in.v[0] = lp.lam;
-    in.v[1] = lp.phi;
-    out = proj_coord_error();
+    PJ_XY xy;
 
-    for (i=0; i<2; i++)
-        out.v[i] = in.v[Q->axis[i]] * Q->sign[i];
-
-    return out.xy;
+    double in[2] = {lp.lam, lp.phi};
+    xy.x = in[Q->axis[0]] * Q->sign[0];
+    xy.y = in[Q->axis[1]] * Q->sign[1];
+    return xy;
 }
 
-
 static PJ_LP reverse_2d(PJ_XY xy, PJ *P) {
     struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
     unsigned int i;
@@ -104,7 +100,6 @@ static PJ_LP reverse_2d(PJ_XY xy, PJ *P) {
     return out.lp;
 }
 
-
 static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
     unsigned int i;
@@ -137,32 +132,31 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     return out.lpz;
 }
 
+static void swap_xy_4d(PJ_COORD& coo, PJ *) {
+    std::swap(coo.xyzt.x, coo.xyzt.y);
+}
+
 
-static PJ_COORD forward_4d(PJ_COORD coo, PJ *P) {
+static void forward_4d(PJ_COORD& coo, PJ *P) {
     struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
     unsigned int i;
     PJ_COORD out;
 
-    out = proj_coord_error();
-
     for (i=0; i<4; i++)
         out.v[i] = coo.v[Q->axis[i]] * Q->sign[i];
-
-    return out;
+    coo = out;
 }
 
 
-static PJ_COORD reverse_4d(PJ_COORD coo, PJ *P) {
+static void reverse_4d(PJ_COORD& coo, PJ *P) {
     struct pj_opaque *Q = (struct pj_opaque *) P->opaque;
     unsigned int i;
     PJ_COORD out;
 
-    out = proj_coord_error();
-
     for (i=0; i<4; i++)
         out.v[Q->axis[i]] = coo.v[i] * Q->sign[i];
 
-    return out;
+    coo = out;
 }
 
 
@@ -278,9 +272,18 @@ PJ *CONVERSION(axisswap,0) {
         P->fwd3d  = forward_3d;
         P->inv3d  = reverse_3d;
     }
-    if (n == 2 && Q->axis[0] < 2 && Q->axis[1] < 2) {
-        P->fwd    = forward_2d;
-        P->inv    = reverse_2d;
+    if (n == 2) {
+        if( Q->axis[0] == 1 && Q->sign[0] == 1 &&
+            Q->axis[1] == 0 && Q->sign[1] == 1 )
+        {
+            P->fwd4d = swap_xy_4d;
+            P->inv4d = swap_xy_4d;
+        }
+        else if( Q->axis[0] < 2 && Q->axis[1] < 2 )
+        {
+            P->fwd    = forward_2d;
+            P->inv    = reverse_2d;
+        }
     }
 
 
diff --git a/src/conversions/geoc.cpp b/src/conversions/geoc.cpp
index b5a6ad7718..19c9b72153 100644
--- a/src/conversions/geoc.cpp
+++ b/src/conversions/geoc.cpp
@@ -36,13 +36,13 @@
 PROJ_HEAD(geoc, "Geocentric Latitude");
 
 /* Geographical to geocentric */
-static PJ_COORD forward(PJ_COORD coo, PJ *P) {
-    return pj_geocentric_latitude (P, PJ_FWD, coo);
+static void forward(PJ_COORD& coo, PJ *P) {
+    coo = pj_geocentric_latitude (P, PJ_FWD, coo);
 }
 
 /* Geocentric to geographical */
-static PJ_COORD inverse(PJ_COORD coo, PJ *P) {
-    return pj_geocentric_latitude (P, PJ_INV, coo);
+static void inverse(PJ_COORD& coo, PJ *P) {
+    coo = pj_geocentric_latitude (P, PJ_INV, coo);
 }
 
 
diff --git a/src/conversions/noop.cpp b/src/conversions/noop.cpp
index f8d1fa5bab..f5e732be8d 100644
--- a/src/conversions/noop.cpp
+++ b/src/conversions/noop.cpp
@@ -4,9 +4,7 @@
 
 PROJ_HEAD(noop,    "No operation");
 
-static PJ_COORD noop(PJ_COORD coord, PJ *P) {
-    (void) P;
-    return coord;
+static void noop(PJ_COORD& , PJ *) {
 }
 
 PJ *CONVERSION(noop, 0) {
diff --git a/src/conversions/set.cpp b/src/conversions/set.cpp
index 1ed0dfc1da..0d22a5ed1d 100644
--- a/src/conversions/set.cpp
+++ b/src/conversions/set.cpp
@@ -19,7 +19,7 @@ struct Set {
 };
 } // anonymous namespace
 
-static PJ_COORD set_fwd_inv(PJ_COORD point, PJ *P) {
+static void set_fwd_inv(PJ_COORD& point, PJ *P) {
 
     struct Set *set = static_cast<struct Set*>(P->opaque);
 
@@ -31,8 +31,6 @@ static PJ_COORD set_fwd_inv(PJ_COORD point, PJ *P) {
         point.v[2] = set->v3_val;
     if (set->v4)
         point.v[3] = set->v4_val;
-
-    return point;
 }
 
 PJ *OPERATION(set, 0) {
diff --git a/src/conversions/topocentric.cpp b/src/conversions/topocentric.cpp
index 2fee149d52..97db380118 100644
--- a/src/conversions/topocentric.cpp
+++ b/src/conversions/topocentric.cpp
@@ -48,28 +48,27 @@ struct pj_opaque {
 } // anonymous namespace
 
 // Convert from geocentric to topocentric
-static PJ_COORD topocentric_fwd(PJ_COORD in, PJ * P)
+static void topocentric_fwd(PJ_COORD& coo, PJ * P)
 {
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    PJ_COORD out;
-    const double dX = in.xyz.x - Q->X0;
-    const double dY = in.xyz.y - Q->Y0;
-    const double dZ = in.xyz.z - Q->Z0;
-    out.xyz.x = -dX * Q->sinlam0              + dY * Q->coslam0;
-    out.xyz.y = -dX * Q->sinphi0 * Q->coslam0 - dY * Q->sinphi0 * Q->sinlam0 + dZ * Q->cosphi0;
-    out.xyz.z =  dX * Q->cosphi0 * Q->coslam0 + dY * Q->cosphi0 * Q->sinlam0 + dZ * Q->sinphi0;
-    return out;
+    const double dX = coo.xyz.x - Q->X0;
+    const double dY = coo.xyz.y - Q->Y0;
+    const double dZ = coo.xyz.z - Q->Z0;
+    coo.xyz.x = -dX * Q->sinlam0              + dY * Q->coslam0;
+    coo.xyz.y = -dX * Q->sinphi0 * Q->coslam0 - dY * Q->sinphi0 * Q->sinlam0 + dZ * Q->cosphi0;
+    coo.xyz.z =  dX * Q->cosphi0 * Q->coslam0 + dY * Q->cosphi0 * Q->sinlam0 + dZ * Q->sinphi0;
 }
 
 // Convert from topocentric to geocentric
-static PJ_COORD topocentric_inv(PJ_COORD in, PJ * P)
+static void topocentric_inv(PJ_COORD& coo, PJ * P)
 {
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    PJ_COORD out;
-    out.xyz.x = Q->X0 - in.xyz.x * Q->sinlam0 - in.xyz.y * Q->sinphi0 * Q->coslam0 + in.xyz.z * Q->cosphi0 * Q->coslam0;
-    out.xyz.y = Q->Y0 + in.xyz.x * Q->coslam0 - in.xyz.y * Q->sinphi0 * Q->sinlam0 + in.xyz.z * Q->cosphi0 * Q->sinlam0;
-    out.xyz.z = Q->Z0                         + in.xyz.y * Q->cosphi0              + in.xyz.z * Q->sinphi0;
-    return out;
+    const double x = coo.xyz.x;
+    const double y = coo.xyz.y;
+    const double z = coo.xyz.z;
+    coo.xyz.x = Q->X0 - x * Q->sinlam0 - y * Q->sinphi0 * Q->coslam0 + z * Q->cosphi0 * Q->coslam0;
+    coo.xyz.y = Q->Y0 + x * Q->coslam0 - y * Q->sinphi0 * Q->sinlam0 + z * Q->cosphi0 * Q->sinlam0;
+    coo.xyz.z = Q->Z0                         + y * Q->cosphi0              + z * Q->sinphi0;
 }
 
 
diff --git a/src/conversions/unitconvert.cpp b/src/conversions/unitconvert.cpp
index 3be0faff24..417d1ae4f0 100644
--- a/src/conversions/unitconvert.cpp
+++ b/src/conversions/unitconvert.cpp
@@ -351,42 +351,36 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
 
 
 /***********************************************************************/
-static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+static void forward_4d(PJ_COORD& coo, PJ *P) {
 /************************************************************************
     Forward conversion of time units
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_COORD out = obs;
 
     /* delegate unit conversion of physical dimensions to the 3D function */
-    out.xyz = forward_3d(obs.lpz, P);
+    coo.xyz = forward_3d(coo.lpz, P);
 
     if (Q->t_in_id >= 0)
-        out.xyzt.t = time_units[Q->t_in_id].t_in( obs.xyzt.t );
+        coo.xyzt.t = time_units[Q->t_in_id].t_in( coo.xyzt.t );
     if (Q->t_out_id >= 0)
-        out.xyzt.t = time_units[Q->t_out_id].t_out( out.xyzt.t );
-
-    return out;
+        coo.xyzt.t = time_units[Q->t_out_id].t_out( coo.xyzt.t );
 }
 
 
 /***********************************************************************/
-static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+static void reverse_4d(PJ_COORD& coo, PJ *P) {
 /************************************************************************
     Reverse conversion of time units
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_COORD out = obs;
 
     /* delegate unit conversion of physical dimensions to the 3D function */
-    out.lpz = reverse_3d(obs.xyz, P);
+    coo.lpz = reverse_3d(coo.xyz, P);
 
     if (Q->t_out_id >= 0)
-        out.xyzt.t = time_units[Q->t_out_id].t_in( obs.xyzt.t );
+        coo.xyzt.t = time_units[Q->t_out_id].t_in( coo.xyzt.t );
     if (Q->t_in_id >= 0)
-        out.xyzt.t = time_units[Q->t_in_id].t_out( out.xyzt.t );
-
-    return out;
+        coo.xyzt.t = time_units[Q->t_in_id].t_out( coo.xyzt.t );
 }
 
 /***********************************************************************/
diff --git a/src/fwd.cpp b/src/fwd.cpp
index 36e2ce81c2..fd893d3ded 100644
--- a/src/fwd.cpp
+++ b/src/fwd.cpp
@@ -203,7 +203,7 @@ PJ_XY pj_fwd(PJ_LP lp, PJ *P) {
         coo.xyz = xyz;
     }
     else if (P->fwd4d)
-        coo = P->fwd4d (coo, P);
+        P->fwd4d (coo, P);
     else {
         proj_errno_set (P, PROJ_ERR_OTHER_NO_INVERSE_OP);
         return proj_coord_error ().xy;
@@ -238,7 +238,7 @@ PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
         coo.xyz = xyz;
     }
     else if (P->fwd4d)
-        coo = P->fwd4d (coo, P);
+        P->fwd4d (coo, P);
     else if (P->fwd)
     {
         const auto xy = P->fwd (coo.lp, P);
@@ -259,7 +259,7 @@ PJ_XYZ pj_fwd3d(PJ_LPZ lpz, PJ *P) {
 
 
 
-PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
+bool pj_fwd4d (PJ_COORD& coo, PJ *P) {
 
     const int last_errno = P->ctx->last_errno;
     P->ctx->last_errno = 0;
@@ -267,11 +267,14 @@ PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
     if (!P->skip_fwd_prepare)
         fwd_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
-        return proj_coord_error ();
+    {
+        coo = proj_coord_error ();
+        return false;
+    }
 
     /* Call the highest dimensional converter available */
     if (P->fwd4d)
-        coo = P->fwd4d (coo, P);
+        P->fwd4d (coo, P);
     else if (P->fwd3d)
     {
         const auto xyz = P->fwd3d (coo.lpz, P);
@@ -284,13 +287,24 @@ PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
     }
     else {
         proj_errno_set (P, PROJ_ERR_OTHER_NO_INVERSE_OP);
-        return proj_coord_error ();
+        coo = proj_coord_error ();
+        return false;
     }
     if (HUGE_VAL==coo.v[0])
-        return proj_coord_error ();
+    {
+        coo = proj_coord_error ();
+        return false;
+    }
 
     if (!P->skip_fwd_finalize)
         fwd_finalize (P, coo);
 
-    return error_or_coord(P, coo, last_errno);
+    if (P->ctx->last_errno)
+    {
+        coo = proj_coord_error();
+        return false;
+    }
+
+    P->ctx->last_errno = last_errno;
+    return true;
 }
diff --git a/src/inv.cpp b/src/inv.cpp
index 92b3c1d62d..14d6233441 100644
--- a/src/inv.cpp
+++ b/src/inv.cpp
@@ -161,7 +161,7 @@ PJ_LP pj_inv(PJ_XY xy, PJ *P) {
         coo.lpz = lpz;
     }
     else if (P->inv4d)
-        coo = P->inv4d (coo, P);
+        P->inv4d (coo, P);
     else {
         proj_errno_set (P, PROJ_ERR_OTHER_NO_INVERSE_OP);
         return proj_coord_error ().lp;
@@ -196,7 +196,7 @@ PJ_LPZ pj_inv3d (PJ_XYZ xyz, PJ *P) {
         coo.lpz = lpz;
     }
     else if (P->inv4d)
-        coo = P->inv4d (coo, P);
+        P->inv4d (coo, P);
     else if (P->inv)
     {
         const auto lp = P->inv (coo.xy, P);
@@ -217,7 +217,7 @@ PJ_LPZ pj_inv3d (PJ_XYZ xyz, PJ *P) {
 
 
 
-PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
+bool pj_inv4d (PJ_COORD& coo, PJ *P) {
 
     const int last_errno = P->ctx->last_errno;
     P->ctx->last_errno = 0;
@@ -225,11 +225,14 @@ PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
     if (!P->skip_inv_prepare)
         inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
-        return proj_coord_error ();
+    {
+        coo = proj_coord_error ();
+        return false;
+    }
 
     /* Call the highest dimensional converter available */
     if (P->inv4d)
-        coo = P->inv4d (coo, P);
+        P->inv4d (coo, P);
     else if (P->inv3d)
     {
         const auto lpz = P->inv3d (coo.xyz, P);
@@ -242,13 +245,24 @@ PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
     }
     else {
         proj_errno_set (P, PROJ_ERR_OTHER_NO_INVERSE_OP);
-        return proj_coord_error ();
+        coo = proj_coord_error ();
+        return false;
     }
     if (HUGE_VAL==coo.v[0])
-        return proj_coord_error ();
+    {
+        coo = proj_coord_error ();
+        return false;
+    }
 
     if (!P->skip_inv_finalize)
         inv_finalize (P, coo);
 
-    return error_or_coord(P, coo, last_errno);
+    if (P->ctx->last_errno)
+    {
+        coo = proj_coord_error();
+        return false;
+    }
+
+    P->ctx->last_errno = last_errno;
+    return true;
 }
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index 117474d09b..7c835c15ab 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -202,6 +202,7 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) {
             ctx->defer_grid_opening = false;
             if (pj) {
                 pj->iso_obj = objIn;
+                pj->iso_obj_is_coordinate_operation = true;
                 return pj;
             }
         } catch (const std::exception &) {
@@ -214,6 +215,7 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) {
         pj->ctx = ctx;
         pj->descr = "ISO-19111 object";
         pj->iso_obj = objIn;
+        pj->iso_obj_is_coordinate_operation = coordop != nullptr;
         try {
             auto crs = dynamic_cast<const CRS *>(objIn.get());
             if (crs) {
diff --git a/src/pipeline.cpp b/src/pipeline.cpp
index e8dc319a53..8351e7c871 100644
--- a/src/pipeline.cpp
+++ b/src/pipeline.cpp
@@ -147,8 +147,8 @@ struct PushPop {
 } // anonymous namespace
 
 
-static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P);
-static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P);
+static void pipeline_forward_4d (PJ_COORD& point, PJ *P);
+static void pipeline_reverse_4d (PJ_COORD& point, PJ *P);
 static PJ_XYZ    pipeline_forward_3d (PJ_LPZ lpz, PJ *P);
 static PJ_LPZ    pipeline_reverse_3d (PJ_XYZ xyz, PJ *P);
 static PJ_XY     pipeline_forward (PJ_LP lp, PJ *P);
@@ -162,24 +162,25 @@ static void pipeline_reassign_context( PJ* P, PJ_CONTEXT* ctx )
 }
 
 
-static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P) {
+static void pipeline_forward_4d (PJ_COORD& point, PJ *P) {
     auto pipeline = static_cast<struct Pipeline*>(P->opaque);
     for( auto& step: pipeline->steps )
     {
         if( !step.omit_fwd )
         {
-            point = proj_trans (step.pj, PJ_FWD, point);
+            if (!step.pj->inverted)
+                pj_fwd4d (point, step.pj);
+            else
+                pj_inv4d (point, step.pj);
             if( point.xyzt.x == HUGE_VAL ) {
                 break;
             }
         }
     }
-
-    return point;
 }
 
 
-static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P) {
+static void pipeline_reverse_4d (PJ_COORD& point, PJ *P) {
     auto pipeline = static_cast<struct Pipeline*>(P->opaque);
     for( auto iterStep = pipeline->steps.rbegin();
               iterStep != pipeline->steps.rend(); ++iterStep )
@@ -187,14 +188,15 @@ static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P) {
         const auto& step = *iterStep;
         if( !step.omit_inv )
         {
-            point = proj_trans (step.pj, PJ_INV, point);
+            if (step.pj->inverted)
+                pj_fwd4d (point, step.pj);
+            else
+                pj_inv4d (point, step.pj);
             if( point.xyzt.x == HUGE_VAL ) {
                 break;
             }
         }
     }
-
-    return point;
 }
 
 
@@ -619,9 +621,9 @@ PJ *OPERATION(pipeline,0) {
     return P;
 }
 
-static PJ_COORD push(PJ_COORD point, PJ *P) {
+static void push(PJ_COORD& point, PJ *P) {
     if (P->parent == nullptr)
-        return point;
+        return;
 
     struct Pipeline *pipeline = static_cast<struct Pipeline*>(P->parent->opaque);
     struct PushPop *pushpop = static_cast<struct PushPop*>(P->opaque);
@@ -634,13 +636,11 @@ static PJ_COORD push(PJ_COORD point, PJ *P) {
         pipeline->stack[2].push(point.v[2]);
     if (pushpop->v4)
         pipeline->stack[3].push(point.v[3]);
-
-    return point;
 }
 
-static PJ_COORD pop(PJ_COORD point, PJ *P) {
+static void pop(PJ_COORD& point, PJ *P) {
     if (P->parent == nullptr)
-        return point;
+        return;
 
     struct Pipeline *pipeline = static_cast<struct Pipeline*>(P->parent->opaque);
     struct PushPop *pushpop = static_cast<struct PushPop*>(P->opaque);
@@ -664,8 +664,6 @@ static PJ_COORD pop(PJ_COORD point, PJ *P) {
             point.v[3] = pipeline->stack[3].top();
             pipeline->stack[3].pop();
     }
-
-    return point;
 }
 
 
diff --git a/src/proj_internal.h b/src/proj_internal.h
index c1eb9c0425..bf5bb797da 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -194,8 +194,8 @@ struct projCppContext;
 /* not sure why we need to export it, but mingw needs it */
 void PROJ_DLL proj_context_delete_cpp_context(struct projCppContext* cppContext);
 
-PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P);
-PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P);
+bool pj_fwd4d (PJ_COORD& coo, PJ *P);
+bool pj_inv4d (PJ_COORD& coo, PJ *P);
 
 PJ_COORD PROJ_DLL pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 PJ_COORD PROJ_DLL pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
@@ -277,13 +277,13 @@ struct PJ_AREA {
 
 PJ_OPERATOR:
 
-    A function taking a PJ_COORD and a pointer-to-PJ as args, applying the
-    PJ to the PJ_COORD, and returning the resulting PJ_COORD.
+    A function taking a reference to a PJ_COORD and a pointer-to-PJ as args, applying the
+    PJ to the PJ_COORD, and modifying in-place the passed PJ_COORD.
 
 *****************************************************************************/
 typedef    PJ       *(* PJ_CONSTRUCTOR) (PJ *);
 typedef    PJ       *(* PJ_DESTRUCTOR)  (PJ *, int);
-typedef    PJ_COORD  (* PJ_OPERATOR)    (PJ_COORD, PJ *);
+typedef    void      (* PJ_OPERATOR)    (PJ_COORD&, PJ *);
 /****************************************************************************/
 
 
@@ -577,6 +577,7 @@ struct PJconsts {
     **************************************************************************************/
 
     NS_PROJ::common::IdentifiedObjectPtr iso_obj{};
+    bool                                 iso_obj_is_coordinate_operation = false;
 
     // cached results
     mutable std::string lastWKT{};
diff --git a/src/projections/latlong.cpp b/src/projections/latlong.cpp
index 8be1521901..48d2291c13 100644
--- a/src/projections/latlong.cpp
+++ b/src/projections/latlong.cpp
@@ -74,15 +74,11 @@ static PJ_LPZ latlong_inverse_3d (PJ_XYZ xyz, PJ *P) {
     return lpz;
 }
 
-static PJ_COORD latlong_forward_4d (PJ_COORD obs, PJ *P) {
-    (void) P;
-    return obs;
+static void latlong_forward_4d (PJ_COORD&, PJ *) {
 }
 
 
-static PJ_COORD latlong_inverse_4d (PJ_COORD obs, PJ *P) {
-    (void) P;
-    return obs;
+static void latlong_inverse_4d (PJ_COORD&, PJ *) {
 }
 
 
diff --git a/src/transformations/affine.cpp b/src/transformations/affine.cpp
index f99c282fe1..297c9b5c44 100644
--- a/src/transformations/affine.cpp
+++ b/src/transformations/affine.cpp
@@ -58,55 +58,58 @@ struct pj_opaque_affine {
 } // anonymous namespace
 
 
-static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
-    PJ_COORD newObs;
+static void forward_4d(PJ_COORD& coo, PJ *P) {
     const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque;
     const struct pj_affine_coeffs *C = &(Q->forward);
-    newObs.xyzt.x = Q->xoff + C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z;
-    newObs.xyzt.y = Q->yoff + C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z;
-    newObs.xyzt.z = Q->zoff + C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z;
-    newObs.xyzt.t = Q->toff + C->tscale * obs.xyzt.t;
-    return newObs;
+    const double x = coo.xyz.x;
+    const double y = coo.xyz.y;
+    const double z = coo.xyz.z;
+    coo.xyzt.x = Q->xoff + C->s11 * x + C->s12 * y + C->s13 * z;
+    coo.xyzt.y = Q->yoff + C->s21 * x + C->s22 * y + C->s23 * z;
+    coo.xyzt.z = Q->zoff + C->s31 * x + C->s32 * y + C->s33 * z;
+    coo.xyzt.t = Q->toff + C->tscale * coo.xyzt.t;
 }
 
 static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
-    return forward_4d(point, P).xyz;
+    forward_4d(point, P);
+    return point.xyz;
 }
 
 
 static PJ_XY forward_2d(PJ_LP lp, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.lp = lp;
-    return forward_4d(point, P).xy;
+    forward_4d(point, P);
+    return point.xy;
 }
 
 
-static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
-    PJ_COORD newObs;
+static void reverse_4d(PJ_COORD& coo, PJ *P) {
     const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque;
     const struct pj_affine_coeffs *C = &(Q->reverse);
-    obs.xyzt.x -= Q->xoff;
-    obs.xyzt.y -= Q->yoff;
-    obs.xyzt.z -= Q->zoff;
-    newObs.xyzt.x = C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z;
-    newObs.xyzt.y = C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z;
-    newObs.xyzt.z = C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z;
-    newObs.xyzt.t = C->tscale * (obs.xyzt.t - Q->toff);
-    return newObs;
+    double x = coo.xyzt.x - Q->xoff;
+    double y = coo.xyzt.y - Q->yoff;
+    double z = coo.xyzt.z - Q->zoff;
+    coo.xyzt.x = C->s11 * x + C->s12 * y + C->s13 * z;
+    coo.xyzt.y = C->s21 * x + C->s22 * y + C->s23 * z;
+    coo.xyzt.z = C->s31 * x + C->s32 * y + C->s33 * z;
+    coo.xyzt.t = C->tscale * (coo.xyzt.t - Q->toff);
 }
 
 static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
-    return reverse_4d(point, P).lpz;
+    reverse_4d(point, P);
+    return point.lpz;
 }
 
 static PJ_LP reverse_2d(PJ_XY xy, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
     point.xy = xy;
-    return reverse_4d(point, P).lp;
+    reverse_4d(point, P);
+    return point.lp;
 }
 
 static struct pj_opaque_affine * initQ() {
diff --git a/src/transformations/defmodel.cpp b/src/transformations/defmodel.cpp
index 9e4590e1d5..03c8fd9cc1 100644
--- a/src/transformations/defmodel.cpp
+++ b/src/transformations/defmodel.cpp
@@ -343,34 +343,24 @@ static PJ *destructor(PJ *P, int errlev) {
     return pj_default_destructor(P, errlev);
 }
 
-static PJ_COORD forward_4d(PJ_COORD in, PJ *P) {
+static void forward_4d(PJ_COORD &coo, PJ *P) {
     auto *Q = (struct defmodelData *)P->opaque;
 
-    PJ_COORD out;
-    out.xyzt.t = in.xyzt.t;
-
-    if (!Q->evaluator->forward(Q->evaluatorIface, in.xyzt.x, in.xyzt.y,
-                               in.xyzt.z, in.xyzt.t, out.xyzt.x, out.xyzt.y,
-                               out.xyzt.z)) {
-        return proj_coord_error();
+    if (!Q->evaluator->forward(Q->evaluatorIface, coo.xyzt.x, coo.xyzt.y,
+                               coo.xyzt.z, coo.xyzt.t, coo.xyzt.x, coo.xyzt.y,
+                               coo.xyzt.z)) {
+        coo = proj_coord_error();
     }
-
-    return out;
 }
 
-static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) {
+static void reverse_4d(PJ_COORD &coo, PJ *P) {
     auto *Q = (struct defmodelData *)P->opaque;
 
-    PJ_COORD out;
-    out.xyzt.t = in.xyzt.t;
-
-    if (!Q->evaluator->inverse(Q->evaluatorIface, in.xyzt.x, in.xyzt.y,
-                               in.xyzt.z, in.xyzt.t, out.xyzt.x, out.xyzt.y,
-                               out.xyzt.z)) {
-        return proj_coord_error();
+    if (!Q->evaluator->inverse(Q->evaluatorIface, coo.xyzt.x, coo.xyzt.y,
+                               coo.xyzt.z, coo.xyzt.t, coo.xyzt.x, coo.xyzt.y,
+                               coo.xyzt.z)) {
+        coo = proj_coord_error();
     }
-
-    return out;
 }
 
 // Function called by proj_assign_context() when a new context is assigned to
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index 20bac9202b..7405ed6fa0 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -278,26 +278,22 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 }
 
 
-static PJ_COORD forward_4d(PJ_COORD in, PJ *P) {
+static void forward_4d(PJ_COORD& coo, PJ *P) {
     struct deformationData *Q = (struct deformationData *) P->opaque;
     double dt;
     PJ_XYZ shift;
-    PJ_COORD out = in;
 
     if (Q->dt != HUGE_VAL) {
         dt = Q->dt;
     } else {
-        dt = in.xyzt.t - Q->t_epoch ;
+        dt = coo.xyzt.t - Q->t_epoch ;
     }
 
-    shift = get_grid_shift(P, in.xyz);
-
-    out.xyzt.x += dt*shift.x;
-    out.xyzt.y += dt*shift.y;
-    out.xyzt.z += dt*shift.z;
-
+    shift = get_grid_shift(P, coo.xyz);
 
-    return out;
+    coo.xyzt.x += dt*shift.x;
+    coo.xyzt.y += dt*shift.y;
+    coo.xyzt.z += dt*shift.z;
 }
 
 
@@ -317,20 +313,17 @@ static PJ_LPZ reverse_3d(PJ_XYZ in, PJ *P) {
     return out.lpz;
 }
 
-static PJ_COORD reverse_4d(PJ_COORD in, PJ *P) {
+static void reverse_4d(PJ_COORD& coo, PJ *P) {
     struct deformationData *Q = (struct deformationData *) P->opaque;
-    PJ_COORD out = in;
     double dt;
 
-
     if (Q->dt != HUGE_VAL) {
             dt = Q->dt;
         } else {
-            dt = in.xyzt.t - Q->t_epoch;
+            dt = coo.xyzt.t - Q->t_epoch;
     }
 
-    out.xyz = reverse_shift(P, in.xyz, dt);
-    return out;
+    coo.xyz = reverse_shift(P, coo.xyz, dt);
 }
 
 static PJ *destructor(PJ *P, int errlev) {
diff --git a/src/transformations/helmert.cpp b/src/transformations/helmert.cpp
index 01ae02fa2f..027cd4b097 100644
--- a/src/transformations/helmert.cpp
+++ b/src/transformations/helmert.cpp
@@ -438,7 +438,7 @@ static PJ_LPZ helmert_reverse_3d (PJ_XYZ xyz, PJ *P) {
 }
 
 
-static PJ_COORD helmert_forward_4d (PJ_COORD point, PJ *P) {
+static void helmert_forward_4d (PJ_COORD &point, PJ *P) {
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
 
     /* We only need to rebuild the rotation matrix if the
@@ -450,14 +450,11 @@ static PJ_COORD helmert_forward_4d (PJ_COORD point, PJ *P) {
         build_rot_matrix(P);
     }
 
-    const auto xyz = helmert_forward_3d (point.lpz, P);
-    point.xyz = xyz;
-
-    return point;
+    point.xyz = helmert_forward_3d (point.lpz, P);
 }
 
 
-static PJ_COORD helmert_reverse_4d (PJ_COORD point, PJ *P) {
+static void helmert_reverse_4d (PJ_COORD& point, PJ *P) {
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
 
     /* We only need to rebuild the rotation matrix if the
@@ -469,10 +466,7 @@ static PJ_COORD helmert_reverse_4d (PJ_COORD point, PJ *P) {
         build_rot_matrix(P);
     }
 
-    const auto lpz = helmert_reverse_3d (point.xyz, P);
-    point.lpz = lpz;
-
-    return point;
+    point.lpz = helmert_reverse_3d (point.xyz, P);
 }
 
 /* Arcsecond to radians */
diff --git a/src/transformations/hgridshift.cpp b/src/transformations/hgridshift.cpp
index 322e29f92a..9ef8f0e1ea 100644
--- a/src/transformations/hgridshift.cpp
+++ b/src/transformations/hgridshift.cpp
@@ -70,39 +70,32 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
     return point.lpz;
 }
 
-static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+static void forward_4d(PJ_COORD& coo, PJ *P) {
     struct hgridshiftData *Q = (struct hgridshiftData *) P->opaque;
-    PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
     if (Q->t_final==0 || Q->t_epoch==0) {
-        point.xyz = forward_3d (obs.lpz, P);
-        return point;
+        coo.xyz = forward_3d (coo.lpz, P);
+        return;
     }
 
     /* Time restricted - only apply transform if within time bracket */
-    if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
-        point.xyz = forward_3d (obs.lpz, P);
-
-
-    return point;
+    if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
+        coo.xyz = forward_3d (coo.lpz, P);
 }
 
-static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+static void reverse_4d(PJ_COORD& coo, PJ *P) {
     struct hgridshiftData *Q = (struct hgridshiftData *) P->opaque;
-    PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
     if (Q->t_final==0 || Q->t_epoch==0) {
-        point.lpz = reverse_3d (obs.xyz, P);
-        return point;
+        coo.lpz = reverse_3d (coo.xyz, P);
+        return;
     }
 
     /* Time restricted - only apply transform if within time bracket */
-    if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
-        point.lpz = reverse_3d (obs.xyz, P);
-
-    return point;
+    if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
+        coo.lpz = reverse_3d (coo.xyz, P);
 }
 
 static PJ *destructor (PJ *P, int errlev) {
diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp
index 2c85151aae..4c94335390 100644
--- a/src/transformations/horner.cpp
+++ b/src/transformations/horner.cpp
@@ -383,22 +383,19 @@ static PJ_UV real_iterative_inverse_impl(PJ *P, const HORNER *transformation, PJ
     }
 }
 
-static PJ_COORD horner_forward_4d (PJ_COORD point, PJ *P) {
+static void horner_forward_4d (PJ_COORD& point, PJ *P) {
     const HORNER *transformation = reinterpret_cast<const HORNER*>(P->opaque);
     point.uv = real_default_impl(P, transformation, PJ_FWD, point.uv);
-    return point;
 }
 
-static PJ_COORD horner_inverse_4d (PJ_COORD point, PJ *P) {
+static void horner_inverse_4d (PJ_COORD& point, PJ *P) {
     const HORNER *transformation = reinterpret_cast<const HORNER*>(P->opaque);
     point.uv = real_default_impl(P, transformation, PJ_INV, point.uv);
-    return point;
 }
 
-static PJ_COORD horner_iterative_inverse_4d(PJ_COORD point, PJ *P) {
+static void horner_iterative_inverse_4d(PJ_COORD& point, PJ *P) {
     const HORNER *transformation = reinterpret_cast<const HORNER*>(P->opaque);
     point.uv = real_iterative_inverse_impl(P, transformation, point.uv);
-    return point;
 }
 
 static PJ_UV complex_default_impl(PJ *P, const HORNER *transformation, PJ_DIRECTION direction, PJ_UV position) {
@@ -481,22 +478,19 @@ static PJ_UV complex_iterative_inverse_impl(PJ *P, const HORNER *transformation,
     }
 }
 
-static PJ_COORD complex_horner_forward_4d (PJ_COORD point, PJ *P) {
+static void complex_horner_forward_4d (PJ_COORD& point, PJ *P) {
     const HORNER *transformation = reinterpret_cast<const HORNER*>(P->opaque);
     point.uv = complex_default_impl(P, transformation, PJ_FWD, point.uv);
-    return point;
 }
 
-static PJ_COORD complex_horner_inverse_4d (PJ_COORD point, PJ *P) {
+static void complex_horner_inverse_4d (PJ_COORD& point, PJ *P) {
     const HORNER *transformation = reinterpret_cast<const HORNER*>(P->opaque);
     point.uv = complex_default_impl(P, transformation, PJ_INV, point.uv);
-    return point;
 }
 
-static PJ_COORD complex_horner_iterative_inverse_4d (PJ_COORD point, PJ *P) {
+static void complex_horner_iterative_inverse_4d (PJ_COORD& point, PJ *P) {
     const HORNER *transformation = reinterpret_cast<const HORNER*>(P->opaque);
     point.uv = complex_iterative_inverse_impl(P, transformation, point.uv);
-    return point;
 }
 
 
diff --git a/src/transformations/molodensky.cpp b/src/transformations/molodensky.cpp
index 5d873cf000..d78a2db20c 100644
--- a/src/transformations/molodensky.cpp
+++ b/src/transformations/molodensky.cpp
@@ -260,10 +260,8 @@ static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
 }
 
 
-static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
-    const auto xyz = forward_3d(obs.lpz, P);
-    obs.xyz = xyz;
-    return obs;
+static void forward_4d(PJ_COORD& obs, PJ *P) {
+    obs.xyz = forward_3d(obs.lpz, P);
 }
 
 
@@ -293,10 +291,8 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
 }
 
 
-static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
-    const auto lpz = reverse_3d(obs.xyz, P);
-    obs.lpz = lpz;
-    return obs;
+static void reverse_4d(PJ_COORD& obs, PJ *P) {
+    obs.lpz = reverse_3d(obs.xyz, P);
 }
 
 
diff --git a/src/transformations/tinshift.cpp b/src/transformations/tinshift.cpp
index 4b45567d60..c780e5d4bd 100644
--- a/src/transformations/tinshift.cpp
+++ b/src/transformations/tinshift.cpp
@@ -60,26 +60,22 @@ static PJ *destructor(PJ *P, int errlev) {
     return pj_default_destructor(P, errlev);
 }
 
-static PJ_COORD tinshift_forward_4d(PJ_COORD in, PJ *P) {
+static void tinshift_forward_4d(PJ_COORD &coo, PJ *P) {
     auto *Q = (struct tinshiftData *)P->opaque;
 
-    PJ_COORD out = in;
-    if (!Q->evaluator->forward(in.xyz.x, in.xyz.y, in.xyz.z, out.xyz.x,
-                               out.xyz.y, out.xyz.z)) {
-        return proj_coord_error();
+    if (!Q->evaluator->forward(coo.xyz.x, coo.xyz.y, coo.xyz.z, coo.xyz.x,
+                               coo.xyz.y, coo.xyz.z)) {
+        coo = proj_coord_error();
     }
-    return out;
 }
 
-static PJ_COORD tinshift_reverse_4d(PJ_COORD in, PJ *P) {
+static void tinshift_reverse_4d(PJ_COORD &coo, PJ *P) {
     auto *Q = (struct tinshiftData *)P->opaque;
 
-    PJ_COORD out = in;
-    if (!Q->evaluator->inverse(in.xyz.x, in.xyz.y, in.xyz.z, out.xyz.x,
-                               out.xyz.y, out.xyz.z)) {
-        return proj_coord_error();
+    if (!Q->evaluator->inverse(coo.xyz.x, coo.xyz.y, coo.xyz.z, coo.xyz.x,
+                               coo.xyz.y, coo.xyz.z)) {
+        coo = proj_coord_error();
     }
-    return out;
 }
 
 PJ *TRANSFORMATION(tinshift, 1) {
diff --git a/src/transformations/vgridshift.cpp b/src/transformations/vgridshift.cpp
index 5ebf74cd03..6218553d93 100644
--- a/src/transformations/vgridshift.cpp
+++ b/src/transformations/vgridshift.cpp
@@ -101,39 +101,32 @@ static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
 }
 
 
-static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+static void forward_4d(PJ_COORD& coo, PJ *P) {
     struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
-    PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
     if (Q->t_final==0 || Q->t_epoch==0) {
-        point.xyz = forward_3d (obs.lpz, P);
-        return point;
+        coo.xyz = forward_3d (coo.lpz, P);
+        return;
     }
 
     /* Time restricted - only apply transform if within time bracket */
-    if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
-        point.xyz = forward_3d (obs.lpz, P);
-
-
-    return point;
+    if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
+        coo.xyz = forward_3d (coo.lpz, P);
 }
 
-static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
+static void reverse_4d(PJ_COORD& coo, PJ *P) {
     struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
-    PJ_COORD point = obs;
 
     /* If transformation is not time restricted, we always call it */
     if (Q->t_final==0 || Q->t_epoch==0) {
-        point.lpz = reverse_3d (obs.xyz, P);
-        return point;
+        coo.lpz = reverse_3d (coo.xyz, P);
+        return;
     }
 
     /* Time restricted - only apply transform if within time bracket */
-    if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
-        point.lpz = reverse_3d (obs.xyz, P);
-
-    return point;
+    if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
+        coo.lpz = reverse_3d (coo.xyz, P);
 }
 
 static PJ *destructor (PJ *P, int errlev) {