Merge pull request #2206 from rouault/deformation_model_for_merge

 Add a +proj=defmodel transformation for multi-component time-based deformation models

13 files changed, 2654 insertions, 181 deletions

diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 2b00011d..8f427412 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -1681,14 +1681,14 @@ PJ_GRID_INFO proj_grid_info(const char *gridname) {
         grinfo.n_lat = grid.height();
 
         /* cell size */
-        grinfo.cs_lon = extent.resLon;
-        grinfo.cs_lat = extent.resLat;
+        grinfo.cs_lon = extent.resX;
+        grinfo.cs_lat = extent.resY;
 
         /* bounds of grid */
-        grinfo.lowerleft.lam  = extent.westLon;
-        grinfo.lowerleft.phi  = extent.southLat;
-        grinfo.upperright.lam = extent.eastLon;
-        grinfo.upperright.phi = extent.northLat;
+        grinfo.lowerleft.lam  = extent.west;
+        grinfo.lowerleft.phi  = extent.south;
+        grinfo.upperright.lam = extent.east;
+        grinfo.upperright.phi = extent.north;
     };
 
     {
diff --git a/src/Makefile.am b/src/Makefile.am
index 3feaa810..9f24a60d 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -196,6 +196,10 @@ libproj_la_SOURCES = \
 	transformations/molodensky.cpp \
 	transformations/vgridshift.cpp \
 	transformations/xyzgridshift.cpp \
+	transformations/defmodel.cpp \
+	transformations/defmodel.hpp \
+	transformations/defmodel_exceptions.hpp \
+	transformations/defmodel_impl.hpp \
 	\
 	aasincos.cpp adjlon.cpp \
 	dmstor.cpp auth.cpp \
diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp
index 6ce0ab18..8940afde 100644
--- a/src/apps/gie.cpp
+++ b/src/apps/gie.cpp
@@ -723,18 +723,27 @@ Attempt to interpret args as a PJ_COORD.
         /* This could be avoided if proj_dmstor used the same proj_strtod()  */
         /* as gie, but that is not the case (yet). When we remove projects.h */
         /* from the public API we can change that.                           */
+
+        // Even Rouault: unsure about the above. Coordinates are not necessarily
+        // geographic coordinates, and the roundtrip through radians for
+        // big projected coordinates cause inaccuracies, that can cause
+        // test failures when testing points at edge of grids.
+        // For example 1501000.0 becomes 1501000.000000000233
         double d = proj_strtod(prev,  (char **) &endp);
-        double dms = PJ_TODEG(proj_dmstor (prev, (char **) &dmsendp));
-       /* TODO: When projects.h is removed, call proj_dmstor() in all cases */
-        if (d != dms && fabs(d) < fabs(dms) && fabs(dms) < fabs(d) + 1) {
-            d = dms;
-            endp = dmsendp;
+        if( *endp != '\0' && !isspace(*endp) )
+        {
+            double dms = PJ_TODEG(proj_dmstor (prev, (char **) &dmsendp));
+            /* TODO: When projects.h is removed, call proj_dmstor() in all cases */
+            if (d != dms && fabs(d) < fabs(dms) && fabs(dms) < fabs(d) + 1) {
+                d = dms;
+                endp = dmsendp;
+            }
+            /* A number like -81d00'00.000 will be parsed correctly by both */
+            /* proj_strtod and proj_dmstor but only the latter will return  */
+            /* the correct end-pointer.                                     */
+            if (d == dms && endp != dmsendp)
+                endp = dmsendp;
         }
-        /* A number like -81d00'00.000 will be parsed correctly by both */
-        /* proj_strtod and proj_dmstor but only the latter will return  */
-        /* the correct end-pointer.                                     */
-        if (d == dms && endp != dmsendp)
-            endp = dmsendp;
 
         /* Break out if there were no more numerals */
         if (prev==endp)
diff --git a/src/grids.cpp b/src/grids.cpp
index 067d3d45..8065813a 100644
--- a/src/grids.cpp
+++ b/src/grids.cpp
@@ -77,21 +77,21 @@ static void swap_words(void *dataIn, size_t word_size, size_t word_count)
 // ---------------------------------------------------------------------------
 
 bool ExtentAndRes::fullWorldLongitude() const {
-    return eastLon - westLon + resLon >= 2 * M_PI - 1e-10;
+    return isGeographic && east - west + resX >= 2 * M_PI - 1e-10;
 }
 
 // ---------------------------------------------------------------------------
 
 bool ExtentAndRes::contains(const ExtentAndRes &other) const {
-    return other.westLon >= westLon && other.eastLon <= eastLon &&
-           other.southLat >= southLat && other.northLat <= northLat;
+    return other.west >= west && other.east <= east && other.south >= south &&
+           other.north <= north;
 }
 
 // ---------------------------------------------------------------------------
 
 bool ExtentAndRes::intersects(const ExtentAndRes &other) const {
-    return other.westLon < eastLon && westLon <= other.westLon &&
-           other.southLat < northLat && southLat <= other.northLat;
+    return other.west < east && west <= other.west && other.south < north &&
+           south <= other.north;
 }
 
 // ---------------------------------------------------------------------------
@@ -119,12 +119,13 @@ VerticalShiftGrid::~VerticalShiftGrid() = default;
 
 static ExtentAndRes globalExtent() {
     ExtentAndRes extent;
-    extent.westLon = -M_PI;
-    extent.southLat = -M_PI / 2;
-    extent.eastLon = M_PI;
-    extent.northLat = M_PI / 2;
-    extent.resLon = M_PI;
-    extent.resLat = M_PI / 2;
+    extent.isGeographic = true;
+    extent.west = -M_PI;
+    extent.south = -M_PI / 2;
+    extent.east = M_PI;
+    extent.north = M_PI / 2;
+    extent.resX = M_PI;
+    extent.resY = M_PI / 2;
     return extent;
 }
 
@@ -238,12 +239,13 @@ GTXVerticalShiftGrid *GTXVerticalShiftGrid::open(PJ_CONTEXT *ctx,
     }
 
     ExtentAndRes extent;
-    extent.westLon = xorigin * DEG_TO_RAD;
-    extent.southLat = yorigin * DEG_TO_RAD;
-    extent.resLon = xstep * DEG_TO_RAD;
-    extent.resLat = ystep * DEG_TO_RAD;
-    extent.eastLon = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD;
-    extent.northLat = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD;
+    extent.isGeographic = true;
+    extent.west = xorigin * DEG_TO_RAD;
+    extent.south = yorigin * DEG_TO_RAD;
+    extent.resX = xstep * DEG_TO_RAD;
+    extent.resY = ystep * DEG_TO_RAD;
+    extent.east = (xorigin + xstep * (columns - 1)) * DEG_TO_RAD;
+    extent.north = (yorigin + ystep * (rows - 1)) * DEG_TO_RAD;
 
     return new GTXVerticalShiftGrid(ctx, std::move(fp), name, columns, rows,
                                     extent);
@@ -931,6 +933,10 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
     unsigned short count = 0;
     unsigned short *geokeys = nullptr;
     bool pixelIsArea = false;
+
+    ExtentAndRes extent;
+    extent.isGeographic = true;
+
     if (!TIFFGetField(m_hTIFF, TIFFTAG_GEOKEYDIRECTORY, &count, &geokeys)) {
         pj_log(m_ctx, PJ_LOG_DEBUG_MINOR, "No GeoKeys tag");
     } else {
@@ -953,6 +959,7 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
 
         for (unsigned int i = 4; i + 3 < count; i += 4) {
             constexpr unsigned short GTModelTypeGeoKey = 1024;
+            constexpr unsigned short ModelTypeProjected = 1;
             constexpr unsigned short ModelTypeGeographic = 2;
 
             constexpr unsigned short GTRasterTypeGeoKey = 1025;
@@ -960,10 +967,13 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
             // constexpr unsigned short RasterPixelIsPoint = 2;
 
             if (geokeys[i] == GTModelTypeGeoKey) {
-                if (geokeys[i + 3] != ModelTypeGeographic) {
-                    pj_log(m_ctx, PJ_LOG_ERROR, "Only GTModelTypeGeoKey = "
-                                                "ModelTypeGeographic is "
-                                                "supported");
+                if (geokeys[i + 3] == ModelTypeProjected) {
+                    extent.isGeographic = false;
+                } else if (geokeys[i + 3] != ModelTypeGeographic) {
+                    pj_log(m_ctx, PJ_LOG_ERROR,
+                           "Only GTModelTypeGeoKey = "
+                           "ModelTypeGeographic or ModelTypeProjected are "
+                           "supported");
                     return nullptr;
                 }
             } else if (geokeys[i] == GTRasterTypeGeoKey) {
@@ -976,8 +986,8 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
 
     double hRes = 0;
     double vRes = 0;
-    double westLon = 0;
-    double northLat = 0;
+    double west = 0;
+    double north = 0;
 
     double *matrix = nullptr;
     if (TIFFGetField(m_hTIFF, TIFFTAG_GEOTRANSMATRIX, &count, &matrix) &&
@@ -991,9 +1001,9 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
             return nullptr;
         }
 
-        westLon = matrix[3];
+        west = matrix[3];
         hRes = matrix[0];
-        northLat = matrix[7];
+        north = matrix[7];
         vRes = -matrix[5]; // negation to simulate GeoPixelScale convention
     } else {
         double *geopixelscale = nullptr;
@@ -1022,34 +1032,33 @@ std::unique_ptr<GTiffGrid> GTiffDataset::nextGrid() {
             return nullptr;
         }
 
-        westLon = geotiepoints[3] - geotiepoints[0] * hRes;
-        northLat = geotiepoints[4] + geotiepoints[1] * vRes;
+        west = geotiepoints[3] - geotiepoints[0] * hRes;
+        north = geotiepoints[4] + geotiepoints[1] * vRes;
     }
 
     if (pixelIsArea) {
-        westLon += 0.5 * hRes;
-        northLat -= 0.5 * vRes;
+        west += 0.5 * hRes;
+        north -= 0.5 * vRes;
     }
 
-    ExtentAndRes extent;
-    extent.westLon = westLon * DEG_TO_RAD;
-    extent.northLat = northLat * DEG_TO_RAD;
-    extent.resLon = hRes * DEG_TO_RAD;
-    extent.resLat = fabs(vRes) * DEG_TO_RAD;
-    extent.eastLon = (westLon + hRes * (width - 1)) * DEG_TO_RAD;
-    extent.southLat = (northLat - vRes * (height - 1)) * DEG_TO_RAD;
+    const double mulFactor = extent.isGeographic ? DEG_TO_RAD : 1;
+    extent.west = west * mulFactor;
+    extent.north = north * mulFactor;
+    extent.resX = hRes * mulFactor;
+    extent.resY = fabs(vRes) * mulFactor;
+    extent.east = (west + hRes * (width - 1)) * mulFactor;
+    extent.south = (north - vRes * (height - 1)) * mulFactor;
 
     if (vRes < 0) {
-        std::swap(extent.northLat, extent.southLat);
+        std::swap(extent.north, extent.south);
     }
 
-    if (!(fabs(extent.westLon) <= 4 * M_PI &&
-          fabs(extent.eastLon) <= 4 * M_PI &&
-          fabs(extent.northLat) <= M_PI + 1e-5 &&
-          fabs(extent.southLat) <= M_PI + 1e-5 &&
-          extent.westLon < extent.eastLon &&
-          extent.southLat < extent.northLat && extent.resLon > 1e-10 &&
-          extent.resLat > 1e-10)) {
+    if (!((!extent.isGeographic ||
+           (fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI &&
+            fabs(extent.north) <= M_PI + 1e-5 &&
+            fabs(extent.south) <= M_PI + 1e-5)) &&
+          extent.west < extent.east && extent.south < extent.north &&
+          extent.resX > 1e-10 && extent.resY > 1e-10)) {
         pj_log(m_ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
                m_filename.c_str());
         return nullptr;
@@ -1408,17 +1417,19 @@ bool VerticalShiftGridSet::reopen(PJ_CONTEXT *ctx) {
 
 // ---------------------------------------------------------------------------
 
-static bool isPointInExtent(double lon, double lat, const ExtentAndRes &extent,
+static bool isPointInExtent(double x, double y, const ExtentAndRes &extent,
                             double eps = 0) {
-    if (!(lat + eps >= extent.southLat && lat - eps <= extent.northLat))
+    if (!(y + eps >= extent.south && y - eps <= extent.north))
         return false;
     if (extent.fullWorldLongitude())
         return true;
-    if (lon + eps < extent.westLon)
-        lon += 2 * M_PI;
-    else if (lon - eps > extent.eastLon)
-        lon -= 2 * M_PI;
-    if (!(lon + eps >= extent.westLon && lon - eps <= extent.eastLon))
+    if (extent.isGeographic) {
+        if (x + eps < extent.west)
+            x += 2 * M_PI;
+        else if (x - eps > extent.east)
+            x -= 2 * M_PI;
+    }
+    if (!(x + eps >= extent.west && x - eps <= extent.east))
         return false;
     return true;
 }
@@ -1584,28 +1595,27 @@ NTv1Grid *NTv1Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
     }
 
     ExtentAndRes extent;
-    extent.westLon = -to_double(header + 72) * DEG_TO_RAD;
-    extent.southLat = to_double(header + 24) * DEG_TO_RAD;
-    extent.eastLon = -to_double(header + 56) * DEG_TO_RAD;
-    extent.northLat = to_double(header + 40) * DEG_TO_RAD;
-    extent.resLon = to_double(header + 104) * DEG_TO_RAD;
-    extent.resLat = to_double(header + 88) * DEG_TO_RAD;
-    if (!(fabs(extent.westLon) <= 4 * M_PI &&
-          fabs(extent.eastLon) <= 4 * M_PI &&
-          fabs(extent.northLat) <= M_PI + 1e-5 &&
-          fabs(extent.southLat) <= M_PI + 1e-5 &&
-          extent.westLon < extent.eastLon &&
-          extent.southLat < extent.northLat && extent.resLon > 1e-10 &&
-          extent.resLat > 1e-10)) {
+    extent.isGeographic = true;
+    extent.west = -to_double(header + 72) * DEG_TO_RAD;
+    extent.south = to_double(header + 24) * DEG_TO_RAD;
+    extent.east = -to_double(header + 56) * DEG_TO_RAD;
+    extent.north = to_double(header + 40) * DEG_TO_RAD;
+    extent.resX = to_double(header + 104) * DEG_TO_RAD;
+    extent.resY = to_double(header + 88) * DEG_TO_RAD;
+    if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI &&
+          fabs(extent.north) <= M_PI + 1e-5 &&
+          fabs(extent.south) <= M_PI + 1e-5 && extent.west < extent.east &&
+          extent.south < extent.north && extent.resX > 1e-10 &&
+          extent.resY > 1e-10)) {
         pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
                filename.c_str());
         pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
         return nullptr;
     }
     const int columns = static_cast<int>(
-        fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1);
+        fabs((extent.east - extent.west) / extent.resX + 0.5) + 1);
     const int rows = static_cast<int>(
-        fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) + 1);
+        fabs((extent.north - extent.south) / extent.resY + 0.5) + 1);
 
     return new NTv1Grid(ctx, std::move(fp), filename, columns, rows, extent);
 }
@@ -1695,17 +1705,17 @@ CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
     }
 
     ExtentAndRes extent;
-    static_assert(sizeof(extent.westLon) == 8, "wrong sizeof");
-    static_assert(sizeof(extent.southLat) == 8, "wrong sizeof");
-    static_assert(sizeof(extent.resLon) == 8, "wrong sizeof");
-    static_assert(sizeof(extent.resLat) == 8, "wrong sizeof");
-    memcpy(&extent.westLon, header + 96, 8);
-    memcpy(&extent.southLat, header + 104, 8);
-    memcpy(&extent.resLon, header + 112, 8);
-    memcpy(&extent.resLat, header + 120, 8);
-    if (!(fabs(extent.westLon) <= 4 * M_PI &&
-          fabs(extent.southLat) <= M_PI + 1e-5 && extent.resLon > 1e-10 &&
-          extent.resLat > 1e-10)) {
+    extent.isGeographic = true;
+    static_assert(sizeof(extent.west) == 8, "wrong sizeof");
+    static_assert(sizeof(extent.south) == 8, "wrong sizeof");
+    static_assert(sizeof(extent.resX) == 8, "wrong sizeof");
+    static_assert(sizeof(extent.resY) == 8, "wrong sizeof");
+    memcpy(&extent.west, header + 96, 8);
+    memcpy(&extent.south, header + 104, 8);
+    memcpy(&extent.resX, header + 112, 8);
+    memcpy(&extent.resY, header + 120, 8);
+    if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.south) <= M_PI + 1e-5 &&
+          extent.resX > 1e-10 && extent.resY > 1e-10)) {
         pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
                filename.c_str());
         pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
@@ -1719,8 +1729,8 @@ CTable2Grid *CTable2Grid::open(PJ_CONTEXT *ctx, std::unique_ptr<File> fp,
         pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
         return nullptr;
     }
-    extent.eastLon = extent.westLon + (width - 1) * extent.resLon;
-    extent.northLat = extent.southLat + (height - 1) * extent.resLon;
+    extent.east = extent.west + (width - 1) * extent.resX;
+    extent.north = extent.south + (height - 1) * extent.resX;
 
     return new CTable2Grid(ctx, std::move(fp), filename, width, height, extent);
 }
@@ -1902,8 +1912,8 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx,
         constexpr int OFFSET_GS_COUNT = 8 + 16 * 10;
         constexpr int OFFSET_SOUTH_LAT = 8 + 16 * 4;
         if (must_swap) {
-            // 6 double values: southLat, northLat, eastLon, westLon, resLat,
-            // resLon
+            // 6 double values: south, north, east, west, resY,
+            // resX
             for (int i = 0; i < 6; i++) {
                 swap_words(header + OFFSET_SOUTH_LAT + 16 * i, sizeof(double),
                            1);
@@ -1915,42 +1925,40 @@ std::unique_ptr<NTv2GridSet> NTv2GridSet::open(PJ_CONTEXT *ctx,
         gridName.append(header + 8, 8);
 
         ExtentAndRes extent;
-        extent.southLat = to_double(header + OFFSET_SOUTH_LAT) * DEG_TO_RAD /
-                          3600.0; /* S_LAT */
-        extent.northLat = to_double(header + OFFSET_SOUTH_LAT + 16) *
-                          DEG_TO_RAD / 3600.0; /* N_LAT */
-        extent.eastLon = -to_double(header + OFFSET_SOUTH_LAT + 16 * 2) *
-                         DEG_TO_RAD / 3600.0; /* E_LONG */
-        extent.westLon = -to_double(header + OFFSET_SOUTH_LAT + 16 * 3) *
-                         DEG_TO_RAD / 3600.0; /* W_LONG */
-        extent.resLat =
+        extent.isGeographic = true;
+        extent.south = to_double(header + OFFSET_SOUTH_LAT) * DEG_TO_RAD /
+                       3600.0; /* S_LAT */
+        extent.north = to_double(header + OFFSET_SOUTH_LAT + 16) * DEG_TO_RAD /
+                       3600.0; /* N_LAT */
+        extent.east = -to_double(header + OFFSET_SOUTH_LAT + 16 * 2) *
+                      DEG_TO_RAD / 3600.0; /* E_LONG */
+        extent.west = -to_double(header + OFFSET_SOUTH_LAT + 16 * 3) *
+                      DEG_TO_RAD / 3600.0; /* W_LONG */
+        extent.resY =
             to_double(header + OFFSET_SOUTH_LAT + 16 * 4) * DEG_TO_RAD / 3600.0;
-        extent.resLon =
+        extent.resX =
             to_double(header + OFFSET_SOUTH_LAT + 16 * 5) * DEG_TO_RAD / 3600.0;
 
-        if (!(fabs(extent.westLon) <= 4 * M_PI &&
-              fabs(extent.eastLon) <= 4 * M_PI &&
-              fabs(extent.northLat) <= M_PI + 1e-5 &&
-              fabs(extent.southLat) <= M_PI + 1e-5 &&
-              extent.westLon < extent.eastLon &&
-              extent.southLat < extent.northLat && extent.resLon > 1e-10 &&
-              extent.resLat > 1e-10)) {
+        if (!(fabs(extent.west) <= 4 * M_PI && fabs(extent.east) <= 4 * M_PI &&
+              fabs(extent.north) <= M_PI + 1e-5 &&
+              fabs(extent.south) <= M_PI + 1e-5 && extent.west < extent.east &&
+              extent.south < extent.north && extent.resX > 1e-10 &&
+              extent.resY > 1e-10)) {
             pj_log(ctx, PJ_LOG_ERROR, "Inconsistent georeferencing for %s",
                    filename.c_str());
             pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
             return nullptr;
         }
         const int columns = static_cast<int>(
-            fabs((extent.eastLon - extent.westLon) / extent.resLon + 0.5) + 1);
+            fabs((extent.east - extent.west) / extent.resX + 0.5) + 1);
         const int rows = static_cast<int>(
-            fabs((extent.northLat - extent.southLat) / extent.resLat + 0.5) +
-            1);
+            fabs((extent.north - extent.south) / extent.resY + 0.5) + 1);
 
         pj_log(ctx, PJ_LOG_DEBUG_MINOR,
                "NTv2 %s %dx%d: LL=(%.9g,%.9g) UR=(%.9g,%.9g)", gridName.c_str(),
-               columns, rows, extent.westLon * RAD_TO_DEG,
-               extent.southLat * RAD_TO_DEG, extent.eastLon * RAD_TO_DEG,
-               extent.northLat * RAD_TO_DEG);
+               columns, rows, extent.west * RAD_TO_DEG,
+               extent.south * RAD_TO_DEG, extent.east * RAD_TO_DEG,
+               extent.north * RAD_TO_DEG);
 
         unsigned int gs_count;
         memcpy(&gs_count, header + OFFSET_GS_COUNT, 4);
@@ -2393,8 +2401,8 @@ const HorizontalShiftGrid *HorizontalShiftGrid::gridAt(double lon,
                                                        double lat) const {
     for (const auto &child : m_children) {
         const auto &extentChild = child->extentAndRes();
-        const double epsilon = (extentChild.resLon + extentChild.resLat) *
-                               REL_TOLERANCE_HGRIDSHIFT;
+        const double epsilon =
+            (extentChild.resX + extentChild.resY) * REL_TOLERANCE_HGRIDSHIFT;
         if (isPointInExtent(lon, lat, extentChild, epsilon)) {
             return child->gridAt(lon, lat);
         }
@@ -2411,7 +2419,7 @@ const HorizontalShiftGrid *HorizontalShiftGridSet::gridAt(double lon,
         }
         const auto &extent = grid->extentAndRes();
         const double epsilon =
-            (extent.resLon + extent.resLat) * REL_TOLERANCE_HGRIDSHIFT;
+            (extent.resX + extent.resY) * REL_TOLERANCE_HGRIDSHIFT;
         if (isPointInExtent(lon, lat, extent, epsilon)) {
             return grid->gridAt(lon, lat);
         }
@@ -2723,11 +2731,11 @@ bool GenericShiftGridSet::reopen(PJ_CONTEXT *ctx) {
 
 // ---------------------------------------------------------------------------
 
-const GenericShiftGrid *GenericShiftGrid::gridAt(double lon, double lat) const {
+const GenericShiftGrid *GenericShiftGrid::gridAt(double x, double y) const {
     for (const auto &child : m_children) {
         const auto &extentChild = child->extentAndRes();
-        if (isPointInExtent(lon, lat, extentChild)) {
-            return child->gridAt(lon, lat);
+        if (isPointInExtent(x, y, extentChild)) {
+            return child->gridAt(x, y);
         }
     }
     return this;
@@ -2735,15 +2743,14 @@ const GenericShiftGrid *GenericShiftGrid::gridAt(double lon, double lat) const {
 
 // ---------------------------------------------------------------------------
 
-const GenericShiftGrid *GenericShiftGridSet::gridAt(double lon,
-                                                    double lat) const {
+const GenericShiftGrid *GenericShiftGridSet::gridAt(double x, double y) const {
     for (const auto &grid : m_grids) {
         if (dynamic_cast<NullGenericShiftGrid *>(grid.get())) {
             return grid.get();
         }
         const auto &extent = grid->extentAndRes();
-        if (isPointInExtent(lon, lat, extent)) {
-            return grid->gridAt(lon, lat);
+        if (isPointInExtent(x, y, extent)) {
+            return grid->gridAt(x, y);
         }
     }
     return nullptr;
@@ -2872,9 +2879,9 @@ static PJ_LP pj_hgrid_interpolate(PJ_LP t, const HorizontalShiftGrid *grid,
     int in;
 
     const auto &extent = grid->extentAndRes();
-    t.lam /= extent.resLon;
+    t.lam /= extent.resX;
     indx.lam = std::isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam));
-    t.phi /= extent.resLat;
+    t.phi /= extent.resY;
     indx.phi = std::isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi));
 
     frct.lam = t.lam - indx.lam;
@@ -2959,13 +2966,13 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in,
     tb = in;
     const auto *extent = &(grid->extentAndRes());
     const double epsilon =
-        (extent->resLon + extent->resLat) * REL_TOLERANCE_HGRIDSHIFT;
-    tb.lam -= extent->westLon;
+        (extent->resX + extent->resY) * REL_TOLERANCE_HGRIDSHIFT;
+    tb.lam -= extent->west;
     if (tb.lam + epsilon < 0)
         tb.lam += 2 * M_PI;
-    else if (tb.lam - epsilon > extent->eastLon - extent->westLon)
+    else if (tb.lam - epsilon > extent->east - extent->west)
         tb.lam -= 2 * M_PI;
-    tb.phi -= extent->southLat;
+    tb.phi -= extent->south;
 
     t = pj_hgrid_interpolate(tb, grid, true);
     if (grid->hasChanged()) {
@@ -2995,8 +3002,8 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in,
         /* to fetch a new grid into which iterate... */
         if (del.lam == HUGE_VAL) {
             PJ_LP lp;
-            lp.lam = t.lam + extent->westLon;
-            lp.phi = t.phi + extent->southLat;
+            lp.lam = t.lam + extent->west;
+            lp.phi = t.phi + extent->south;
             auto newGrid = findGrid(grids, lp, gridset);
             if (newGrid == nullptr || newGrid == grid || newGrid->isNullGrid())
                 break;
@@ -3004,15 +3011,15 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in,
                    grid->name().c_str(), newGrid->name().c_str());
             grid = newGrid;
             extent = &(grid->extentAndRes());
-            t.lam = lp.lam - extent->westLon;
-            t.phi = lp.phi - extent->southLat;
+            t.lam = lp.lam - extent->west;
+            t.phi = lp.phi - extent->south;
             tb = in;
-            tb.lam -= extent->westLon;
+            tb.lam -= extent->west;
             if (tb.lam + epsilon < 0)
                 tb.lam += 2 * M_PI;
-            else if (tb.lam - epsilon > extent->eastLon - extent->westLon)
+            else if (tb.lam - epsilon > extent->east - extent->west)
                 tb.lam -= 2 * M_PI;
-            tb.phi -= extent->southLat;
+            tb.phi -= extent->south;
             dif.lam = std::numeric_limits<double>::max();
             dif.phi = std::numeric_limits<double>::max();
             continue;
@@ -3041,8 +3048,8 @@ static PJ_LP pj_hgrid_apply_internal(projCtx ctx, PJ_LP in,
         fprintf(stderr, "Inverse grid shift iteration failed, presumably at "
                         "grid edge.\nUsing first approximation.\n");
 
-    in.lam = adjlon(t.lam + extent->westLon);
-    in.phi = t.phi + extent->southLat;
+    in.lam = adjlon(t.lam + extent->west);
+    in.phi = t.phi + extent->south;
     return in;
 }
 
@@ -3097,14 +3104,21 @@ PJ_LP pj_hgrid_value(PJ *P, const ListOfHGrids &grids, PJ_LP lp) {
 
     /* normalize input to ll origin */
     const auto &extent = grid->extentAndRes();
+    if (!extent.isGeographic) {
+        pj_log(P->ctx, PJ_LOG_ERROR,
+               "Can only handle grids referenced in a geographic CRS");
+        pj_ctx_set_errno(P->ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return out;
+    }
+
     const double epsilon =
-        (extent.resLon + extent.resLat) * REL_TOLERANCE_HGRIDSHIFT;
-    lp.lam -= extent.westLon;
+        (extent.resX + extent.resY) * REL_TOLERANCE_HGRIDSHIFT;
+    lp.lam -= extent.west;
     if (lp.lam + epsilon < 0)
         lp.lam += 2 * M_PI;
-    else if (lp.lam - epsilon > extent.eastLon - extent.westLon)
+    else if (lp.lam - epsilon > extent.east - extent.west)
         lp.lam -= 2 * M_PI;
-    lp.phi -= extent.southLat;
+    lp.phi -= extent.south;
 
     out = pj_hgrid_interpolate(lp, grid, false);
     if (grid->hasChanged()) {
@@ -3151,10 +3165,16 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
     }
 
     const auto &extent = grid->extentAndRes();
+    if (!extent.isGeographic) {
+        pj_log(ctx, PJ_LOG_ERROR,
+               "Can only handle grids referenced in a geographic CRS");
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return HUGE_VAL;
+    }
 
     /* Interpolation of a location within the grid */
-    double grid_x = (input.lam - extent.westLon) / extent.resLon;
-    if (input.lam < extent.westLon) {
+    double grid_x = (input.lam - extent.west) / extent.resX;
+    if (input.lam < extent.west) {
         if (extent.fullWorldLongitude()) {
             // The first fmod goes to ]-lim, lim[ range
             // So we add lim again to be in ]0, 2*lim[ and fmod again
@@ -3162,9 +3182,9 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
                               grid->width(),
                           grid->width());
         } else {
-            grid_x = (input.lam + 2 * M_PI - extent.westLon) / extent.resLon;
+            grid_x = (input.lam + 2 * M_PI - extent.west) / extent.resX;
         }
-    } else if (input.lam > extent.eastLon) {
+    } else if (input.lam > extent.east) {
         if (extent.fullWorldLongitude()) {
             // The first fmod goes to ]-lim, lim[ range
             // So we add lim again to be in ]0, 2*lim[ and fmod again
@@ -3172,10 +3192,10 @@ static double read_vgrid_value(PJ_CONTEXT *ctx, const ListOfVGrids &grids,
                               grid->width(),
                           grid->width());
         } else {
-            grid_x = (input.lam - 2 * M_PI - extent.westLon) / extent.resLon;
+            grid_x = (input.lam - 2 * M_PI - extent.west) / extent.resX;
         }
     }
-    double grid_y = (input.phi - extent.southLat) / extent.resLat;
+    double grid_y = (input.phi - extent.south) / extent.resY;
     int grid_ix = static_cast<int>(lround(floor(grid_x)));
     if (!(grid_ix >= 0 && grid_ix < grid->width())) {
         // in the unlikely case we end up here...
@@ -3342,11 +3362,9 @@ const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids,
 
 // Used by +proj=deformation and +proj=xyzgridshift to do bilinear interpolation
 // on 3 sample values per node.
-bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid,
-                                             const PJ_LP &lp, int idx1,
-                                             int idx2, int idx3, double &v1,
-                                             double &v2, double &v3,
-                                             bool &must_retry) {
+bool pj_bilinear_interpolation_three_samples(
+    PJ_CONTEXT *ctx, const GenericShiftGrid *grid, const PJ_LP &lp, int idx1,
+    int idx2, int idx3, double &v1, double &v2, double &v3, bool &must_retry) {
     must_retry = false;
     if (grid->isNullGrid()) {
         v1 = 0.0;
@@ -3356,13 +3374,25 @@ bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid,
     }
 
     const auto &extent = grid->extentAndRes();
-    double grid_x = (lp.lam - extent.westLon) / extent.resLon;
-    if (lp.lam < extent.westLon) {
-        grid_x = (lp.lam + 2 * M_PI - extent.westLon) / extent.resLon;
-    } else if (lp.lam > extent.eastLon) {
-        grid_x = (lp.lam - 2 * M_PI - extent.westLon) / extent.resLon;
+    if (!extent.isGeographic) {
+        pj_log(ctx, PJ_LOG_ERROR,
+               "Can only handle grids referenced in a geographic CRS");
+        pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
+        return false;
+    }
+
+    // From a input location lp, determine the grid cell into which it falls,
+    // by identifying the lower-left x,y of it (ix, iy), and the upper-right
+    // (ix2, iy2)
+
+    double grid_x = (lp.lam - extent.west) / extent.resX;
+    // Special case for grids with world extent, and dealing with wrap-around
+    if (lp.lam < extent.west) {
+        grid_x = (lp.lam + 2 * M_PI - extent.west) / extent.resX;
+    } else if (lp.lam > extent.east) {
+        grid_x = (lp.lam - 2 * M_PI - extent.west) / extent.resX;
     }
-    double grid_y = (lp.phi - extent.southLat) / extent.resLat;
+    double grid_y = (lp.phi - extent.south) / extent.resY;
     int ix = static_cast<int>(grid_x);
     int iy = static_cast<int>(grid_y);
     int ix2 = std::min(ix + 1, grid->width() - 1);
@@ -3392,6 +3422,7 @@ bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid,
         return false;
     }
 
+    // Bilinear interpolation
     double frct_lam = grid_x - ix;
     double frct_phi = grid_y - iy;
     double m10 = frct_lam;
diff --git a/src/grids.hpp b/src/grids.hpp
index 0fd1b7b0..d060fc95 100644
--- a/src/grids.hpp
+++ b/src/grids.hpp
@@ -37,12 +37,14 @@
 NS_PROJ_START
 
 struct ExtentAndRes {
-    double westLon;  // in radian
-    double southLat; // in radian
-    double eastLon;  // in radian
-    double northLat; // in radian
-    double resLon;   // in radian
-    double resLat;   // in radian
+    bool isGeographic; // whether extent and resolutions are in a geographic or
+                       // projected CRS
+    double west;       // in radian for geographic, in CRS units otherwise
+    double south;      // in radian for geographic, in CRS units otherwise
+    double east;       // in radian for geographic, in CRS units otherwise
+    double north;      // in radian for geographic, in CRS units otherwise
+    double resX;       // in radian for geographic, in CRS units otherwise
+    double resY;       // in radian for geographic, in CRS units otherwise
 
     bool fullWorldLongitude() const;
     bool contains(const ExtentAndRes &other) const;
@@ -188,7 +190,7 @@ class PROJ_GCC_DLL GenericShiftGrid : public Grid {
 
     PROJ_FOR_TEST ~GenericShiftGrid() override;
 
-    PROJ_FOR_TEST const GenericShiftGrid *gridAt(double lon, double lat) const;
+    PROJ_FOR_TEST const GenericShiftGrid *gridAt(double x, double y) const;
 
     PROJ_FOR_TEST virtual std::string unit(int sample) const = 0;
 
@@ -228,7 +230,7 @@ class PROJ_GCC_DLL GenericShiftGridSet {
     grids() const {
         return m_grids;
     }
-    PROJ_FOR_TEST const GenericShiftGrid *gridAt(double lon, double lat) const;
+    PROJ_FOR_TEST const GenericShiftGrid *gridAt(double x, double y) const;
 
     PROJ_FOR_TEST virtual void reassign_context(PJ_CONTEXT *ctx);
     PROJ_FOR_TEST virtual bool reopen(PJ_CONTEXT *ctx);
@@ -253,11 +255,9 @@ PJ_LP pj_hgrid_apply(PJ_CONTEXT *ctx, const ListOfHGrids &grids, PJ_LP lp,
 const GenericShiftGrid *pj_find_generic_grid(const ListOfGenericGrids &grids,
                                              const PJ_LP &input,
                                              GenericShiftGridSet *&gridSetOut);
-bool pj_bilinear_interpolation_three_samples(const GenericShiftGrid *grid,
-                                             const PJ_LP &lp, int idx1,
-                                             int idx2, int idx3, double &v1,
-                                             double &v2, double &v3,
-                                             bool &must_retry);
+bool pj_bilinear_interpolation_three_samples(
+    PJ_CONTEXT *ctx, const GenericShiftGrid *grid, const PJ_LP &lp, int idx1,
+    int idx2, int idx3, double &v1, double &v2, double &v3, bool &must_retry);
 
 NS_PROJ_END
 
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 71561fcd..c318a266 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -184,6 +184,7 @@ set(SRC_LIBPROJ_TRANSFORMATIONS
   transformations/molodensky.cpp
   transformations/vgridshift.cpp
   transformations/xyzgridshift.cpp
+  transformations/defmodel.cpp
 )
 
 set(SRC_LIBPROJ_ISO19111
diff --git a/src/pj_list.h b/src/pj_list.h
index b8790b45..cf81ae0e 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -32,6 +32,7 @@ PROJ_HEAD(chamb, "Chamberlin Trimetric")
 PROJ_HEAD(collg, "Collignon")
 PROJ_HEAD(comill, "Compact Miller")
 PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)")
+PROJ_HEAD(defmodel, "Deformation model")
 PROJ_HEAD(deformation, "Kinematic grid shift")
 PROJ_HEAD(denoy, "Denoyer Semi-Elliptical")
 PROJ_HEAD(eck1, "Eckert I")
diff --git a/src/transformations/defmodel.cpp b/src/transformations/defmodel.cpp
new file mode 100644
index 00000000..3d0f2a58
--- /dev/null
+++ b/src/transformations/defmodel.cpp
@@ -0,0 +1,451 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to deformation model
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#define PJ_LIB__
+#define PROJ_COMPILATION
+
+#include "defmodel.hpp"
+#include "filemanager.hpp"
+#include "grids.hpp"
+#include "proj_internal.h"
+
+#include <assert.h>
+
+#include <map>
+#include <memory>
+#include <utility>
+
+PROJ_HEAD(defmodel, "Deformation model");
+
+using namespace DeformationModel;
+
+namespace {
+
+struct Grid : public GridPrototype {
+    PJ_CONTEXT *ctx;
+    const NS_PROJ::GenericShiftGrid *realGrid;
+    mutable bool checkedHorizontal = false;
+    mutable bool checkedVertical = false;
+    mutable int sampleX = 0;
+    mutable int sampleY = 1;
+    mutable int sampleZ = 2;
+
+    Grid(PJ_CONTEXT *ctxIn, const NS_PROJ::GenericShiftGrid *realGridIn)
+        : ctx(ctxIn), realGrid(realGridIn) {
+        minx = realGridIn->extentAndRes().west;
+        miny = realGridIn->extentAndRes().south;
+        resx = realGridIn->extentAndRes().resX;
+        resy = realGridIn->extentAndRes().resY;
+        width = realGridIn->width();
+        height = realGridIn->height();
+    }
+
+    bool checkHorizontal(const std::string &expectedUnit) const {
+        if (!checkedHorizontal) {
+            const auto samplesPerPixel = realGrid->samplesPerPixel();
+            if (samplesPerPixel < 2) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "defmodel: grid %s has not enough samples",
+                       realGrid->name().c_str());
+                return false;
+            }
+            bool foundDescX = false;
+            bool foundDescY = false;
+            bool foundDesc = false;
+            for (int i = 0; i < samplesPerPixel; i++) {
+                const auto desc = realGrid->description(i);
+                if (desc == "east_offset") {
+                    sampleX = i;
+                    foundDescX = true;
+                } else if (desc == "north_offset") {
+                    sampleY = i;
+                    foundDescY = true;
+                }
+                if (!desc.empty()) {
+                    foundDesc = true;
+                }
+            }
+            if (foundDesc && (!foundDescX || !foundDescY)) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "defmodel: grid %s : Found band description, "
+                       "but not the ones expected",
+                       realGrid->name().c_str());
+                return false;
+            }
+            const auto unit = realGrid->unit(sampleX);
+            if (!unit.empty() && unit != expectedUnit) {
+                pj_log(ctx, PJ_LOG_ERROR, "defmodel: grid %s : Only unit=%s "
+                                          "currently handled for this mode",
+                       realGrid->name().c_str(), expectedUnit.c_str());
+                return false;
+            }
+            checkedHorizontal = true;
+        }
+        return true;
+    }
+
+    bool getLonLatOffset(int ix, int iy, double &lonOffsetRadian,
+                         double &latOffsetRadian) const {
+        if (!checkHorizontal(STR_DEGREE)) {
+            return false;
+        }
+        float lonOffsetDeg;
+        float latOffsetDeg;
+        if (!realGrid->valueAt(ix, iy, sampleX, lonOffsetDeg) ||
+            !realGrid->valueAt(ix, iy, sampleY, latOffsetDeg)) {
+            return false;
+        }
+        lonOffsetRadian = lonOffsetDeg * DEG_TO_RAD;
+        latOffsetRadian = latOffsetDeg * DEG_TO_RAD;
+        return true;
+    }
+
+    bool getZOffset(int ix, int iy, double &zOffset) const {
+        if (!checkedVertical) {
+            const auto samplesPerPixel = realGrid->samplesPerPixel();
+            if (samplesPerPixel == 1) {
+                sampleZ = 0;
+            } else if (samplesPerPixel < 3) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "defmodel: grid %s has not enough samples",
+                       realGrid->name().c_str());
+                return false;
+            }
+            bool foundDesc = false;
+            bool foundDescZ = false;
+            for (int i = 0; i < samplesPerPixel; i++) {
+                const auto desc = realGrid->description(i);
+                if (desc == "vertical_offset") {
+                    sampleZ = i;
+                    foundDescZ = true;
+                }
+                if (!desc.empty()) {
+                    foundDesc = true;
+                }
+            }
+            if (foundDesc && !foundDescZ) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "defmodel: grid %s : Found band description, "
+                       "but not the ones expected",
+                       realGrid->name().c_str());
+                return false;
+            }
+            const auto unit = realGrid->unit(sampleZ);
+            if (!unit.empty() && unit != STR_METRE) {
+                pj_log(ctx, PJ_LOG_ERROR,
+                       "defmodel: grid %s : Only unit=metre currently "
+                       "handled for this mode",
+                       realGrid->name().c_str());
+                return false;
+            }
+            checkedVertical = true;
+        }
+        float zOffsetFloat = 0.0f;
+        const bool ret = realGrid->valueAt(ix, iy, sampleZ, zOffsetFloat);
+        zOffset = zOffsetFloat;
+        return ret;
+    }
+
+    bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset,
+                                  double &northingOffset) const {
+        if (!checkHorizontal(STR_METRE)) {
+            return false;
+        }
+        float eastingOffsetFloat = 0.0f;
+        float northingOffsetFloat = 0.0f;
+        const bool ret =
+            realGrid->valueAt(ix, iy, sampleX, eastingOffsetFloat) &&
+            realGrid->valueAt(ix, iy, sampleY, northingOffsetFloat);
+        eastingOffset = eastingOffsetFloat;
+        northingOffset = northingOffsetFloat;
+        return ret;
+    }
+
+    bool getLonLatZOffset(int ix, int iy, double &lonOffsetRadian,
+                          double &latOffsetRadian, double &zOffset) const {
+        return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) &&
+               getZOffset(ix, iy, zOffset);
+    }
+
+    bool getEastingNorthingZOffset(int ix, int iy, double &eastingOffset,
+                                   double &northingOffset,
+                                   double &zOffset) const {
+        return getEastingNorthingOffset(ix, iy, eastingOffset,
+                                        northingOffset) &&
+               getZOffset(ix, iy, zOffset);
+    }
+
+#ifdef DEBUG_DEFMODEL
+    std::string name() const { return realGrid->name(); }
+#endif
+
+  private:
+    Grid(const Grid &) = delete;
+    Grid &operator=(const Grid &) = delete;
+};
+
+struct GridSet : public GridSetPrototype<Grid> {
+
+    PJ_CONTEXT *ctx;
+    std::unique_ptr<NS_PROJ::GenericShiftGridSet> realGridSet;
+    std::map<const NS_PROJ::GenericShiftGrid *, std::unique_ptr<Grid>>
+        mapGrids{};
+
+    GridSet(PJ_CONTEXT *ctxIn,
+            std::unique_ptr<NS_PROJ::GenericShiftGridSet> &&realGridSetIn)
+        : ctx(ctxIn), realGridSet(std::move(realGridSetIn)) {}
+
+    const Grid *gridAt(double x, double y) {
+        const NS_PROJ::GenericShiftGrid *realGrid = realGridSet->gridAt(x, y);
+        if (!realGrid) {
+            return nullptr;
+        }
+        auto iter = mapGrids.find(realGrid);
+        if (iter == mapGrids.end()) {
+            iter = mapGrids
+                       .insert(std::pair<const NS_PROJ::GenericShiftGrid *,
+                                         std::unique_ptr<Grid>>(
+                           realGrid,
+                           std::unique_ptr<Grid>(new Grid(ctx, realGrid))))
+                       .first;
+        }
+        return iter->second.get();
+    }
+
+  private:
+    GridSet(const GridSet &) = delete;
+    GridSet &operator=(const GridSet &) = delete;
+};
+
+struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
+
+    PJ_CONTEXT *ctx;
+    PJ *cart;
+
+    EvaluatorIface(PJ_CONTEXT *ctxIn, PJ *cartIn) : ctx(ctxIn), cart(cartIn) {}
+
+    ~EvaluatorIface() {
+        if (cart)
+            cart->destructor(cart, 0);
+    }
+
+    std::unique_ptr<GridSet> open(const std::string &filename) {
+        auto realGridSet = NS_PROJ::GenericShiftGridSet::open(ctx, filename);
+        if (!realGridSet) {
+            pj_log(ctx, PJ_LOG_ERROR, "defmodel: cannot open %s",
+                   filename.c_str());
+            return nullptr;
+        }
+        return std::unique_ptr<GridSet>(
+            new GridSet(ctx, std::move(realGridSet)));
+    }
+
+    bool isGeographicCRS(const std::string &crsDef) {
+        PJ *P = proj_create(ctx, crsDef.c_str());
+        if (P == nullptr) {
+            return true; // reasonable default value
+        }
+        const auto type = proj_get_type(P);
+        bool ret = (type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
+                    type == PJ_TYPE_GEOGRAPHIC_3D_CRS);
+        proj_destroy(P);
+        return ret;
+    }
+
+    void geographicToGeocentric(double lam, double phi, double height, double a,
+                                double b, double /*es*/, double &X, double &Y,
+                                double &Z) {
+        (void)a;
+        (void)b;
+        assert(cart->a == a);
+        assert(cart->b == b);
+        PJ_LPZ lpz;
+        lpz.lam = lam;
+        lpz.phi = phi;
+        lpz.z = height;
+        PJ_XYZ xyz = cart->fwd3d(lpz, cart);
+        X = xyz.x;
+        Y = xyz.y;
+        Z = xyz.z;
+    }
+
+    void geocentricToGeographic(double X, double Y, double Z, double a,
+                                double b, double /*es*/, double &lam,
+                                double &phi, double &height) {
+        (void)a;
+        (void)b;
+        assert(cart->a == a);
+        assert(cart->b == b);
+        PJ_XYZ xyz;
+        xyz.x = X;
+        xyz.y = Y;
+        xyz.z = Z;
+        PJ_LPZ lpz = cart->inv3d(xyz, cart);
+        lam = lpz.lam;
+        phi = lpz.phi;
+        height = lpz.z;
+    }
+
+#ifdef DEBUG_DEFMODEL
+    void log(const std::string &msg) {
+        pj_log(ctx, PJ_LOG_TRACE, "%s", msg.c_str());
+    }
+#endif
+
+  private:
+    EvaluatorIface(const EvaluatorIface &) = delete;
+    EvaluatorIface &operator=(const EvaluatorIface &) = delete;
+};
+
+struct defmodelData {
+    std::unique_ptr<Evaluator<Grid, GridSet, EvaluatorIface>> evaluator{};
+    EvaluatorIface evaluatorIface;
+
+    explicit defmodelData(PJ_CONTEXT *ctx, PJ *cart)
+        : evaluatorIface(ctx, cart) {}
+
+    defmodelData(const defmodelData &) = delete;
+    defmodelData &operator=(const defmodelData &) = delete;
+};
+
+} // namespace
+
+static PJ *destructor(PJ *P, int errlev) {
+    if (nullptr == P)
+        return nullptr;
+
+    auto Q = static_cast<struct defmodelData *>(P->opaque);
+    delete Q;
+    P->opaque = nullptr;
+
+    return pj_default_destructor(P, errlev);
+}
+
+static PJ_COORD forward_4d(PJ_COORD in, 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();
+    }
+
+    return out;
+}
+
+static PJ_COORD reverse_4d(PJ_COORD in, 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();
+    }
+
+    return out;
+}
+
+// Function called by proj_assign_context() when a new context is assigned to
+// an existing PJ object. Mostly to deal with objects being passed between
+// threads.
+static void reassign_context(PJ *P, PJ_CONTEXT *ctx) {
+    auto *Q = (struct defmodelData *)P->opaque;
+    if (Q->evaluatorIface.ctx != ctx) {
+        Q->evaluator->clearGridCache();
+        Q->evaluatorIface.ctx = ctx;
+    }
+}
+
+PJ *TRANSFORMATION(defmodel, 1) {
+    // Pass a dummy ellipsoid definition that will be overridden just afterwards
+    auto cart = proj_create(P->ctx, "+proj=cart +a=1");
+    if (cart == nullptr)
+        return destructor(P, ENOMEM);
+
+    /* inherit ellipsoid definition from P to Q->cart */
+    pj_inherit_ellipsoid_def(P, cart);
+
+    auto Q = new defmodelData(P->ctx, cart);
+    P->opaque = (void *)Q;
+    P->destructor = destructor;
+    P->reassign_context = reassign_context;
+
+    const char *model = pj_param(P->ctx, P->params, "smodel").s;
+    if (!model) {
+        proj_log_error(P, "defmodel: +model= should be specified.");
+        return destructor(P, PJD_ERR_NO_ARGS);
+    }
+
+    auto file = NS_PROJ::FileManager::open_resource_file(P->ctx, model);
+    if (nullptr == file) {
+        proj_log_error(P, "defmodel: Cannot open %s", model);
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+    file->seek(0, SEEK_END);
+    unsigned long long size = file->tell();
+    // Arbitrary threshold to avoid ingesting an arbitrarily large JSON file,
+    // that could be a denial of service risk. 10 MB should be sufficiently
+    // large for any valid use !
+    if (size > 10 * 1024 * 1024) {
+        proj_log_error(P, "defmodel: File %s too large", model);
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+    file->seek(0);
+    std::string jsonStr;
+    jsonStr.resize(static_cast<size_t>(size));
+    if (file->read(&jsonStr[0], jsonStr.size()) != jsonStr.size()) {
+        proj_log_error(P, "defmodel: Cannot read %s", model);
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+
+    try {
+        Q->evaluator.reset(new Evaluator<Grid, GridSet, EvaluatorIface>(
+            MasterFile::parse(jsonStr), Q->evaluatorIface, P->a, P->b));
+    } catch (const std::exception &e) {
+        proj_log_error(P, "defmodel: invalid model: %s", e.what());
+        return destructor(P, PJD_ERR_INVALID_ARG);
+    }
+
+    P->fwd4d = forward_4d;
+    P->inv4d = reverse_4d;
+
+    if (Q->evaluator->isGeographicCRS()) {
+        P->left = PJ_IO_UNITS_RADIANS;
+        P->right = PJ_IO_UNITS_RADIANS;
+    } else {
+        P->left = PJ_IO_UNITS_PROJECTED;
+        P->right = PJ_IO_UNITS_PROJECTED;
+    }
+
+    return P;
+}
diff --git a/src/transformations/defmodel.hpp b/src/transformations/defmodel.hpp
new file mode 100644
index 00000000..9b28a7d8
--- /dev/null
+++ b/src/transformations/defmodel.hpp
@@ -0,0 +1,630 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to deformation model
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+/** This file implements the gridded deformation model proposol of
+ * https://docs.google.com/document/d/1wiyrAmzqh8MZlzHSp3wf594Ob_M1LeFtDA5swuzvLZY
+ * It is written in a generic way, independent of the rest of PROJ
+ * infrastructure.
+ *
+ * Verbose debugging info can be turned on by setting the DEBUG_DEFMODEL macro
+ */
+
+#ifndef DEFMODEL_HPP
+#define DEFMODEL_HPP
+
+#ifdef PROJ_COMPILATION
+#include "proj/internal/include_nlohmann_json.hpp"
+#else
+#include "nlohmann/json.hpp"
+#endif
+
+#include <algorithm>
+#include <cmath>
+#include <exception>
+#include <limits>
+#include <memory>
+#include <string>
+#include <vector>
+
+#ifndef DEFORMATON_MODEL_NAMESPACE
+#define DEFORMATON_MODEL_NAMESPACE DeformationModel
+#endif
+
+#include "defmodel_exceptions.hpp"
+
+namespace DEFORMATON_MODEL_NAMESPACE {
+
+using json = nlohmann::json;
+
+// ---------------------------------------------------------------------------
+
+/** Spatial extent as a bounding box. */
+class SpatialExtent {
+  public:
+    /** Parse the provided object as an extent.
+    *
+    * @throws ParsingException
+    */
+    static SpatialExtent parse(const json &j);
+
+    double minx() const { return mMinx; }
+    double miny() const { return mMiny; }
+    double maxx() const { return mMaxx; }
+    double maxy() const { return mMaxy; }
+
+    double minxNormalized(bool bIsGeographic) const {
+        return bIsGeographic ? mMinxRad : mMinx;
+    }
+    double minyNormalized(bool bIsGeographic) const {
+        return bIsGeographic ? mMinyRad : mMiny;
+    }
+    double maxxNormalized(bool bIsGeographic) const {
+        return bIsGeographic ? mMaxxRad : mMaxx;
+    }
+    double maxyNormalized(bool bIsGeographic) const {
+        return bIsGeographic ? mMaxyRad : mMaxy;
+    }
+
+  protected:
+    friend class MasterFile;
+    friend class Component;
+    SpatialExtent() = default;
+
+  private:
+    double mMinx = std::numeric_limits<double>::quiet_NaN();
+    double mMiny = std::numeric_limits<double>::quiet_NaN();
+    double mMaxx = std::numeric_limits<double>::quiet_NaN();
+    double mMaxy = std::numeric_limits<double>::quiet_NaN();
+    double mMinxRad = std::numeric_limits<double>::quiet_NaN();
+    double mMinyRad = std::numeric_limits<double>::quiet_NaN();
+    double mMaxxRad = std::numeric_limits<double>::quiet_NaN();
+    double mMaxyRad = std::numeric_limits<double>::quiet_NaN();
+};
+
+// ---------------------------------------------------------------------------
+
+/** Epoch */
+class Epoch {
+  public:
+    /** Constructor from a ISO 8601 date-time. May throw ParsingException */
+    explicit Epoch(const std::string &dt = std::string());
+
+    /** Return ISO 8601 date-time */
+    const std::string &toString() const { return mDt; }
+
+    /** Return decimal year */
+    double toDecimalYear() const;
+
+  private:
+    std::string mDt{};
+    double mDecimalYear = 0;
+};
+
+// ---------------------------------------------------------------------------
+
+/** Component of a deformation model. */
+class Component {
+  public:
+    /** Parse the provided object as a component.
+    *
+    * @throws ParsingException
+    */
+    static Component parse(const json &j);
+
+    /** Get a text description of the component. */
+    const std::string &description() const { return mDescription; }
+
+    /** Get the region within which the component is defined. Outside this
+     * region the component evaluates to 0. */
+    const SpatialExtent &extent() const { return mSpatialExtent; }
+
+    /** Get the displacement parameters defined by the component, one of
+     * "none", "horizontal", "vertical", and "3d".  The "none" option allows
+     *  for a component which defines uncertainty with different grids to those
+     *  defining displacement. */
+    const std::string &displacementType() const { return mDisplacementType; }
+
+    /** Get the uncertainty parameters defined by the component,
+     * one of "none", "horizontal", "vertical", "3d". */
+    const std::string &uncertaintyType() const { return mUncertaintyType; }
+
+    /** Get the horizontal uncertainty to use if it is not defined explicitly
+     * in the spatial model. */
+    double horizontalUncertainty() const { return mHorizontalUncertainty; }
+
+    /** Get the vertical uncertainty to use if it is not defined explicitly in
+     * the spatial model. */
+    double verticalUncertainty() const { return mVerticalUncertainty; }
+
+    struct SpatialModel {
+        /** Specifies the type of the spatial model data file.  Initially it
+         * is proposed that only "GeoTIFF" is supported. */
+        std::string type{};
+
+        /** How values in model should be interpolated. This proposal will
+         * support "bilinear" and "geocentric_bilinear".  */
+        std::string interpolationMethod{};
+
+        /** Specifies location of the spatial model GeoTIFF file relative to
+         * the master JSON file. */
+        std::string filename{};
+
+        /** A hex encoded MD5 checksum of the grid file that can be used to
+         * validate that it is the correct version of the file. */
+        std::string md5Checksum{};
+    };
+
+    /** Get the spatial model. */
+    const SpatialModel &spatialModel() const { return mSpatialModel; }
+
+    /** Generic type for a type functon */
+    struct TimeFunction {
+        std::string type{};
+
+        virtual ~TimeFunction();
+
+        virtual double evaluateAt(double dt) const = 0;
+
+      protected:
+        TimeFunction() = default;
+    };
+    struct ConstantTimeFunction : public TimeFunction {
+
+        virtual double evaluateAt(double dt) const override;
+    };
+    struct VelocityTimeFunction : public TimeFunction {
+        /** Date/time at which the velocity function is zero. */
+        Epoch referenceEpoch{};
+
+        virtual double evaluateAt(double dt) const override;
+    };
+
+    struct StepTimeFunction : public TimeFunction {
+        /** Epoch at which the step function transitions from 0 to 1. */
+        Epoch stepEpoch{};
+
+        virtual double evaluateAt(double dt) const override;
+    };
+
+    struct ReverseStepTimeFunction : public TimeFunction {
+        /** Epoch at which the reverse step function transitions from 1. to 0 */
+        Epoch stepEpoch{};
+
+        virtual double evaluateAt(double dt) const override;
+    };
+
+    struct PiecewiseTimeFunction : public TimeFunction {
+        /** One of "zero", "constant", and "linear", defines the behaviour of
+         * the function before the first defined epoch */
+        std::string beforeFirst{};
+
+        /** One of "zero", "constant", and "linear", defines the behaviour of
+         * the function after the last defined epoch */
+        std::string afterLast{};
+
+        struct EpochScaleFactorTuple {
+            /** Defines the date/time of the data point. */
+            Epoch epoch{};
+
+            /** Function value at the above epoch */
+            double scaleFactor = std::numeric_limits<double>::quiet_NaN();
+        };
+
+        /** A sorted array data points each defined by two elements.
+         * The array is sorted in order of increasing epoch.
+         * Note: where the time function includes a step it is represented by
+         * two consecutive data points with the same epoch.  The first defines
+         * the scale factor that applies before the epoch and the second the
+         * scale factor that applies after the epoch. */
+        std::vector<EpochScaleFactorTuple> model{};
+
+        virtual double evaluateAt(double dt) const override;
+    };
+
+    struct ExponentialTimeFunction : public TimeFunction {
+        /** The date/time at which the exponential decay starts. */
+        Epoch referenceEpoch{};
+
+        /** The date/time at which the exponential decay ends. */
+        Epoch endEpoch{};
+
+        /** The relaxation constant in years. */
+        double relaxationConstant = std::numeric_limits<double>::quiet_NaN();
+
+        /** The scale factor that applies before the reference epoch. */
+        double beforeScaleFactor = std::numeric_limits<double>::quiet_NaN();
+
+        /** Initial scale factor. */
+        double initialScaleFactor = std::numeric_limits<double>::quiet_NaN();
+
+        /** The scale factor the exponential function approaches. */
+        double finalScaleFactor = std::numeric_limits<double>::quiet_NaN();
+
+        virtual double evaluateAt(double dt) const override;
+    };
+
+    /** Get the time function. */
+    const TimeFunction *timeFunction() const { return mTimeFunction.get(); }
+
+  private:
+    Component() = default;
+
+    std::string mDescription{};
+    SpatialExtent mSpatialExtent{};
+    std::string mDisplacementType{};
+    std::string mUncertaintyType{};
+    double mHorizontalUncertainty = std::numeric_limits<double>::quiet_NaN();
+    double mVerticalUncertainty = std::numeric_limits<double>::quiet_NaN();
+    SpatialModel mSpatialModel{};
+    std::unique_ptr<TimeFunction> mTimeFunction{};
+};
+
+Component::TimeFunction::~TimeFunction() = default;
+
+// ---------------------------------------------------------------------------
+
+/** Master file of a deformation model. */
+class MasterFile {
+  public:
+    /** Parse the provided serialized JSON content and return an object.
+     *
+     * @throws ParsingException
+     */
+    static std::unique_ptr<MasterFile> parse(const std::string &text);
+
+    /** Get file type. Should always be "deformation_model_master_file" */
+    const std::string &fileType() const { return mFileType; }
+
+    /** Get the version of the format. At time of writing, only "1.0" is known
+     */
+    const std::string &formatVersion() const { return mFormatVersion; }
+
+    /** Get brief descriptive name of the deformation model. */
+    const std::string &name() const { return mName; }
+
+    /** Get a string identifying the version of the deformation model.
+     * The format for specifying version is defined by the agency
+     * responsible for the deformation model. */
+    const std::string &version() const { return mVersion; }
+
+    /** Get a string identifying the license of the file.
+     * e.g "Create Commons Attribution 4.0 International" */
+    const std::string &license() const { return mLicense; }
+
+    /** Get a text description of the model. Intended to be longer than name()
+     */
+    const std::string &description() const { return mDescription; }
+
+    /** Get a text description of the model. Intended to be longer than name()
+     */
+    const std::string &publicationDate() const { return mPublicationDate; }
+
+    /** Basic information on the agency responsible for the model. */
+    struct Authority {
+        std::string name{};
+        std::string url{};
+        std::string address{};
+        std::string email{};
+    };
+
+    /** Get basic information on the agency responsible for the model. */
+    const Authority &authority() const { return mAuthority; }
+
+    /** Hyperlink related to the model. */
+    struct Link {
+        /** URL holding the information */
+        std::string href{};
+
+        /** Relationship to the dataset. e.g. "about", "source", "license",
+         * "metadata" */
+        std::string rel{};
+
+        /** Mime type */
+        std::string type{};
+
+        /** Description of the link */
+        std::string title{};
+    };
+
+    /** Get links to related information. */
+    const std::vector<Link> links() const { return mLinks; }
+
+    /** Get a string identifying the source CRS. That is the coordinate
+     * reference system to which the deformation model applies. Typically
+     * "EPSG:XXXX" */
+    const std::string &sourceCRS() const { return mSourceCRS; }
+
+    /** Get a string identifying the target CRS. That is, for a time
+     * dependent coordinate transformation, the coordinate reference
+     * system resulting from applying the deformation.
+     * Typically "EPSG:XXXX" */
+    const std::string &targetCRS() const { return mTargetCRS; }
+
+    /** Get a string identifying the definition CRS. That is, the
+     * coordinate reference system used to define the component spatial
+     * models. Typically "EPSG:XXXX" */
+    const std::string &definitionCRS() const { return mDefinitionCRS; }
+
+    /** Get the nominal reference epoch of the deformation model. Formated
+     * as a ISO-8601 date-time. This is not necessarily used to calculate
+     * the deformation model - each component defines its own time function. */
+    const std::string &referenceEpoch() const { return mReferenceEpoch; }
+
+    /** Get the epoch at which the uncertainties of the deformation model
+     * are calculated. Formated as a ISO-8601 date-time. */
+    const std::string &uncertaintyReferenceEpoch() const {
+        return mUncertaintyReferenceEpoch;
+    }
+
+    /** Unit of horizontal offsets. Only "metre" and "degree" are supported. */
+    const std::string &horizontalOffsetUnit() const {
+        return mHorizontalOffsetUnit;
+    }
+
+    /** Unit of vertical offsets. Only "metre" is supported. */
+    const std::string &verticalOffsetUnit() const {
+        return mVerticalOffsetUnit;
+    }
+
+    /** Type of horizontal uncertainty. e.g "circular 95% confidence limit" */
+    const std::string &horizontalUncertaintyType() const {
+        return mHorizontalUncertaintyType;
+    }
+
+    /** Unit of horizontal uncertainty. Only "metre" is supported. */
+    const std::string &horizontalUncertaintyUnit() const {
+        return mHorizontalUncertaintyUnit;
+    }
+
+    /** Type of vertical uncertainty. e.g "circular 95% confidence limit" */
+    const std::string &verticalUncertaintyType() const {
+        return mVerticalUncertaintyType;
+    }
+
+    /** Unit of vertical uncertainty. Only "metre" is supported. */
+    const std::string &verticalUncertaintyUnit() const {
+        return mVerticalUncertaintyUnit;
+    }
+
+    /** Defines how the horizontal offsets are applied to geographic
+     * coordinates. Only "addition" and "geocentric" are supported */
+    const std::string &horizontalOffsetMethod() const {
+        return mHorizontalOffsetMethod;
+    }
+
+    /** Get the region within which the deformation model is defined.
+     * It cannot be calculated outside this region */
+    const SpatialExtent &extent() const { return mSpatialExtent; }
+
+    /** Defines the range of times for which the model is valid, specified
+     * by a first and a last value. The deformation model is undefined for
+     * dates outside this range. */
+    struct TimeExtent {
+        Epoch first{};
+        Epoch last{};
+    };
+
+    /** Get the range of times for which the model is valid. */
+    const TimeExtent &timeExtent() const { return mTimeExtent; }
+
+    /** Get an array of the components comprising the deformation model. */
+    const std::vector<Component> &components() const { return mComponents; }
+
+  private:
+    MasterFile() = default;
+
+    std::string mFileType{};
+    std::string mFormatVersion{};
+    std::string mName{};
+    std::string mVersion{};
+    std::string mLicense{};
+    std::string mDescription{};
+    std::string mPublicationDate{};
+    Authority mAuthority{};
+    std::vector<Link> mLinks{};
+    std::string mSourceCRS{};
+    std::string mTargetCRS{};
+    std::string mDefinitionCRS{};
+    std::string mReferenceEpoch{};
+    std::string mUncertaintyReferenceEpoch{};
+    std::string mHorizontalOffsetUnit{};
+    std::string mVerticalOffsetUnit{};
+    std::string mHorizontalUncertaintyType{};
+    std::string mHorizontalUncertaintyUnit{};
+    std::string mVerticalUncertaintyType{};
+    std::string mVerticalUncertaintyUnit{};
+    std::string mHorizontalOffsetMethod{};
+    SpatialExtent mSpatialExtent{};
+    TimeExtent mTimeExtent{};
+    std::vector<Component> mComponents{};
+};
+
+// ---------------------------------------------------------------------------
+
+/** Prototype for a Grid used by GridSet. Intended to be implemented
+ * by user code */
+struct GridPrototype {
+    double minx = 0;
+    double miny = 0;
+    double resx = 0;
+    double resy = 0;
+    int width = 0;
+    int height = 0;
+
+    bool getLonLatOffset(int /*ix*/, int /*iy*/, double & /*lonOffsetRadian*/,
+                         double & /*latOffsetRadian*/) const {
+        throw UnimplementedException("getLonLatOffset unimplemented");
+    }
+
+    bool getZOffset(int /*ix*/, int /*iy*/, double & /*zOffset*/) const {
+        throw UnimplementedException("getZOffset unimplemented");
+    }
+
+    bool getEastingNorthingOffset(int /*ix*/, int /*iy*/,
+                                  double & /*eastingOffset*/,
+                                  double & /*northingOffset*/) const {
+        throw UnimplementedException("getEastingNorthingOffset unimplemented");
+    }
+
+    bool getLonLatZOffset(int /*ix*/, int /*iy*/, double & /*lonOffsetRadian*/,
+                          double & /*latOffsetRadian*/,
+                          double & /*zOffset*/) const {
+        throw UnimplementedException("getLonLatZOffset unimplemented");
+#if 0
+        return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) &&
+               getZOffset(ix, iy, zOffset);
+#endif
+    }
+
+    bool getEastingNorthingZOffset(int /*ix*/, int /*iy*/,
+                                   double & /*eastingOffset*/,
+                                   double & /*northingOffset*/,
+                                   double & /*zOffset*/) const {
+        throw UnimplementedException("getEastingNorthingOffset unimplemented");
+#if 0
+        return getEastingNorthingOffset(ix, iy, eastingOffset,
+                                        northingOffset) &&
+               getZOffset(ix, iy, zOffset);
+#endif
+    }
+
+#ifdef DEBUG_DEFMODEL
+    std::string name() const {
+        throw UnimplementedException("name() unimplemented");
+    }
+#endif
+};
+
+// ---------------------------------------------------------------------------
+
+/** Prototype for a GridSet used by EvaluatorIface. Intended to be implemented
+ * by user code */
+template <class Grid = GridPrototype> struct GridSetPrototype {
+    // The return pointer should remain "stable" over time for a given grid
+    // of a GridSet.
+    const Grid *gridAt(double /*x */, double /* y */) {
+        throw UnimplementedException("gridAt unimplemented");
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+/** Prototype for a EvaluatorIface used by Evaluator. Intended to be implemented
+ * by user code */
+template <class Grid = GridPrototype, class GridSet = GridSetPrototype<>>
+struct EvaluatorIfacePrototype {
+
+    std::unique_ptr<GridSet> open(const std::string & /* filename*/) {
+        throw UnimplementedException("open unimplemented");
+    }
+
+    void geographicToGeocentric(double /* lam */, double /* phi */,
+                                double /* height*/, double /* a */,
+                                double /* b */, double /*es*/, double & /* X */,
+                                double & /* Y */, double & /* Z */) {
+        throw UnimplementedException("geographicToGeocentric unimplemented");
+    }
+
+    void geocentricToGeographic(double /* X */, double /* Y */, double /* Z */,
+                                double /* a */, double /* b */, double /*es*/,
+                                double & /* lam */, double & /* phi */,
+                                double & /* height*/) {
+        throw UnimplementedException("geocentricToGeographic unimplemented");
+    }
+
+    bool isGeographicCRS(const std::string & /* crsDef */) {
+        throw UnimplementedException("isGeographicCRS unimplemented");
+    }
+
+#ifdef DEBUG_DEFMODEL
+    void log(const std::string & /* msg */) {
+        throw UnimplementedException("log unimplemented");
+    }
+#endif
+};
+
+// ---------------------------------------------------------------------------
+
+/** Internal class to offer caching services over a Component */
+template <class Grid, class GridSet> struct ComponentEx;
+
+// ---------------------------------------------------------------------------
+
+/** Class to evaluate the transformation of a coordinate */
+template <class Grid = GridPrototype, class GridSet = GridSetPrototype<>,
+          class EvaluatorIface = EvaluatorIfacePrototype<>>
+class Evaluator {
+  public:
+    /** Constructor. May throw EvaluatorException */
+    explicit Evaluator(std::unique_ptr<MasterFile> &&model,
+                       EvaluatorIface &iface, double a, double b);
+
+    /** Evaluate displacement of a position given by (x,y,z,t) and
+     * return it in (x_out,y_out_,z_out).
+     * For geographic CRS (only supported at that time), x must be a
+     * longitude, and y a latitude.
+     */
+    bool forward(EvaluatorIface &iface, double x, double y, double z, double t,
+                 double &x_out, double &y_out, double &z_out) {
+        return forward(iface, x, y, z, t, false, x_out, y_out, z_out);
+    }
+
+    /** Apply inverse transformation. */
+    bool inverse(EvaluatorIface &iface, double x, double y, double z, double t,
+                 double &x_out, double &y_out, double &z_out);
+
+    /** Clear grid cache */
+    void clearGridCache();
+
+    /** Return whether the definition CRS is a geographic CRS */
+    bool isGeographicCRS() const { return mIsGeographicCRS; }
+
+  private:
+    std::unique_ptr<MasterFile> mModel;
+    const double mA;
+    const double mB;
+    const double mEs;
+    const bool mIsHorizontalUnitDegree; /* degree vs metre */
+    const bool mIsAddition;             /* addition vs geocentric */
+    const bool mIsGeographicCRS;
+
+    bool forward(EvaluatorIface &iface, double x, double y, double z, double t,
+                 bool forInverseComputation, double &x_out, double &y_out,
+                 double &z_out);
+
+    std::vector<std::unique_ptr<ComponentEx<Grid, GridSet>>> mComponents{};
+};
+
+// ---------------------------------------------------------------------------
+
+} // namespace DeformationModel
+
+// ---------------------------------------------------------------------------
+
+#include "defmodel_impl.hpp"
+
+#endif // DEFMODEL_HPP
diff --git a/src/transformations/defmodel_exceptions.hpp b/src/transformations/defmodel_exceptions.hpp
new file mode 100644
index 00000000..8addf60f
--- /dev/null
+++ b/src/transformations/defmodel_exceptions.hpp
@@ -0,0 +1,81 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to deformation model
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef DEFORMATON_MODEL_NAMESPACE
+#error "Should be included only by defmodel.hpp"
+#endif
+
+#include <exception>
+
+namespace DEFORMATON_MODEL_NAMESPACE {
+
+// ---------------------------------------------------------------------------
+
+/** Parsing exception. */
+class ParsingException : public std::exception {
+  public:
+    ParsingException(const std::string &msg) : msg_(msg) {}
+    const char *what() const noexcept override;
+
+  private:
+    std::string msg_;
+};
+
+const char *ParsingException::what() const noexcept { return msg_.c_str(); }
+
+// ---------------------------------------------------------------------------
+
+class UnimplementedException : public std::exception {
+  public:
+    UnimplementedException(const std::string &msg) : msg_(msg) {}
+    const char *what() const noexcept override;
+
+  private:
+    std::string msg_;
+};
+
+const char *UnimplementedException::what() const noexcept {
+    return msg_.c_str();
+}
+
+// ---------------------------------------------------------------------------
+
+/** Evaluator exception. */
+class EvaluatorException : public std::exception {
+  public:
+    EvaluatorException(const std::string &msg) : msg_(msg) {}
+    const char *what() const noexcept override;
+
+  private:
+    std::string msg_;
+};
+
+const char *EvaluatorException::what() const noexcept { return msg_.c_str(); }
+
+// ---------------------------------------------------------------------------
+
+} // namespace DEFORMATON_MODEL_NAMESPACE
diff --git a/src/transformations/defmodel_impl.hpp b/src/transformations/defmodel_impl.hpp
new file mode 100644
index 00000000..a15137d7
--- /dev/null
+++ b/src/transformations/defmodel_impl.hpp
@@ -0,0 +1,1265 @@
+/******************************************************************************
+ * Project:  PROJ
+ * Purpose:  Functionality related to deformation model
+ * Author:   Even Rouault, <even.rouault at spatialys.com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#ifndef DEFORMATON_MODEL_NAMESPACE
+#error "Should be included only by defmodel.hpp"
+#endif
+
+namespace DEFORMATON_MODEL_NAMESPACE {
+
+// ---------------------------------------------------------------------------
+
+static const std::string STR_DEGREE("degree");
+static const std::string STR_METRE("metre");
+
+static const std::string STR_ADDITION("addition");
+static const std::string STR_GEOCENTRIC("geocentric");
+
+static const std::string STR_BILINEAR("bilinear");
+static const std::string STR_GEOCENTRIC_BILINEAR("geocentric_bilinear");
+
+static const std::string STR_NONE("none");
+static const std::string STR_HORIZONTAL("horizontal");
+static const std::string STR_VERTICAL("vertical");
+static const std::string STR_3D("3d");
+
+constexpr double DEFMODEL_PI = 3.14159265358979323846;
+constexpr double DEG_TO_RAD_CONSTANT = 3.14159265358979323846 / 180.;
+inline constexpr double DegToRad(double d) { return d * DEG_TO_RAD_CONSTANT; }
+
+// ---------------------------------------------------------------------------
+
+enum class DisplacementType { NONE, HORIZONTAL, VERTICAL, THREE_D };
+
+// ---------------------------------------------------------------------------
+
+/** Internal class to offer caching services over a Grid */
+template <class Grid> struct GridEx {
+    const Grid *grid;
+
+    bool smallResx; // lesser than one degree
+
+    double sinhalfresx;
+    double coshalfresx;
+
+    double sinresy;
+    double cosresy;
+
+    int last_ix0 = -1;
+    int last_iy0 = -1;
+    double dX00 = 0;
+    double dY00 = 0;
+    double dZ00 = 0;
+    double dX01 = 0;
+    double dY01 = 0;
+    double dZ01 = 0;
+    double dX10 = 0;
+    double dY10 = 0;
+    double dZ10 = 0;
+    double dX11 = 0;
+    double dY11 = 0;
+    double dZ11 = 0;
+    double sinphi0 = 0;
+    double cosphi0 = 0;
+    double sinphi1 = 0;
+    double cosphi1 = 0;
+
+    explicit GridEx(const Grid *gridIn)
+        : grid(gridIn), smallResx(grid->resx < DegToRad(1)),
+          sinhalfresx(sin(grid->resx / 2)), coshalfresx(cos(grid->resx / 2)),
+          sinresy(sin(grid->resy)), cosresy(cos(grid->resy)) {}
+
+    // Return geocentric offset (dX, dY, dZ) relative to a point
+    // where x0 = -resx / 2
+    inline void getBilinearGeocentric(int ix0, int iy0, double de00,
+                                      double dn00, double de01, double dn01,
+                                      double de10, double dn10, double de11,
+                                      double dn11, double m00, double m01,
+                                      double m10, double m11, double &dX,
+                                      double &dY, double &dZ) {
+        // If interpolating in the same cell as before, then we can skip
+        // the recomputation of dXij, dYij and dZij
+        if (ix0 != last_ix0 || iy0 != last_iy0) {
+
+            last_ix0 = ix0;
+            if (iy0 != last_iy0) {
+                const double y0 = grid->miny + iy0 * grid->resy;
+                sinphi0 = sin(y0);
+                cosphi0 = cos(y0);
+
+                // Use trigonometric formulas to avoid new calls to sin/cos
+                sinphi1 =
+                    /* sin(y0+grid->resyRad) = */ sinphi0 * cosresy +
+                    cosphi0 * sinresy;
+                cosphi1 =
+                    /* cos(y0+grid->resyRad) = */ cosphi0 * cosresy -
+                    sinphi0 * sinresy;
+
+                last_iy0 = iy0;
+            }
+
+            // Using "traditional" formulas to convert from easting, northing
+            // offsets to geocentric offsets
+            const double sinlam00 = -sinhalfresx;
+            const double coslam00 = coshalfresx;
+            const double sinphi00 = sinphi0;
+            const double cosphi00 = cosphi0;
+            const double dn00sinphi00 = dn00 * sinphi00;
+            dX00 = -de00 * sinlam00 - dn00sinphi00 * coslam00;
+            dY00 = de00 * coslam00 - dn00sinphi00 * sinlam00;
+            dZ00 = dn00 * cosphi00;
+
+            const double sinlam01 = -sinhalfresx;
+            const double coslam01 = coshalfresx;
+            const double sinphi01 = sinphi1;
+            const double cosphi01 = cosphi1;
+            const double dn01sinphi01 = dn01 * sinphi01;
+            dX01 = -de01 * sinlam01 - dn01sinphi01 * coslam01;
+            dY01 = de01 * coslam01 - dn01sinphi01 * sinlam01;
+            dZ01 = dn01 * cosphi01;
+
+            const double sinlam10 = sinhalfresx;
+            const double coslam10 = coshalfresx;
+            const double sinphi10 = sinphi0;
+            const double cosphi10 = cosphi0;
+            const double dn10sinphi10 = dn10 * sinphi10;
+            dX10 = -de10 * sinlam10 - dn10sinphi10 * coslam10;
+            dY10 = de10 * coslam10 - dn10sinphi10 * sinlam10;
+            dZ10 = dn10 * cosphi10;
+
+            const double sinlam11 = sinhalfresx;
+            const double coslam11 = coshalfresx;
+            const double sinphi11 = sinphi1;
+            const double cosphi11 = cosphi1;
+            const double dn11sinphi11 = dn11 * sinphi11;
+            dX11 = -de11 * sinlam11 - dn11sinphi11 * coslam11;
+            dY11 = de11 * coslam11 - dn11sinphi11 * sinlam11;
+            dZ11 = dn11 * cosphi11;
+        }
+
+        dX = m00 * dX00 + m01 * dX01 + m10 * dX10 + m11 * dX11;
+        dY = m00 * dY00 + m01 * dY01 + m10 * dY10 + m11 * dY11;
+        dZ = m00 * dZ00 + m01 * dZ01 + m10 * dZ10 + m11 * dZ11;
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+/** Internal class to offer caching services over a Component */
+template <class Grid, class GridSet> struct ComponentEx {
+    const Component &component;
+
+    const bool isBilinearInterpolation; /* bilinear vs geocentric_bilinear */
+
+    const DisplacementType displacementType;
+
+    // Cache
+    std::unique_ptr<GridSet> gridSet{};
+    std::map<const Grid *, GridEx<Grid>> mapGrids{};
+
+  private:
+    mutable double mCachedDt = 0;
+    mutable double mCachedValue = 0;
+
+    static DisplacementType getDisplacementType(const std::string &s) {
+        if (s == STR_HORIZONTAL)
+            return DisplacementType::HORIZONTAL;
+        if (s == STR_VERTICAL)
+            return DisplacementType::VERTICAL;
+        if (s == STR_3D)
+            return DisplacementType::THREE_D;
+        return DisplacementType::NONE;
+    }
+
+  public:
+    explicit ComponentEx(const Component &componentIn)
+        : component(componentIn),
+          isBilinearInterpolation(
+              componentIn.spatialModel().interpolationMethod == STR_BILINEAR),
+          displacementType(getDisplacementType(component.displacementType())) {}
+
+    double evaluateAt(double dt) const {
+        if (dt == mCachedDt)
+            return mCachedValue;
+        mCachedDt = dt;
+        mCachedValue = component.timeFunction()->evaluateAt(dt);
+        return mCachedValue;
+    }
+
+    void clearGridCache() {
+        gridSet.reset();
+        mapGrids.clear();
+    }
+};
+
+// ---------------------------------------------------------------------------
+
+/** Converts a ISO8601 date-time string, formatted as "YYYY-MM-DDTHH:MM:SSZ",
+ * into a decimal year.
+ * Leap years are taken into account, but not leap seconds.
+ */
+static double ISO8601ToDecimalYear(const std::string &dt) {
+    int year, month, day, hour, min, sec;
+    if (sscanf(dt.c_str(), "%04d-%02d-%02dT%02d:%02d:%02dZ", &year, &month,
+               &day, &hour, &min, &sec) != 6 ||
+        year < 1582 || // Start of Gregorian calendar
+        month < 1 || month > 12 || day < 1 || day > 31 || hour < 0 ||
+        hour >= 24 || min < 0 || min >= 60 || sec < 0 || sec >= 61) {
+        throw ParsingException("Wrong formatting / invalid date-time for " +
+                               dt);
+    }
+    const bool isLeapYear =
+        (((year % 4) == 0 && (year % 100) != 0) || (year % 400) == 0);
+    // Given the intended use, we omit leap seconds...
+    const int month_table[2][12] = {
+        {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
+        {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}};
+    int dayInYear = day - 1;
+    for (int m = 1; m < month; m++) {
+        dayInYear += month_table[isLeapYear ? 1 : 0][m - 1];
+    }
+    if (day > month_table[isLeapYear ? 1 : 0][month - 1]) {
+        throw ParsingException("Wrong formatting / invalid date-time for " +
+                               dt);
+    }
+    return year +
+           (dayInYear * 86400 + hour * 3600 + min * 60 + sec) /
+               (isLeapYear ? 86400. * 366 : 86400. * 365);
+}
+
+// ---------------------------------------------------------------------------
+
+Epoch::Epoch(const std::string &dt) : mDt(dt) {
+    if (!dt.empty()) {
+        mDecimalYear = ISO8601ToDecimalYear(dt);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+double Epoch::toDecimalYear() const { return mDecimalYear; }
+
+// ---------------------------------------------------------------------------
+
+static std::string getString(const json &j, const char *key, bool optional) {
+    if (!j.contains(key)) {
+        if (optional) {
+            return std::string();
+        }
+        throw ParsingException(std::string("Missing \"") + key + "\" key");
+    }
+    const json v = j[key];
+    if (!v.is_string()) {
+        throw ParsingException(std::string("The value of \"") + key +
+                               "\" should be a string");
+    }
+    return v.get<std::string>();
+}
+
+static std::string getReqString(const json &j, const char *key) {
+    return getString(j, key, false);
+}
+
+static std::string getOptString(const json &j, const char *key) {
+    return getString(j, key, true);
+}
+
+// ---------------------------------------------------------------------------
+
+static double getDouble(const json &j, const char *key, bool optional) {
+    if (!j.contains(key)) {
+        if (optional) {
+            return std::numeric_limits<double>::quiet_NaN();
+        }
+        throw ParsingException(std::string("Missing \"") + key + "\" key");
+    }
+    const json v = j[key];
+    if (!v.is_number()) {
+        throw ParsingException(std::string("The value of \"") + key +
+                               "\" should be a number");
+    }
+    return v.get<double>();
+}
+
+static double getReqDouble(const json &j, const char *key) {
+    return getDouble(j, key, false);
+}
+static double getOptDouble(const json &j, const char *key) {
+    return getDouble(j, key, true);
+}
+
+// ---------------------------------------------------------------------------
+
+static json getObjectMember(const json &j, const char *key) {
+    if (!j.contains(key)) {
+        throw ParsingException(std::string("Missing \"") + key + "\" key");
+    }
+    const json obj = j[key];
+    if (!obj.is_object()) {
+        throw ParsingException(std::string("The value of \"") + key +
+                               "\" should be a object");
+    }
+    return obj;
+}
+
+// ---------------------------------------------------------------------------
+
+static json getArrayMember(const json &j, const char *key) {
+    if (!j.contains(key)) {
+        throw ParsingException(std::string("Missing \"") + key + "\" key");
+    }
+    const json obj = j[key];
+    if (!obj.is_array()) {
+        throw ParsingException(std::string("The value of \"") + key +
+                               "\" should be a array");
+    }
+    return obj;
+}
+
+// ---------------------------------------------------------------------------
+
+std::unique_ptr<MasterFile> MasterFile::parse(const std::string &text) {
+    std::unique_ptr<MasterFile> dmmf(new MasterFile());
+    json j;
+    try {
+        j = json::parse(text);
+    } catch (const std::exception &e) {
+        throw ParsingException(e.what());
+    }
+    if (!j.is_object()) {
+        throw ParsingException("Not an object");
+    }
+    dmmf->mFileType = getReqString(j, "file_type");
+    dmmf->mFormatVersion = getReqString(j, "format_version");
+    dmmf->mName = getOptString(j, "name");
+    dmmf->mVersion = getOptString(j, "version");
+    dmmf->mLicense = getOptString(j, "license");
+    dmmf->mDescription = getOptString(j, "description");
+    dmmf->mPublicationDate = getOptString(j, "publication_date");
+
+    if (j.contains("authority")) {
+        const json jAuthority = j["authority"];
+        if (!jAuthority.is_object()) {
+            throw ParsingException("authority is not a object");
+        }
+        dmmf->mAuthority.name = getOptString(jAuthority, "name");
+        dmmf->mAuthority.url = getOptString(jAuthority, "url");
+        dmmf->mAuthority.address = getOptString(jAuthority, "address");
+        dmmf->mAuthority.email = getOptString(jAuthority, "email");
+    }
+
+    if (j.contains("links")) {
+        const json jLinks = j["links"];
+        if (!jLinks.is_array()) {
+            throw ParsingException("links is not an array");
+        }
+        for (const json &jLink : jLinks) {
+            if (!jLink.is_object()) {
+                throw ParsingException("links[] item is not an object");
+            }
+            Link link;
+            link.href = getOptString(jLink, "href");
+            link.rel = getOptString(jLink, "rel");
+            link.type = getOptString(jLink, "type");
+            link.title = getOptString(jLink, "title");
+            dmmf->mLinks.emplace_back(std::move(link));
+        }
+    }
+    dmmf->mSourceCRS = getReqString(j, "source_crs");
+    dmmf->mTargetCRS = getReqString(j, "target_crs");
+    dmmf->mDefinitionCRS = getReqString(j, "definition_crs");
+    if (dmmf->mSourceCRS != dmmf->mDefinitionCRS) {
+        throw ParsingException(
+            "source_crs != definition_crs not currently supported");
+    }
+    dmmf->mReferenceEpoch = getOptString(j, "reference_epoch");
+    dmmf->mUncertaintyReferenceEpoch =
+        getOptString(j, "uncertainty_reference_epoch");
+    dmmf->mHorizontalOffsetUnit = getOptString(j, "horizontal_offset_unit");
+    if (!dmmf->mHorizontalOffsetUnit.empty() &&
+        dmmf->mHorizontalOffsetUnit != STR_METRE &&
+        dmmf->mHorizontalOffsetUnit != STR_DEGREE) {
+        throw ParsingException("Unsupported value for horizontal_offset_unit");
+    }
+    dmmf->mVerticalOffsetUnit = getOptString(j, "vertical_offset_unit");
+    if (!dmmf->mVerticalOffsetUnit.empty() &&
+        dmmf->mVerticalOffsetUnit != STR_METRE) {
+        throw ParsingException("Unsupported value for vertical_offset_unit");
+    }
+    dmmf->mHorizontalUncertaintyType =
+        getOptString(j, "horizontal_uncertainty_type");
+    dmmf->mHorizontalUncertaintyUnit =
+        getOptString(j, "horizontal_uncertainty_unit");
+    dmmf->mVerticalUncertaintyType =
+        getOptString(j, "vertical_uncertainty_type");
+    dmmf->mVerticalUncertaintyUnit =
+        getOptString(j, "vertical_uncertainty_unit");
+    dmmf->mHorizontalOffsetMethod = getOptString(j, "horizontal_offset_method");
+    if (!dmmf->mHorizontalOffsetMethod.empty() &&
+        dmmf->mHorizontalOffsetMethod != STR_ADDITION &&
+        dmmf->mHorizontalOffsetMethod != STR_GEOCENTRIC) {
+        throw ParsingException(
+            "Unsupported value for horizontal_offset_method");
+    }
+    dmmf->mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent"));
+
+    const json jTimeExtent = getObjectMember(j, "time_extent");
+    dmmf->mTimeExtent.first = Epoch(getReqString(jTimeExtent, "first"));
+    dmmf->mTimeExtent.last = Epoch(getReqString(jTimeExtent, "last"));
+
+    const json jComponents = getArrayMember(j, "components");
+    for (const json &jComponent : jComponents) {
+        dmmf->mComponents.emplace_back(Component::parse(jComponent));
+        const auto &comp = dmmf->mComponents.back();
+        if (comp.displacementType() == STR_HORIZONTAL ||
+            comp.displacementType() == STR_3D) {
+            if (dmmf->mHorizontalOffsetUnit.empty()) {
+                throw ParsingException("horizontal_offset_unit should be "
+                                       "defined as there is a component with "
+                                       "displacement_type = horizontal/3d");
+            }
+            if (dmmf->mHorizontalOffsetMethod.empty()) {
+                throw ParsingException("horizontal_offset_method should be "
+                                       "defined as there is a component with "
+                                       "displacement_type = horizontal/3d");
+            }
+        }
+        if (comp.displacementType() == STR_VERTICAL ||
+            comp.displacementType() == STR_3D) {
+            if (dmmf->mVerticalOffsetUnit.empty()) {
+                throw ParsingException("vertical_offset_unit should be defined "
+                                       "as there is a component with "
+                                       "displacement_type = vertical/3d");
+            }
+        }
+        if (dmmf->mHorizontalOffsetUnit == STR_DEGREE &&
+            comp.spatialModel().interpolationMethod != STR_BILINEAR) {
+            throw ParsingException("horizontal_offset_unit = degree can only "
+                                   "be used with interpolation_method = "
+                                   "bilinear");
+        }
+    }
+
+    if (dmmf->mHorizontalOffsetUnit == STR_DEGREE &&
+        dmmf->mHorizontalOffsetMethod != STR_ADDITION) {
+        throw ParsingException("horizontal_offset_unit = degree can only be "
+                               "used with horizontal_offset_method = addition");
+    }
+
+    return dmmf;
+}
+
+// ---------------------------------------------------------------------------
+
+SpatialExtent SpatialExtent::parse(const json &j) {
+    SpatialExtent ex;
+
+    const std::string type = getReqString(j, "type");
+    if (type != "bbox") {
+        throw ParsingException("unsupported type of extent");
+    }
+
+    const json jParameter = getObjectMember(j, "parameters");
+    const json jBbox = getArrayMember(jParameter, "bbox");
+    if (jBbox.size() != 4) {
+        throw ParsingException("bbox is not an array of 4 numeric elements");
+    }
+    for (int i = 0; i < 4; i++) {
+        if (!jBbox[i].is_number()) {
+            throw ParsingException(
+                "bbox is not an array of 4 numeric elements");
+        }
+    }
+    ex.mMinx = jBbox[0].get<double>();
+    ex.mMiny = jBbox[1].get<double>();
+    ex.mMaxx = jBbox[2].get<double>();
+    ex.mMaxy = jBbox[3].get<double>();
+
+    ex.mMinxRad = DegToRad(ex.mMinx);
+    ex.mMinyRad = DegToRad(ex.mMiny);
+    ex.mMaxxRad = DegToRad(ex.mMaxx);
+    ex.mMaxyRad = DegToRad(ex.mMaxy);
+
+    return ex;
+}
+
+// ---------------------------------------------------------------------------
+
+Component Component::parse(const json &j) {
+    Component comp;
+    if (!j.is_object()) {
+        throw ParsingException("component is not an object");
+    }
+    comp.mDescription = getOptString(j, "description");
+    comp.mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent"));
+    comp.mDisplacementType = getReqString(j, "displacement_type");
+    if (comp.mDisplacementType != STR_NONE &&
+        comp.mDisplacementType != STR_HORIZONTAL &&
+        comp.mDisplacementType != STR_VERTICAL &&
+        comp.mDisplacementType != STR_3D) {
+        throw ParsingException("Unsupported value for displacement_type");
+    }
+    comp.mUncertaintyType = getReqString(j, "uncertainty_type");
+    comp.mHorizontalUncertainty = getOptDouble(j, "horizontal_uncertainty");
+    comp.mVerticalUncertainty = getOptDouble(j, "vertical_uncertainty");
+
+    const json jSpatialModel = getObjectMember(j, "spatial_model");
+    comp.mSpatialModel.type = getReqString(jSpatialModel, "type");
+    comp.mSpatialModel.interpolationMethod =
+        getReqString(jSpatialModel, "interpolation_method");
+    if (comp.mSpatialModel.interpolationMethod != STR_BILINEAR &&
+        comp.mSpatialModel.interpolationMethod != STR_GEOCENTRIC_BILINEAR) {
+        throw ParsingException("Unsupported value for interpolation_method");
+    }
+    comp.mSpatialModel.filename = getReqString(jSpatialModel, "filename");
+    comp.mSpatialModel.md5Checksum =
+        getOptString(jSpatialModel, "md5_checksum");
+
+    const json jTimeFunction = getObjectMember(j, "time_function");
+    const std::string timeFunctionType = getReqString(jTimeFunction, "type");
+    const json jParameters = timeFunctionType == "constant"
+                                 ? json()
+                                 : getObjectMember(jTimeFunction, "parameters");
+
+    if (timeFunctionType == "constant") {
+        std::unique_ptr<ConstantTimeFunction> tf(new ConstantTimeFunction());
+        tf->type = timeFunctionType;
+        comp.mTimeFunction = std::move(tf);
+    } else if (timeFunctionType == "velocity") {
+        std::unique_ptr<VelocityTimeFunction> tf(new VelocityTimeFunction());
+        tf->type = timeFunctionType;
+        tf->referenceEpoch =
+            Epoch(getReqString(jParameters, "reference_epoch"));
+        comp.mTimeFunction = std::move(tf);
+    } else if (timeFunctionType == "step") {
+        std::unique_ptr<StepTimeFunction> tf(new StepTimeFunction());
+        tf->type = timeFunctionType;
+        tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch"));
+        comp.mTimeFunction = std::move(tf);
+    } else if (timeFunctionType == "reverse_step") {
+        std::unique_ptr<ReverseStepTimeFunction> tf(
+            new ReverseStepTimeFunction());
+        tf->type = timeFunctionType;
+        tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch"));
+        comp.mTimeFunction = std::move(tf);
+    } else if (timeFunctionType == "piecewise") {
+        std::unique_ptr<PiecewiseTimeFunction> tf(new PiecewiseTimeFunction());
+        tf->type = timeFunctionType;
+        tf->beforeFirst = getReqString(jParameters, "before_first");
+        if (tf->beforeFirst != "zero" && tf->beforeFirst != "constant" &&
+            tf->beforeFirst != "linear") {
+            throw ParsingException("Unsupported value for before_first");
+        }
+        tf->afterLast = getReqString(jParameters, "after_last");
+        if (tf->afterLast != "zero" && tf->afterLast != "constant" &&
+            tf->afterLast != "linear") {
+            throw ParsingException("Unsupported value for afterLast");
+        }
+        const json jModel = getArrayMember(jParameters, "model");
+        for (const json &jModelElt : jModel) {
+            if (!jModelElt.is_object()) {
+                throw ParsingException("model[] element is not an object");
+            }
+            PiecewiseTimeFunction::EpochScaleFactorTuple tuple;
+            tuple.epoch = Epoch(getReqString(jModelElt, "epoch"));
+            tuple.scaleFactor = getReqDouble(jModelElt, "scale_factor");
+            tf->model.emplace_back(std::move(tuple));
+        }
+
+        comp.mTimeFunction = std::move(tf);
+    } else if (timeFunctionType == "exponential") {
+        std::unique_ptr<ExponentialTimeFunction> tf(
+            new ExponentialTimeFunction());
+        tf->type = timeFunctionType;
+        tf->referenceEpoch =
+            Epoch(getReqString(jParameters, "reference_epoch"));
+        tf->endEpoch = Epoch(getOptString(jParameters, "end_epoch"));
+        tf->relaxationConstant =
+            getReqDouble(jParameters, "relaxation_constant");
+        if (tf->relaxationConstant <= 0.0) {
+            throw ParsingException("Invalid value for relaxation_constant");
+        }
+        tf->beforeScaleFactor =
+            getReqDouble(jParameters, "before_scale_factor");
+        tf->initialScaleFactor =
+            getReqDouble(jParameters, "initial_scale_factor");
+        tf->finalScaleFactor = getReqDouble(jParameters, "final_scale_factor");
+        comp.mTimeFunction = std::move(tf);
+    } else {
+        throw ParsingException("Unsupported type of time function: " +
+                               timeFunctionType);
+    }
+
+    return comp;
+}
+
+// ---------------------------------------------------------------------------
+
+double Component::ConstantTimeFunction::evaluateAt(double) const { return 1.0; }
+
+// ---------------------------------------------------------------------------
+
+double Component::VelocityTimeFunction::evaluateAt(double dt) const {
+    return dt - referenceEpoch.toDecimalYear();
+}
+
+// ---------------------------------------------------------------------------
+
+double Component::StepTimeFunction::evaluateAt(double dt) const {
+    if (dt < stepEpoch.toDecimalYear())
+        return 0.0;
+    return 1.0;
+}
+
+// ---------------------------------------------------------------------------
+
+double Component::ReverseStepTimeFunction::evaluateAt(double dt) const {
+    if (dt < stepEpoch.toDecimalYear())
+        return -1.0;
+    return 0.0;
+}
+
+// ---------------------------------------------------------------------------
+
+double Component::PiecewiseTimeFunction::evaluateAt(double dt) const {
+    if (model.empty()) {
+        return 0.0;
+    }
+
+    const double dt1 = model[0].epoch.toDecimalYear();
+    if (dt < dt1) {
+        if (beforeFirst == "zero")
+            return 0.0;
+        if (beforeFirst == "constant" || model.size() == 1)
+            return model[0].scaleFactor;
+
+        // linear
+        const double f1 = model[0].scaleFactor;
+        const double dt2 = model[1].epoch.toDecimalYear();
+        const double f2 = model[1].scaleFactor;
+        if (dt1 == dt2)
+            return f1;
+        return (f1 * (dt2 - dt) + f2 * (dt - dt1)) / (dt2 - dt1);
+    }
+    for (size_t i = 1; i < model.size(); i++) {
+        const double dtip1 = model[i].epoch.toDecimalYear();
+        if (dt < dtip1) {
+            const double dti = model[i - 1].epoch.toDecimalYear();
+            const double fip1 = model[i].scaleFactor;
+            const double fi = model[i - 1].scaleFactor;
+            return (fi * (dtip1 - dt) + fip1 * (dt - dti)) / (dtip1 - dti);
+        }
+    }
+    if (afterLast == "zero") {
+        return 0.0;
+    }
+    if (afterLast == "constant" || model.size() == 1)
+        return model.back().scaleFactor;
+
+    // linear
+    const double dtnm1 = model[model.size() - 2].epoch.toDecimalYear();
+    const double fnm1 = model[model.size() - 2].scaleFactor;
+    const double dtn = model.back().epoch.toDecimalYear();
+    const double fn = model.back().scaleFactor;
+    if (dtnm1 == dtn)
+        return fn;
+    return (fnm1 * (dtn - dt) + fn * (dt - dtnm1)) / (dtn - dtnm1);
+}
+
+// ---------------------------------------------------------------------------
+
+double Component::ExponentialTimeFunction::evaluateAt(double dt) const {
+    const double t0 = referenceEpoch.toDecimalYear();
+    if (dt < t0)
+        return beforeScaleFactor;
+    if (!endEpoch.toString().empty()) {
+        dt = std::min(dt, endEpoch.toDecimalYear());
+    }
+    return initialScaleFactor +
+           (finalScaleFactor - initialScaleFactor) *
+               (1.0 - std::exp(-(dt - t0) / relaxationConstant));
+}
+
+// ---------------------------------------------------------------------------
+
+inline void DeltaEastingNorthingToLongLat(double cosphi, double de, double dn,
+                                          double a, double b, double es,
+                                          double &dlam, double &dphi) {
+    const double oneMinuX = es * (1 - cosphi * cosphi);
+    const double X = 1 - oneMinuX;
+    const double sqrtX = sqrt(X);
+#if 0
+    // With es of Earth, absolute/relative error is at most 2e-8
+    // const double sqrtX = 1 - 0.5 * oneMinuX * (1 + 0.25 * oneMinuX);
+#endif
+    dlam = de * sqrtX / (a * cosphi);
+    dphi = dn * a * sqrtX * X / (b * b);
+}
+
+// ---------------------------------------------------------------------------
+
+template <class Grid, class GridSet, class EvaluatorIface>
+Evaluator<Grid, GridSet, EvaluatorIface>::Evaluator(
+    std::unique_ptr<MasterFile> &&model, EvaluatorIface &iface, double a,
+    double b)
+    : mModel(std::move(model)), mA(a), mB(b), mEs(1 - (b * b) / (a * a)),
+      mIsHorizontalUnitDegree(mModel->horizontalOffsetUnit() == STR_DEGREE),
+      mIsAddition(mModel->horizontalOffsetMethod() == STR_ADDITION),
+      mIsGeographicCRS(iface.isGeographicCRS(mModel->definitionCRS())) {
+    if (!mIsGeographicCRS && mIsHorizontalUnitDegree) {
+        throw EvaluatorException(
+            "definition_crs = projected CRS and "
+            "horizontal_offset_unit = degree are incompatible");
+    }
+    if (!mIsGeographicCRS && !mIsAddition) {
+        throw EvaluatorException(
+            "definition_crs = projected CRS and "
+            "horizontal_offset_method = geocentric are incompatible");
+    }
+    mComponents.reserve(mModel->components().size());
+    for (const auto &comp : mModel->components()) {
+        mComponents.emplace_back(std::unique_ptr<ComponentEx<Grid, GridSet>>(
+            new ComponentEx<Grid, GridSet>(comp)));
+        if (!mIsGeographicCRS && !mComponents.back()->isBilinearInterpolation) {
+            throw EvaluatorException(
+                "definition_crs = projected CRS and "
+                "interpolation_method = geocentric_bilinear are incompatible");
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+template <class Grid, class GridSet, class EvaluatorIface>
+void Evaluator<Grid, GridSet, EvaluatorIface>::clearGridCache() {
+    for (auto &comp : mComponents) {
+        comp->clearGridCache();
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+#ifdef DEBUG_DEFMODEL
+
+static std::string shortName(const Component &comp) {
+    const auto &desc = comp.description();
+    return desc.substr(0, desc.find('\n')) + " (" +
+           comp.spatialModel().filename + ")";
+}
+
+static std::string toString(double val) {
+    char buffer[32];
+    snprintf(buffer, sizeof(buffer), "%.9g", val);
+    return buffer;
+}
+
+#endif
+
+// ---------------------------------------------------------------------------
+
+static bool bboxCheck(double &x, double &y, bool forInverseComputation,
+                      const double minx, const double miny, const double maxx,
+                      const double maxy, const double EPS,
+                      const double extraMarginForInverse) {
+    if (x < minx - EPS || x > maxx + EPS || y < miny - EPS || y > maxy + EPS) {
+        if (!forInverseComputation) {
+            return false;
+        }
+        // In case of iterative computation for inverse, allow to be a
+        // slightly bit outside of the grid and clamp to the edges
+        bool xOk = false;
+        if (x >= minx - EPS && x <= maxx + EPS) {
+            xOk = true;
+        } else if (x > minx - extraMarginForInverse && x < minx) {
+            x = minx;
+            xOk = true;
+        } else if (x < maxx + extraMarginForInverse && x > maxx) {
+            x = maxx;
+            xOk = true;
+        }
+
+        bool yOk = false;
+        if (y >= miny - EPS && y <= maxy + EPS) {
+            yOk = true;
+        } else if (y > miny - extraMarginForInverse && y < miny) {
+            y = miny;
+            yOk = true;
+        } else if (y < maxy + extraMarginForInverse && y > maxy) {
+            y = maxy;
+            yOk = true;
+        }
+
+        return xOk && yOk;
+    }
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+template <class Grid, class GridSet, class EvaluatorIface>
+bool Evaluator<Grid, GridSet, EvaluatorIface>::forward(
+    EvaluatorIface &iface, double x, double y, double z, double t,
+    bool forInverseComputation, double &x_out, double &y_out, double &z_out)
+
+{
+    x_out = x;
+    y_out = y;
+    z_out = z;
+
+    const double EPS = mIsGeographicCRS ? 1e-10 : 1e-5;
+
+    // Check against global model spatial extent, potentially wrapping
+    // longitude to match
+    {
+        const auto &extent = mModel->extent();
+        const double minx = extent.minxNormalized(mIsGeographicCRS);
+        const double maxx = extent.maxxNormalized(mIsGeographicCRS);
+        if (mIsGeographicCRS) {
+            while (x < minx - EPS) {
+                x += 2.0 * DEFMODEL_PI;
+            }
+            while (x > maxx + EPS) {
+                x -= 2.0 * DEFMODEL_PI;
+            }
+        }
+        const double miny = extent.minyNormalized(mIsGeographicCRS);
+        const double maxy = extent.maxyNormalized(mIsGeographicCRS);
+        const double extraMarginForInverse =
+            mIsGeographicCRS ? DegToRad(0.1) : 10000;
+        if (!bboxCheck(x, y, forInverseComputation, minx, miny, maxx, maxy, EPS,
+                       extraMarginForInverse)) {
+#ifdef DEBUG_DEFMODEL
+            iface.log("Calculation point " + toString(x) + "," + toString(y) +
+                      " is outside the extents of the deformation model");
+#endif
+            return false;
+        }
+    }
+
+    // Check against global model temporal extent
+    {
+        const auto &timeExtent = mModel->timeExtent();
+        if (t < timeExtent.first.toDecimalYear() ||
+            t > timeExtent.last.toDecimalYear()) {
+#ifdef DEBUG_DEFMODEL
+            iface.log("Calculation epoch " + toString(t) +
+                      " is not valid for the deformation model");
+#endif
+            return false;
+        }
+    }
+
+    // For mIsHorizontalUnitDegree
+    double dlam = 0;
+    double dphi = 0;
+
+    // For !mIsHorizontalUnitDegree
+    double de = 0;
+    double dn = 0;
+
+    double dz = 0;
+
+    bool sincosphiInitialized = false;
+    double sinphi = 0;
+    double cosphi = 0;
+
+    for (auto &compEx : mComponents) {
+        const auto &comp = compEx->component;
+        if (compEx->displacementType == DisplacementType::NONE) {
+            continue;
+        }
+        const auto &extent = comp.extent();
+        double xForGrid = x;
+        double yForGrid = y;
+        const double minx = extent.minxNormalized(mIsGeographicCRS);
+        const double maxx = extent.maxxNormalized(mIsGeographicCRS);
+        const double miny = extent.minyNormalized(mIsGeographicCRS);
+        const double maxy = extent.maxyNormalized(mIsGeographicCRS);
+        const double extraMarginForInverse = 0;
+        if (!bboxCheck(xForGrid, yForGrid, forInverseComputation, minx, miny,
+                       maxx, maxy, EPS, extraMarginForInverse)) {
+#ifdef DEBUG_DEFMODEL
+            iface.log(
+                "Skipping component " + shortName(comp) +
+                " due to point being outside of its declared spatial extent.");
+#endif
+            continue;
+        }
+        xForGrid = std::max(xForGrid, minx);
+        yForGrid = std::max(yForGrid, miny);
+        xForGrid = std::min(xForGrid, maxx);
+        yForGrid = std::min(yForGrid, maxy);
+        const auto tfactor = compEx->evaluateAt(t);
+        if (tfactor == 0.0) {
+#ifdef DEBUG_DEFMODEL
+            iface.log("Skipping component " + shortName(comp) +
+                      " due to time function evaluating to 0.");
+#endif
+            continue;
+        }
+
+#ifdef DEBUG_DEFMODEL
+        iface.log("Entering component " + shortName(comp) +
+                  " with time function evaluating to " + toString(tfactor) +
+                  ".");
+#endif
+
+        if (compEx->gridSet == nullptr) {
+            compEx->gridSet = iface.open(comp.spatialModel().filename);
+            if (compEx->gridSet == nullptr) {
+                return false;
+            }
+        }
+        const Grid *grid = compEx->gridSet->gridAt(xForGrid, yForGrid);
+        if (grid == nullptr) {
+#ifdef DEBUG_DEFMODEL
+            iface.log("Skipping component " + shortName(comp) +
+                      " due to no grid found for this point in the grid set.");
+#endif
+            continue;
+        }
+        if (grid->width < 2 || grid->height < 2) {
+            return false;
+        }
+        const double ix_d = (xForGrid - grid->minx) / grid->resx;
+        const double iy_d = (yForGrid - grid->miny) / grid->resy;
+        if (ix_d < -EPS || iy_d < -EPS || ix_d + 1 >= grid->width + EPS ||
+            iy_d + 1 >= grid->height + EPS) {
+#ifdef DEBUG_DEFMODEL
+            iface.log("Skipping component " + shortName(comp) +
+                      " due to point being outside of actual spatial extent of "
+                      "grid " +
+                      grid->name() + ".");
+#endif
+            continue;
+        }
+        const int ix0 = std::min(static_cast<int>(ix_d), grid->width - 2);
+        const int iy0 = std::min(static_cast<int>(iy_d), grid->height - 2);
+        const int ix1 = ix0 + 1;
+        const int iy1 = iy0 + 1;
+        const double frct_x = ix_d - ix0;
+        const double frct_y = iy_d - iy0;
+        const double one_minus_frct_x = 1. - frct_x;
+        const double one_minus_frct_y = 1. - frct_y;
+        const double m00 = one_minus_frct_x * one_minus_frct_y;
+        const double m10 = frct_x * one_minus_frct_y;
+        const double m01 = one_minus_frct_x * frct_y;
+        const double m11 = frct_x * frct_y;
+
+        if (compEx->displacementType == DisplacementType::VERTICAL) {
+            double dz00 = 0;
+            double dz01 = 0;
+            double dz10 = 0;
+            double dz11 = 0;
+            if (!grid->getZOffset(ix0, iy0, dz00) ||
+                !grid->getZOffset(ix1, iy0, dz10) ||
+                !grid->getZOffset(ix0, iy1, dz01) ||
+                !grid->getZOffset(ix1, iy1, dz11)) {
+                return false;
+            }
+            const double dzInterp =
+                dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11;
+#ifdef DEBUG_DEFMODEL
+            iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " +
+                      toString(dzInterp) + ".");
+#endif
+            dz += tfactor * dzInterp;
+        } else if (mIsHorizontalUnitDegree) {
+            double dx00 = 0;
+            double dy00 = 0;
+            double dx01 = 0;
+            double dy01 = 0;
+            double dx10 = 0;
+            double dy10 = 0;
+            double dx11 = 0;
+            double dy11 = 0;
+            if (compEx->displacementType == DisplacementType::HORIZONTAL) {
+                if (!grid->getLonLatOffset(ix0, iy0, dx00, dy00) ||
+                    !grid->getLonLatOffset(ix1, iy0, dx10, dy10) ||
+                    !grid->getLonLatOffset(ix0, iy1, dx01, dy01) ||
+                    !grid->getLonLatOffset(ix1, iy1, dx11, dy11)) {
+                    return false;
+                }
+            } else /* if (compEx->displacementType == DisplacementType::THREE_D) */ {
+                double dz00 = 0;
+                double dz01 = 0;
+                double dz10 = 0;
+                double dz11 = 0;
+                if (!grid->getLonLatZOffset(ix0, iy0, dx00, dy00, dz00) ||
+                    !grid->getLonLatZOffset(ix1, iy0, dx10, dy10, dz10) ||
+                    !grid->getLonLatZOffset(ix0, iy1, dx01, dy01, dz01) ||
+                    !grid->getLonLatZOffset(ix1, iy1, dx11, dy11, dz11)) {
+                    return false;
+                }
+                const double dzInterp =
+                    dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11;
+#ifdef DEBUG_DEFMODEL
+                iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " +
+                          toString(dzInterp) + ".");
+#endif
+                dz += tfactor * dzInterp;
+            }
+            const double dlamInterp =
+                dx00 * m00 + dx01 * m01 + dx10 * m10 + dx11 * m11;
+            const double dphiInterp =
+                dy00 * m00 + dy01 * m01 + dy10 * m10 + dy11 * m11;
+#ifdef DEBUG_DEFMODEL
+            iface.log("tfactor * dlamInterp = " + toString(tfactor) + " * " +
+                      toString(dlamInterp) + ".");
+            iface.log("tfactor * dphiInterp = " + toString(tfactor) + " * " +
+                      toString(dphiInterp) + ".");
+#endif
+            dlam += tfactor * dlamInterp;
+            dphi += tfactor * dphiInterp;
+        } else /* horizontal unit is metre */ {
+            double de00 = 0;
+            double dn00 = 0;
+            double de01 = 0;
+            double dn01 = 0;
+            double de10 = 0;
+            double dn10 = 0;
+            double de11 = 0;
+            double dn11 = 0;
+            if (compEx->displacementType == DisplacementType::HORIZONTAL) {
+                if (!grid->getEastingNorthingOffset(ix0, iy0, de00, dn00) ||
+                    !grid->getEastingNorthingOffset(ix1, iy0, de10, dn10) ||
+                    !grid->getEastingNorthingOffset(ix0, iy1, de01, dn01) ||
+                    !grid->getEastingNorthingOffset(ix1, iy1, de11, dn11)) {
+                    return false;
+                }
+            } else /* if (compEx->displacementType == DisplacementType::THREE_D) */ {
+                double dz00 = 0;
+                double dz01 = 0;
+                double dz10 = 0;
+                double dz11 = 0;
+                if (!grid->getEastingNorthingZOffset(ix0, iy0, de00, dn00,
+                                                     dz00) ||
+                    !grid->getEastingNorthingZOffset(ix1, iy0, de10, dn10,
+                                                     dz10) ||
+                    !grid->getEastingNorthingZOffset(ix0, iy1, de01, dn01,
+                                                     dz01) ||
+                    !grid->getEastingNorthingZOffset(ix1, iy1, de11, dn11,
+                                                     dz11)) {
+                    return false;
+                }
+                const double dzInterp =
+                    dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11;
+#ifdef DEBUG_DEFMODEL
+                iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " +
+                          toString(dzInterp) + ".");
+#endif
+                dz += tfactor * dzInterp;
+            }
+            if (compEx->isBilinearInterpolation) {
+                const double deInterp =
+                    de00 * m00 + de01 * m01 + de10 * m10 + de11 * m11;
+                const double dnInterp =
+                    dn00 * m00 + dn01 * m01 + dn10 * m10 + dn11 * m11;
+#ifdef DEBUG_DEFMODEL
+                iface.log("tfactor * deInterp = " + toString(tfactor) + " * " +
+                          toString(deInterp) + ".");
+                iface.log("tfactor * dnInterp = " + toString(tfactor) + " * " +
+                          toString(dnInterp) + ".");
+#endif
+                de += tfactor * deInterp;
+                dn += tfactor * dnInterp;
+            } else /* geocentric_bilinear */ {
+                double dX;
+                double dY;
+                double dZ;
+
+                auto iter = compEx->mapGrids.find(grid);
+                if (iter == compEx->mapGrids.end()) {
+                    GridEx<Grid> gridWithCache(grid);
+                    iter = compEx->mapGrids
+                               .insert(std::pair<const Grid *, GridEx<Grid>>(
+                                   grid, std::move(gridWithCache)))
+                               .first;
+                }
+                GridEx<Grid> &gridwithCacheRef = iter->second;
+
+                gridwithCacheRef.getBilinearGeocentric(
+                    ix0, iy0, de00, dn00, de01, dn01, de10, dn10, de11, dn11,
+                    m00, m01, m10, m11, dX, dY, dZ);
+                if (!sincosphiInitialized) {
+                    sincosphiInitialized = true;
+                    sinphi = sin(y);
+                    cosphi = cos(y);
+                }
+                const double lam_rel_to_cell_center =
+                    (frct_x - 0.5) * grid->resx;
+                // Use small-angle approximation of sin/cos when reasonable
+                // Max abs/rel error on cos is 3.9e-9 and on sin 1.3e-11
+                const double sinlam =
+                    gridwithCacheRef.smallResx
+                        ? lam_rel_to_cell_center *
+                              (1 -
+                               (1. / 6) * (lam_rel_to_cell_center *
+                                           lam_rel_to_cell_center))
+                        : sin(lam_rel_to_cell_center);
+                const double coslam = gridwithCacheRef.smallResx
+                                          ? (1 -
+                                             0.5 * (lam_rel_to_cell_center *
+                                                    lam_rel_to_cell_center))
+                                          : cos(lam_rel_to_cell_center);
+
+                // Convert back from geocentric deltas to easting, northing
+                // deltas
+                const double deInterp = -dX * sinlam + dY * coslam;
+                const double dnInterp =
+                    (-dX * coslam - dY * sinlam) * sinphi + dZ * cosphi;
+#ifdef DEBUG_DEFMODEL
+                iface.log("After geocentric_bilinear interpolation: tfactor * "
+                          "deInterp = " +
+                          toString(tfactor) + " * " + toString(deInterp) + ".");
+                iface.log("After geocentric_bilinear interpolation: tfactor * "
+                          "dnInterp = " +
+                          toString(tfactor) + " * " + toString(dnInterp) + ".");
+#endif
+                de += tfactor * deInterp;
+                dn += tfactor * dnInterp;
+            }
+        }
+    }
+
+    // Apply shifts depending on horizontal_offset_unit and
+    // horizontal_offset_method
+    if (mIsHorizontalUnitDegree) {
+        x_out += dlam;
+        y_out += dphi;
+    } else {
+#ifdef DEBUG_DEFMODEL
+        iface.log("Total sum of de: " + toString(de));
+        iface.log("Total sum of dn: " + toString(dn));
+#endif
+        if (mIsAddition && !mIsGeographicCRS) {
+            x_out += de;
+            y_out += dn;
+        } else if (mIsAddition) {
+            // Simple way of adding the offset
+            if (!sincosphiInitialized) {
+                cosphi = cos(y);
+            }
+            DeltaEastingNorthingToLongLat(cosphi, de, dn, mA, mB, mEs, dlam,
+                                          dphi);
+#ifdef DEBUG_DEFMODEL
+            iface.log("Result dlam: " + toString(dlam));
+            iface.log("Result dphi: " + toString(dphi));
+#endif
+            x_out += dlam;
+            y_out += dphi;
+        } else {
+            // Geocentric way of adding the offset
+            if (!sincosphiInitialized) {
+                sinphi = sin(y);
+                cosphi = cos(y);
+            }
+            const double sinlam = sin(x);
+            const double coslam = cos(x);
+            const double dnsinphi = dn * sinphi;
+            const double dX = -de * sinlam - dnsinphi * coslam;
+            const double dY = de * coslam - dnsinphi * sinlam;
+            const double dZ = dn * cosphi;
+
+            double X;
+            double Y;
+            double Z;
+            iface.geographicToGeocentric(x, y, 0, mA, mB, mEs, X, Y, Z);
+#ifdef DEBUG_DEFMODEL
+            iface.log("Geocentric coordinate before: " + toString(X) + "," +
+                      toString(Y) + "," + toString(Z));
+            iface.log("Geocentric shift: " + toString(dX) + "," + toString(dY) +
+                      "," + toString(dZ));
+#endif
+            X += dX;
+            Y += dY;
+            Z += dZ;
+#ifdef DEBUG_DEFMODEL
+            iface.log("Geocentric coordinate after: " + toString(X) + "," +
+                      toString(Y) + "," + toString(Z));
+#endif
+
+            double h_out_ignored;
+            iface.geocentricToGeographic(X, Y, Z, mA, mB, mEs, x_out, y_out,
+                                         h_out_ignored);
+        }
+    }
+#ifdef DEBUG_DEFMODEL
+    iface.log("Total sum of dz: " + toString(dz));
+#endif
+    z_out += dz;
+
+    return true;
+}
+
+// ---------------------------------------------------------------------------
+
+template <class Grid, class GridSet, class EvaluatorIface>
+bool Evaluator<Grid, GridSet, EvaluatorIface>::inverse(
+    EvaluatorIface &iface, double x, double y, double z, double t,
+    double &x_out, double &y_out, double &z_out)
+
+{
+    x_out = x;
+    y_out = y;
+    z_out = z;
+    constexpr double EPS_HORIZ = 1e-12;
+    constexpr double EPS_VERT = 1e-3;
+    constexpr bool forInverseComputation = true;
+    for (int i = 0; i < 10; i++) {
+#ifdef DEBUG_DEFMODEL
+        iface.log("Iteration " + std::to_string(i) + ": before foward: x=" +
+                  toString(x_out) + ", y=" + toString(y_out));
+#endif
+        double x_new;
+        double y_new;
+        double z_new;
+        if (!forward(iface, x_out, y_out, z_out, t, forInverseComputation,
+                     x_new, y_new, z_new)) {
+            return false;
+        }
+#ifdef DEBUG_DEFMODEL
+        iface.log("After foward: x=" + toString(x_new) + ", y=" +
+                  toString(y_new));
+#endif
+        const double dx = x_new - x;
+        const double dy = y_new - y;
+        const double dz = z_new - z;
+        x_out -= dx;
+        y_out -= dy;
+        z_out -= dz;
+        if (std::max(std::fabs(dx), std::fabs(dy)) < EPS_HORIZ &&
+            std::fabs(dz) < EPS_VERT) {
+            return true;
+        }
+    }
+    return false;
+}
+
+// ---------------------------------------------------------------------------
+
+} // namespace DEFORMATON_MODEL_NAMESPACE
diff --git a/src/transformations/deformation.cpp b/src/transformations/deformation.cpp
index 8aee50c9..8ce02bee 100644
--- a/src/transformations/deformation.cpp
+++ b/src/transformations/deformation.cpp
@@ -124,7 +124,7 @@ static bool get_grid_values(PJ* P,
     }
 
     bool must_retry = false;
-    if( !pj_bilinear_interpolation_three_samples(grid, lp,
+    if( !pj_bilinear_interpolation_three_samples(P->ctx, grid, lp,
                                                  sampleE, sampleN, sampleU,
                                                  vx, vy, vz,
                                                  must_retry) )
diff --git a/src/transformations/xyzgridshift.cpp b/src/transformations/xyzgridshift.cpp
index e1c76518..e37e874d 100644
--- a/src/transformations/xyzgridshift.cpp
+++ b/src/transformations/xyzgridshift.cpp
@@ -105,7 +105,7 @@ static bool get_grid_values(PJ* P,
     }
 
     bool must_retry = false;
-    if( !pj_bilinear_interpolation_three_samples(grid, lp,
+    if( !pj_bilinear_interpolation_three_samples(P->ctx, grid, lp,
                                                  sampleX, sampleY, sampleZ,
                                                  dx, dy, dz,
                                                  must_retry) )