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

Fixes #1001

Co-authored-by: Chris Crook <ccrook@linz.govt.nz>

7 files changed, 1674 insertions, 11 deletions

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index dcb4c26d..9fc29103 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -16,6 +16,7 @@ proj_add_gie_test("adams_ws1" "gie/adams_ws1.gie")
 proj_add_gie_test("adams_ws2" "gie/adams_ws2.gie")
 proj_add_gie_test("guyou" "gie/guyou.gie")
 proj_add_gie_test("peirce_q" "gie/peirce_q.gie")
+proj_add_gie_test("defmodel" "gie/defmodel.gie")
 
 # GIGS tests. Uncommented tests are expected to fail due to issues with
 # various projections. Should be investigated further and fixed.
diff --git a/test/gie/Makefile.am b/test/gie/Makefile.am
index 1637a38e..e6bb9582 100644
--- a/test/gie/Makefile.am
+++ b/test/gie/Makefile.am
@@ -15,8 +15,8 @@ EXTRA_DIST = 4D-API_cs2cs-style.gie \
 	adams_ws1.gie \
 	adams_ws2.gie \
 	guyou.gie \
-	peirce_q.gie
-
+	peirce_q.gie \
+	defmodel.gie
 
 PROJ_LIB ?= ../../data/for_tests
 
@@ -65,4 +65,7 @@ guyou: guyou.gie
 peirce_q: peirce_q.gie
 	PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $<
 
-check-local: 4D-API-cs2cs-style GDA axisswap builtins deformation ellipsoid more_builtins unitconvert DHDN_ETRS89 geotiff_grids adams_hemi adams_ws1 adams_ws2 guyou peirce_q
+defmodel: defmodel.gie
+	PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) $(GIEEXE) $<
+
+check-local: 4D-API-cs2cs-style GDA axisswap builtins deformation ellipsoid more_builtins unitconvert DHDN_ETRS89 geotiff_grids adams_hemi adams_ws1 adams_ws2 guyou peirce_q defmodel
diff --git a/test/gie/defmodel.gie b/test/gie/defmodel.gie
new file mode 100644
index 00000000..94dd5784
--- /dev/null
+++ b/test/gie/defmodel.gie
@@ -0,0 +1,124 @@
+
+-------------------------------------------------------------------------------
+===============================================================================
+Test +proj=defmodel
+===============================================================================
+
+<gie-strict>
+
+# Missing +model
+operation   +proj=defmodel
+expect failure errno no_args
+
+# +model doesn't point to an existing file
+operation   +proj=defmodel +model=i_do_not_exist
+expect failure errno invalid_arg
+
+# Not a JSON file
+operation   +proj=defmodel +model=proj.ini
+expect failure errno invalid_arg
+
+# Horizontal deformation with horizontal unit = degree
+operation   +proj=defmodel +model=tests/simple_model_degree_horizontal.json
+tolerance 0.1 mm
+accept 2 49 30 2020
+expect 3 51 30 2020
+roundtrip 1
+
+# 3D deformation with horizontal unit = degree
+operation   +proj=defmodel +model=tests/simple_model_degree_3d.json
+tolerance 0.1 mm
+accept 2 49 30 2020
+expect 3 51 33 2020
+roundtrip 1
+
+# Horizontal deformation with horizontal unit = metre
+operation   +proj=pipeline +step +inv +proj=merc +step +proj=defmodel +model=tests/simple_model_metre_horizontal.json +step +proj=merc
+tolerance 0.1 mm
+accept 10 20 30 2020
+expect 11 22 30 2020
+roundtrip 1
+
+# 3D deformation with horizontal unit = metre
+operation   +proj=pipeline +step +inv +proj=merc +step +proj=defmodel +model=tests/simple_model_metre_3d.json +step +proj=merc
+tolerance 0.1 mm
+accept 10 20 30 2020
+expect 11 22 33 2020
+roundtrip 1
+
+# 3D deformation with horizontal unit = metre and a projeced grid
+operation   +proj=pipeline +step +proj=defmodel +model=tests/simple_model_projected.json
+tolerance 0.1 mm
+
+accept 1500200.0 5400400.0 30 2020
+expect 1500200.588 5400399.722 30.6084 2020
+roundtrip 1
+
+# South-west corner
+accept 1500000.0 5400000.0 30 2020
+expect 1500000.4 5399999.8 30.84 2020
+roundtrip 1
+
+# South-east corner
+accept 1501000.0 5400000.0 30 2020
+expect 1501000.5 5399999.75 30.75 2020
+roundtrip 1
+
+# North-west corner
+accept 1500000.0 5401000.0 30 2020
+expect 1500000.8 5400999.6 30.36 2020
+roundtrip 1
+
+# North-east corner
+accept 1501000.0 5401000.0 30 2020
+expect 1501001.0 5400999.7 30 2020
+roundtrip 1
+
+# Test geocentric addition of components
+operation   +proj=pipeline +step +inv +proj=merc +step +proj=defmodel +model=tests/simple_model_metre_3d_geocentric.json +step +proj=merc
+tolerance 0.1 mm
+accept 10 20 30 2020
+expect 11 22 33 2020
+roundtrip 1
+
+# Vertical deformation with vertical unit = metre
+operation   +proj=defmodel +model=tests/simple_model_metre_vertical.json
+tolerance 0.1 mm
+accept 2 49 30 2020
+expect 2 49 33 2020
+roundtrip 1
+
+# Adjust for 360 degree longitude offsets
+operation   +proj=defmodel +model=tests/simple_model_metre_vertical.json
+tolerance 0.1 mm
+
+accept 362 49 30 2020
+expect 2 49 33 2020
+
+operation +proj=defmodel +model=tests/simple_model_wrap_east.json
+
+accept 165.9 -37.3 10 2020
+expect 165.9 -37.3 10.4525 2020
+
+operation +proj=defmodel +model=tests/simple_model_wrap_west.json
+
+accept 165.9 -37.3 10 2020
+expect 165.9 -37.3 10.4525 2020
+
+# Test geocentric bilinear interpolation method
+operation   +proj=defmodel +model=tests/simple_model_polar.json
+tolerance 0.1 mm
+
+accept  20 -90 15 2020
+expect 27.4743245365 -89.9999747721 18.0000 2020
+
+accept  120 -90 15 2020
+expect 27.4737934098 -89.9999747718 18.0000 2020
+
+accept  235 -89.5 15 2020
+expect -124.9986638571 -89.5000223708 17.3750 2020
+
+accept  45 -89.5 15 2020
+expect 44.9991295392 -89.4999759438 18.5469 2020
+
+</gie-strict>
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 62932aae..04c6acb6 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -170,3 +170,14 @@ target_link_libraries(test_network
 add_test(NAME test_network COMMAND test_network)
 set_property(TEST test_network
   PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT})
+
+
+add_executable(test_defmodel
+  main.cpp
+  test_defmodel.cpp)
+target_link_libraries(test_defmodel
+  GTest::gtest
+  ${PROJ_LIBRARIES})
+add_test(NAME test_defmodel COMMAND test_defmodel)
+set_property(TEST test_defmodel
+  PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT})
diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am
index 11a6473c..2dd0dbff 100644
--- a/test/unit/Makefile.am
+++ b/test/unit/Makefile.am
@@ -18,6 +18,7 @@ noinst_PROGRAMS += test_cpp_api
 noinst_PROGRAMS += gie_self_tests
 noinst_PROGRAMS += include_proj_h_from_c
 noinst_PROGRAMS += test_network
+noinst_PROGRAMS += test_defmodel
 
 pj_transform_test_SOURCES = pj_transform_test.cpp main.cpp
 pj_transform_test_LDADD = ../../src/libproj.la @GTEST_LIBS@
@@ -70,4 +71,12 @@ test_network_LDADD = ../../src/libproj.la @GTEST_LIBS@ @SQLITE3_LIBS@ @CURL_LIBS
 test_network-check: test_network
 	PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) PROJ_SOURCE_DATA=$(PROJ_LIB) ./test_network
 
-check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check proj_angular_io_test-check proj_context_test-check test_cpp_api-check gie_self_tests-check test_network-check
+test_defmodel_SOURCES = test_defmodel.cpp main.cpp
+test_defmodel_LDADD = ../../src/libproj.la @GTEST_LIBS@
+
+test_defmodel-check: test_defmodel
+	PROJ_LIB=$(PROJ_LIB) ./test_defmodel
+
+check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check \
+		proj_angular_io_test-check proj_context_test-check test_cpp_api-check \
+		gie_self_tests-check test_network-check test_defmodel-check
diff --git a/test/unit/test_defmodel.cpp b/test/unit/test_defmodel.cpp
new file mode 100644
index 00000000..b3228fd0
--- /dev/null
+++ b/test/unit/test_defmodel.cpp
@@ -0,0 +1,1515 @@
+/******************************************************************************
+ *
+ * Project:  PROJ
+ * Purpose:  Test deformation model
+ * Author:   Even Rouault <even dot rouault at spatialys dot com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2020, Even Rouault <even dot rouault at spatialys dot 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.
+ ****************************************************************************/
+
+#include "gtest_include.h"
+
+// to be able to use internal::toString
+#define FROM_PROJ_CPP
+#include "proj/internal/internal.hpp"
+
+#include "proj.h"
+#include "proj_internal.h"
+
+// Silence C4702 (unreachable code) due to some dummy implementation of the
+// interfaces of defmodel.hpp
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable : 4702)
+#endif
+
+#define PROJ_COMPILATION
+#define DEFORMATON_MODEL_NAMESPACE TestDeformationModel
+#include "transformations/defmodel.hpp"
+
+using namespace DEFORMATON_MODEL_NAMESPACE;
+
+namespace {
+
+constexpr double modelMinX = 158;
+constexpr double modelMinY = -58;
+constexpr double modelMaxX = 194;
+constexpr double modelMaxY = -25;
+
+static json getMinValidContent() {
+    json j;
+    j["file_type"] = "GeoTIFF";
+    j["format_version"] = "1.0";
+    j["source_crs"] = "EPSG:4959";
+    j["target_crs"] = "EPSG:7907";
+    j["definition_crs"] = "EPSG:4959";
+    j["extent"]["type"] = "bbox";
+    j["extent"]
+     ["parameters"] = {{"bbox", {modelMinX, modelMinY, modelMaxX, modelMaxY}}};
+    j["time_extent"]["first"] = "1900-01-01T00:00:00Z";
+    j["time_extent"]["last"] = "2050-01-01T00:00:00Z";
+    j["components"] = json::array();
+
+    return j;
+}
+
+// ---------------------------------------------------------------------------
+
+constexpr int IDX_CONSTANT = 0;
+constexpr int IDX_VELOCITY = 1;
+constexpr int IDX_STEP = 2;
+constexpr int IDX_REVERSE_STEP = 3;
+constexpr int IDX_PIECEWISE = 4;
+constexpr int IDX_EXPONENTIAL = 5;
+
+static json getFullValidContent() {
+    json j(getMinValidContent());
+    j["name"] = "name";
+    j["version"] = "version";
+    j["publication_date"] = "2018-07-01T00:00:00Z";
+    j["license"] = "license";
+    j["description"] = "description";
+    j["authority"]["name"] = "authority_name";
+    j["authority"]["url"] = "authority_url";
+    j["authority"]["address"] = "authority_address";
+    j["authority"]["email"] = "authority_email";
+    j["links"] = {{{"href", "href"},
+                   {"rel", "rel"},
+                   {"type", "type"},
+                   {"title", "title"}}};
+    j["reference_epoch"] = "2000-01-01T00:00:00Z";
+    j["uncertainty_reference_epoch"] = "2018-12-15T00:00:00Z";
+    j["horizontal_offset_method"] = "addition";
+    j["horizontal_offset_unit"] = "metre";
+    j["vertical_offset_unit"] = "metre";
+    j["horizontal_uncertainty_type"] = "circular 95% confidence limit";
+    j["horizontal_uncertainty_unit"] = "metre";
+    j["vertical_uncertainty_type"] = "95% confidence limit";
+    j["vertical_uncertainty_unit"] = "metre";
+    j["components"] = {{
+        {"description", "description"},
+        {"displacement_type", "horizontal"},
+        {"uncertainty_type", "none"},
+        {"horizontal_uncertainty", 0.01},
+        {"vertical_uncertainty", 0.02},
+        {"extent",
+         {{"type", "bbox"},
+          {"parameters",
+           {{"bbox", {modelMinX, modelMinY, modelMaxX, modelMaxY}}}}}},
+        {"spatial_model",
+         {
+             {"type", "GeoTIFF"},
+             {"interpolation_method", "bilinear"},
+             {"filename", "nzgd2000-ndm-grid02.tif"},
+             {"md5_checksum", "49fce8ab267be2c8d00d43683060a032"},
+         }},
+        {"time_function",
+         {
+             {"type", "constant"}, {"parameters", json::object()},
+         }},
+    }};
+
+    j["components"].push_back(j["components"][0]);
+    j["components"][IDX_VELOCITY]["time_function"] = {
+        {"type", "velocity"},
+        {"parameters", {{"reference_epoch", "2000-01-01T00:00:00Z"}}},
+    };
+
+    j["components"].push_back(j["components"][0]);
+    j["components"][IDX_STEP]["time_function"] = {
+        {"type", "step"},
+        {"parameters", {{"step_epoch", "2000-01-01T00:00:00Z"}}},
+    };
+
+    j["components"].push_back(j["components"][0]);
+    j["components"][IDX_REVERSE_STEP]["time_function"] = {
+        {"type", "reverse_step"},
+        {"parameters", {{"step_epoch", "2000-01-01T00:00:00Z"}}},
+    };
+
+    j["components"].push_back(j["components"][0]);
+    j["components"][IDX_PIECEWISE]["time_function"] = {
+        {"type", "piecewise"},
+        {"parameters",
+         {{"before_first", "zero"},
+          {"after_last", "constant"},
+          {"model",
+           {{{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}},
+            {{"epoch", "2017-01-01T00:00:00Z"}, {"scale_factor", 1.0}},
+            {{"epoch", "2017-01-01T00:00:00Z"}, {"scale_factor", 2.0}},
+            {{"epoch", "2018-01-01T00:00:00Z"}, {"scale_factor", 1.0}}}}}}};
+
+    j["components"].push_back(j["components"][0]);
+    j["components"][IDX_EXPONENTIAL]["time_function"] = {
+        {"type", "exponential"},
+        {"parameters",
+         {
+             {"reference_epoch", "2000-01-01T00:00:00Z"},
+             {"end_epoch", "2001-01-01T00:00:00Z"},
+             {"relaxation_constant", 2.0},
+             {"before_scale_factor", 0.0},
+             {"initial_scale_factor", 1.0},
+             {"final_scale_factor", 3.0},
+         }},
+    };
+
+    return j;
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, basic) {
+    EXPECT_THROW(MasterFile::parse("foo"), ParsingException);
+    EXPECT_THROW(MasterFile::parse("null"), ParsingException);
+    EXPECT_THROW(MasterFile::parse("{}"), ParsingException);
+
+    const auto jMinValid(getMinValidContent());
+    {
+        auto mf = MasterFile::parse(jMinValid.dump());
+        EXPECT_EQ(mf->fileType(), "GeoTIFF");
+        EXPECT_EQ(mf->formatVersion(), "1.0");
+        EXPECT_EQ(mf->sourceCRS(), "EPSG:4959");
+        EXPECT_EQ(mf->targetCRS(), "EPSG:7907");
+        EXPECT_EQ(mf->definitionCRS(), "EPSG:4959");
+        EXPECT_EQ(mf->extent().minx(), modelMinX);
+        EXPECT_EQ(mf->extent().miny(), modelMinY);
+        EXPECT_EQ(mf->extent().maxx(), modelMaxX);
+        EXPECT_EQ(mf->extent().maxy(), modelMaxY);
+        EXPECT_EQ(mf->timeExtent().first.toString(), "1900-01-01T00:00:00Z");
+        EXPECT_EQ(mf->timeExtent().last.toString(), "2050-01-01T00:00:00Z");
+    }
+
+    // Check that removing one of each required key causes an exception
+    for (const auto &kv : jMinValid.items()) {
+        json jcopy(jMinValid);
+        jcopy.erase(kv.key());
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["definition_crs"] = "EPSG:4326";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["file_type"] = 1;
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"].erase("type");
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"].erase("parameters");
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"]["parameters"].clear();
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"]["parameters"].erase("bbox");
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"]["parameters"]["bbox"] = "foo";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"]["parameters"]["bbox"] = {0, 1, 2};
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["extent"]["parameters"]["bbox"] = {0, 1, 2, "foo"};
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["time_extent"] = "foo";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["time_extent"].erase("first");
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jMinValid);
+        jcopy["time_extent"].erase("last");
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, full) {
+    const auto jFullValid(getFullValidContent());
+    auto mf = MasterFile::parse(jFullValid.dump());
+    EXPECT_EQ(mf->name(), "name");
+    EXPECT_EQ(mf->version(), "version");
+    EXPECT_EQ(mf->publicationDate(), "2018-07-01T00:00:00Z");
+    EXPECT_EQ(mf->license(), "license");
+    EXPECT_EQ(mf->description(), "description");
+    EXPECT_EQ(mf->authority().name, "authority_name");
+    EXPECT_EQ(mf->authority().url, "authority_url");
+    EXPECT_EQ(mf->authority().address, "authority_address");
+    EXPECT_EQ(mf->authority().email, "authority_email");
+    EXPECT_EQ(mf->links().size(), 1U);
+    EXPECT_EQ(mf->links()[0].href, "href");
+    EXPECT_EQ(mf->links()[0].rel, "rel");
+    EXPECT_EQ(mf->links()[0].type, "type");
+    EXPECT_EQ(mf->links()[0].title, "title");
+    EXPECT_EQ(mf->referenceEpoch(), "2000-01-01T00:00:00Z");
+    EXPECT_EQ(mf->uncertaintyReferenceEpoch(), "2018-12-15T00:00:00Z");
+    EXPECT_EQ(mf->horizontalOffsetUnit(), "metre");
+    EXPECT_EQ(mf->verticalOffsetUnit(), "metre");
+    EXPECT_EQ(mf->horizontalUncertaintyType(), "circular 95% confidence limit");
+    EXPECT_EQ(mf->horizontalUncertaintyUnit(), "metre");
+    EXPECT_EQ(mf->verticalUncertaintyType(), "95% confidence limit");
+    EXPECT_EQ(mf->verticalUncertaintyUnit(), "metre");
+    EXPECT_EQ(mf->horizontalOffsetMethod(), "addition");
+    ASSERT_EQ(mf->components().size(), 6U);
+    {
+        const auto &comp = mf->components()[IDX_CONSTANT];
+        EXPECT_EQ(comp.description(), "description");
+        EXPECT_EQ(comp.displacementType(), "horizontal");
+        EXPECT_EQ(comp.uncertaintyType(), "none");
+        EXPECT_EQ(comp.horizontalUncertainty(), 0.01);
+        EXPECT_EQ(comp.verticalUncertainty(), 0.02);
+        EXPECT_EQ(comp.extent().minx(), modelMinX);
+        EXPECT_EQ(comp.extent().miny(), modelMinY);
+        EXPECT_EQ(comp.extent().maxx(), modelMaxX);
+        EXPECT_EQ(comp.extent().maxy(), modelMaxY);
+        EXPECT_EQ(comp.spatialModel().type, "GeoTIFF");
+        EXPECT_EQ(comp.spatialModel().interpolationMethod, "bilinear");
+        EXPECT_EQ(comp.spatialModel().filename, "nzgd2000-ndm-grid02.tif");
+        EXPECT_EQ(comp.spatialModel().md5Checksum,
+                  "49fce8ab267be2c8d00d43683060a032");
+        ASSERT_NE(comp.timeFunction(), nullptr);
+        ASSERT_EQ(comp.timeFunction()->type, "constant");
+    }
+    {
+        const auto &comp = mf->components()[IDX_VELOCITY];
+        ASSERT_NE(comp.timeFunction(), nullptr);
+        ASSERT_EQ(comp.timeFunction()->type, "velocity");
+        const auto velocity =
+            static_cast<const Component::VelocityTimeFunction *>(
+                comp.timeFunction());
+        EXPECT_EQ(velocity->referenceEpoch.toString(), "2000-01-01T00:00:00Z");
+    }
+    {
+        const auto &comp = mf->components()[IDX_STEP];
+        ASSERT_NE(comp.timeFunction(), nullptr);
+        ASSERT_EQ(comp.timeFunction()->type, "step");
+        const auto step = static_cast<const Component::StepTimeFunction *>(
+            comp.timeFunction());
+        EXPECT_EQ(step->stepEpoch.toString(), "2000-01-01T00:00:00Z");
+    }
+    {
+        const auto &comp = mf->components()[IDX_REVERSE_STEP];
+        ASSERT_NE(comp.timeFunction(), nullptr);
+        ASSERT_EQ(comp.timeFunction()->type, "reverse_step");
+        const auto step =
+            static_cast<const Component::ReverseStepTimeFunction *>(
+                comp.timeFunction());
+        EXPECT_EQ(step->stepEpoch.toString(), "2000-01-01T00:00:00Z");
+    }
+    {
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        ASSERT_NE(comp.timeFunction(), nullptr);
+        ASSERT_EQ(comp.timeFunction()->type, "piecewise");
+        const auto piecewise =
+            static_cast<const Component::PiecewiseTimeFunction *>(
+                comp.timeFunction());
+        EXPECT_EQ(piecewise->beforeFirst, "zero");
+        EXPECT_EQ(piecewise->afterLast, "constant");
+        EXPECT_EQ(piecewise->model.size(), 4U);
+        EXPECT_EQ(piecewise->model[0].epoch.toString(), "2016-01-01T00:00:00Z");
+        EXPECT_EQ(piecewise->model[0].scaleFactor, 0.5);
+    }
+    {
+        const auto &comp = mf->components()[IDX_EXPONENTIAL];
+        ASSERT_NE(comp.timeFunction(), nullptr);
+        ASSERT_EQ(comp.timeFunction()->type, "exponential");
+        const auto exponential =
+            static_cast<const Component::ExponentialTimeFunction *>(
+                comp.timeFunction());
+        EXPECT_EQ(exponential->referenceEpoch.toString(),
+                  "2000-01-01T00:00:00Z");
+        EXPECT_EQ(exponential->endEpoch.toString(), "2001-01-01T00:00:00Z");
+        EXPECT_EQ(exponential->relaxationConstant, 2.0);
+        EXPECT_EQ(exponential->beforeScaleFactor, 0.0);
+        EXPECT_EQ(exponential->initialScaleFactor, 1.0);
+        EXPECT_EQ(exponential->finalScaleFactor, 3.0);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, error_cases) {
+    const auto jFullValid(getFullValidContent());
+    {
+        json jcopy(jFullValid);
+        jcopy["horizontal_offset_method"] = "unsupported";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["horizontal_offset_unit"] = "unsupported";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["vertical_offset_unit"] = "unsupported";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_CONSTANT]["spatial_model"]
+             ["interpolation_method"] = "unsupported";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_CONSTANT]["displacement_type"] = "unsupported";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["model"] = "foo";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["before_first"] = "illegal";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["after_last"] = "illegal";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][0]["time_function"]["type"] = "unknown";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    // Unsupported combination
+    {
+        json jcopy(jFullValid);
+        jcopy["horizontal_offset_method"] = "geocentric";
+        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
+    }
+    {
+        json jcopy(jFullValid);
+        jcopy["horizontal_offset_unit"] = "degree";
+        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
+    }
+    {
+        json jcopy(jFullValid);
+        jcopy["horizontal_offset_method"] = "geocentric";
+        jcopy["horizontal_offset_unit"] = "degree";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+
+    // Unsupported combination
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_VELOCITY]["spatial_model"]
+             ["interpolation_method"] = "geocentric_bilinear";
+        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
+    }
+    {
+        json jcopy(jFullValid);
+        jcopy["horizontal_offset_unit"] = "degree";
+        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
+    }
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_VELOCITY]["spatial_model"]
+             ["interpolation_method"] = "geocentric_bilinear";
+        jcopy["horizontal_offset_unit"] = "degree";
+        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, ISO8601ToDecimalYear) {
+    EXPECT_EQ(ISO8601ToDecimalYear("2000-01-01T00:00:00Z"), 2000.0);
+    EXPECT_EQ(ISO8601ToDecimalYear("2000-02-29T12:00:00Z"),
+              2000.0 + ((31 + 28) * 86400. + 12 * 3600) / (366 * 86400));
+    EXPECT_EQ(ISO8601ToDecimalYear("2000-12-31T23:59:59Z"),
+              2000.0 + (366 * 86400 - 1.) / (366 * 86400));
+    EXPECT_EQ(ISO8601ToDecimalYear("2001-01-01T00:00:00Z"), 2001.0);
+    EXPECT_EQ(ISO8601ToDecimalYear("2001-12-31T23:59:59Z"),
+              2001.0 + (365 * 86400 - 1.) / (365 * 86400));
+    EXPECT_THROW(ISO8601ToDecimalYear(""), ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("0000-01-01T00:00:00Z"),
+                 ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("2001-02-29T00:00:00Z"),
+                 ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("2000-13-01T00:00:00Z"),
+                 ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-32T00:00:00Z"),
+                 ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T24:00:00Z"),
+                 ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T00:60:00Z"),
+                 ParsingException);
+    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T00:00:61Z"),
+                 ParsingException);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluate_constant) {
+    const auto jFullValid(getFullValidContent());
+    const auto mf = MasterFile::parse(jFullValid.dump());
+    const auto &comp = mf->components()[IDX_CONSTANT];
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.0), 1.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.0), 1.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.0), 1.0);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluate_velocity) {
+    const auto jFullValid(getFullValidContent());
+    const auto mf = MasterFile::parse(jFullValid.dump());
+    const auto &comp = mf->components()[IDX_VELOCITY];
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.0), -1.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.0), 0.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.0), 1.0);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluate_step) {
+    const auto jFullValid(getFullValidContent());
+    const auto mf = MasterFile::parse(jFullValid.dump());
+    const auto &comp = mf->components()[IDX_STEP];
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), 0.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 1.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.01), 1.0);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluate_reverse_step) {
+    const auto jFullValid(getFullValidContent());
+    const auto mf = MasterFile::parse(jFullValid.dump());
+    const auto &comp = mf->components()[IDX_REVERSE_STEP];
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), -1.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 0.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.01), 0.0);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluate_piecewise) {
+    const auto jFullValid(getFullValidContent());
+    {
+        const auto mf = MasterFile::parse(jFullValid.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.99), 0.0);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.00), 0.5);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 0.75);
+        EXPECT_NEAR(comp.timeFunction()->evaluateAt(2017 - 1e-9), 1.0, 1e-9);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2017.0), 2.0);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2017.5), 1.5);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.0), 1.0);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2019.0), 1.0);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["before_first"] = "zero";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.0);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["before_first"] = "constant";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["before_first"] = "linear";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.25);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["after_last"] = "zero";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 0.0);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["after_last"] = "constant";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 1.0);
+    }
+
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["after_last"] = "linear";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 0.5);
+    }
+
+    // No epoch
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["model"]
+                 .clear();
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.0);
+    }
+
+    // Just one epoch
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["model"] = {
+                 {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}}};
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["before_first"] = "linear";
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["after_last"] = "linear";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 0.5);
+    }
+
+    // Two identical epochs
+    {
+        json jcopy(jFullValid);
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["model"] = {
+                 {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}},
+                 {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 1.0}}};
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["before_first"] = "linear";
+        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
+             ["after_last"] = "linear";
+        const auto mf = MasterFile::parse(jcopy.dump());
+        const auto &comp = mf->components()[IDX_PIECEWISE];
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5);
+        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 1.0);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluate_exponential) {
+    const auto jFullValid(getFullValidContent());
+    const auto mf = MasterFile::parse(jFullValid.dump());
+    const auto &comp = mf->components()[IDX_EXPONENTIAL];
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), 0.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 1.0);
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.50),
+              1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2000.50 - 2000.00) / 2.0)));
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.00),
+              1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2001.00 - 2000.00) / 2.0)));
+    EXPECT_EQ(comp.timeFunction()->evaluateAt(2002.00),
+              1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2001.00 - 2000.00) / 2.0)));
+}
+
+// ---------------------------------------------------------------------------
+
+inline double RadToDeg(double d) { return d / DEG_TO_RAD_CONSTANT; }
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluator_horizontal_unit_degree) {
+    json j(getMinValidContent());
+    j["horizontal_offset_method"] = "addition";
+    j["horizontal_offset_unit"] = "degree";
+    constexpr double tFactor = 0.5;
+    constexpr double gridMinX = 160;
+    constexpr double gridMinY = -50;
+    constexpr double gridMaxX = 190;
+    constexpr double gridMaxY = -30;
+    j["components"] = {
+        {{"displacement_type", "horizontal"},
+         {"uncertainty_type", "none"},
+         {"extent",
+          {{"type", "bbox"},
+           {"parameters",
+            {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}}}},
+         {"spatial_model",
+          {
+              {"type", "GeoTIFF"},
+              {"interpolation_method", "bilinear"},
+              {"filename", "bla.tif"},
+          }},
+         {"time_function",
+          {{"type", "piecewise"},
+           {"parameters",
+            {{"before_first", "zero"},
+             {"after_last", "zero"},
+             {"model",
+              {{{"epoch", "2010-01-01T00:00:00Z"}, {"scale_factor", tFactor}},
+               {{"epoch", "2020-01-01T00:00:00Z"},
+                {"scale_factor", tFactor}}}}}}}}}};
+
+    constexpr int iQueriedX = 1;
+    constexpr int iQueriedY = 3;
+    constexpr double lonOffsetQueriedX = 0.01;
+    constexpr double lonOffsetQueriedXp1 = 0.02;
+    constexpr double latOffsetQueriedY = 0.03;
+    constexpr double latOffsetQueriedYp1 = 0.04;
+    constexpr double zOffsetQueriedXY = 10.;
+    constexpr double zOffsetQueriedXp1Y = 11.;
+    constexpr double zOffsetQueriedXYp1 = 11.;
+    constexpr double zOffsetQueriedXp1Yp1 = 12.;
+    constexpr double gridResX = 2;
+    constexpr double gridResY = 0.5;
+
+    struct Grid : public GridPrototype {
+        bool getLonLatOffset(int ix, int iy, double &lonOffsetRadian,
+                             double &latOffsetRadian) const {
+            if (ix == iQueriedX) {
+                lonOffsetRadian = DegToRad(lonOffsetQueriedX);
+            } else if (ix == iQueriedX + 1) {
+                lonOffsetRadian = DegToRad(lonOffsetQueriedXp1);
+            } else {
+                return false;
+            }
+            if (iy == iQueriedY) {
+                latOffsetRadian = DegToRad(latOffsetQueriedY);
+            } else if (iy == iQueriedY + 1) {
+                latOffsetRadian = DegToRad(latOffsetQueriedYp1);
+            } else {
+                return false;
+            }
+            return true;
+        }
+
+        bool getZOffset(int ix, int iy, double &zOffset) const {
+            if (ix == iQueriedX && iy == iQueriedY) {
+                zOffset = zOffsetQueriedXY;
+            } else if (ix == iQueriedX + 1 && iy == iQueriedY) {
+                zOffset = zOffsetQueriedXp1Y;
+            } else if (ix == iQueriedX && iy == iQueriedY + 1) {
+                zOffset = zOffsetQueriedXYp1;
+            } else if (ix == iQueriedX + 1 && iy == iQueriedY + 1) {
+                zOffset = zOffsetQueriedXp1Yp1;
+            } else {
+                return false;
+            }
+            return true;
+        }
+
+        bool getLonLatZOffset(int ix, int iy, double &lonOffsetRadian,
+                              double &latOffsetRadian, double &zOffset) const {
+            return getLonLatOffset(ix, iy, lonOffsetRadian, latOffsetRadian) &&
+                   getZOffset(ix, iy, zOffset);
+        }
+
+#ifdef DEBUG_DEFMODEL
+        std::string name() const { return std::string(); }
+#endif
+    };
+
+    struct GridSet : public GridSetPrototype<Grid> {
+
+        Grid grid{};
+
+        GridSet() {
+            grid.minx = DegToRad(gridMinX);
+            grid.miny = DegToRad(gridMinY);
+            grid.resx = DegToRad(gridResX);
+            grid.resy = DegToRad(gridResY);
+            grid.width =
+                1 + static_cast<int>(0.5 + (gridMaxX - gridMinX) / gridResX);
+            grid.height =
+                1 + static_cast<int>(0.5 + (gridMaxY - gridMinY) / gridResY);
+        }
+
+        const Grid *gridAt(double /*x */, double /* y */) { return &grid; }
+    };
+
+    struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
+        std::unique_ptr<GridSet> open(const std::string &filename) {
+            if (filename != "bla.tif")
+                return nullptr;
+            return std::unique_ptr<GridSet>(new GridSet());
+        }
+
+        bool isGeographicCRS(const std::string & /* crsDef */) { return true; }
+
+#ifdef DEBUG_DEFMODEL
+        void log(const std::string & /* msg */) {}
+#endif
+    };
+
+    EvaluatorIface iface;
+
+    Evaluator<Grid, GridSet, EvaluatorIface> eval(MasterFile::parse(j.dump()),
+                                                  iface, 1, 1);
+    double newLon;
+    double newLat;
+    double newZ;
+    constexpr double tValid = 2018;
+    constexpr double EPS = 1e-9;
+    constexpr double zVal = 100;
+
+    // Query on exact grid intersection
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon + tFactor * lonOffsetQueriedX, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat + tFactor * latOffsetQueriedY, EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Query between grid points
+    {
+        constexpr double alphaX = 0.25;
+        constexpr double alphaY = 0.125;
+        const double lon = gridMinX + iQueriedX * gridResX + alphaX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(
+            RadToDeg(newLon),
+            lon +
+                tFactor * (lonOffsetQueriedX +
+                           alphaX * (lonOffsetQueriedXp1 - lonOffsetQueriedX)),
+            EPS);
+        EXPECT_NEAR(
+            RadToDeg(newLat),
+            lat +
+                tFactor * (latOffsetQueriedY +
+                           alphaY * (latOffsetQueriedYp1 - latOffsetQueriedY)),
+            EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Longitude < model min
+    {
+        const double lon = modelMinX - 1e-1;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                  tValid, newLon, newLat, newZ));
+    }
+
+    // Longitude > model max
+    {
+        const double lon = modelMaxX + 1e-1;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                  tValid, newLon, newLat, newZ));
+    }
+
+    // Latitude < model min
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = modelMinY - 1e-1;
+        EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                  tValid, newLon, newLat, newZ));
+    }
+
+    // Latitude > model max
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = modelMaxY + 1e-1;
+        EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                  tValid, newLon, newLat, newZ));
+    }
+
+    // Before timeExtent.first
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                  1000, newLon, newLat, newZ));
+    }
+
+    // After timeExtent.last
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_FALSE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                  3000, newLon, newLat, newZ));
+    }
+
+    // Longitude < grid min
+    {
+        const double lon = gridMinX - 1e-1;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Longitude > grid max
+    {
+        const double lon = gridMaxX + 1e-1;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Latitude < grid min
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = gridMinY - 1e-1;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Latitude > grid max
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = gridMaxY + 1e-1;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Time function values to zero
+    {
+        const double lon = gridMinX + iQueriedX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 2000, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
+        EXPECT_EQ(newZ, zVal);
+    }
+
+    // Test vertical
+    j["components"][0]["displacement_type"] = "vertical";
+    j["vertical_offset_unit"] = "metre";
+    Evaluator<Grid, GridSet, EvaluatorIface> evalVertical(
+        MasterFile::parse(j.dump()), iface, 1, 1);
+    {
+        constexpr double alphaX = 0.25;
+        constexpr double alphaY = 0.125;
+        const double lon = gridMinX + iQueriedX * gridResX + alphaX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY;
+        EXPECT_TRUE(evalVertical.forward(iface, DegToRad(lon), DegToRad(lat),
+                                         zVal, tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(RadToDeg(newLon), lon, EPS);
+        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
+
+        const double zBottom =
+            zOffsetQueriedXY + alphaX * (zOffsetQueriedXp1Y - zOffsetQueriedXY);
+        const double zTop =
+            zOffsetQueriedXYp1 +
+            alphaX * (zOffsetQueriedXp1Yp1 - zOffsetQueriedXYp1);
+        EXPECT_NEAR(
+            newZ, zVal + tFactor * (zBottom + alphaY * (zTop - zBottom)), EPS);
+    }
+
+    // Test 3d
+    j["components"][0]["displacement_type"] = "3d";
+    j["vertical_offset_unit"] = "metre";
+    Evaluator<Grid, GridSet, EvaluatorIface> eval3d(MasterFile::parse(j.dump()),
+                                                    iface, 1, 1);
+    {
+        constexpr double alphaX = 0.25;
+        constexpr double alphaY = 0.125;
+        const double lon = gridMinX + iQueriedX * gridResX + alphaX * gridResX;
+        const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY;
+        EXPECT_TRUE(eval3d.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                   tValid, newLon, newLat, newZ));
+        EXPECT_NEAR(
+            RadToDeg(newLon),
+            lon +
+                tFactor * (lonOffsetQueriedX +
+                           alphaX * (lonOffsetQueriedXp1 - lonOffsetQueriedX)),
+            EPS);
+        EXPECT_NEAR(
+            RadToDeg(newLat),
+            lat +
+                tFactor * (latOffsetQueriedY +
+                           alphaY * (latOffsetQueriedYp1 - latOffsetQueriedY)),
+            EPS);
+
+        const double zBottom =
+            zOffsetQueriedXY + alphaX * (zOffsetQueriedXp1Y - zOffsetQueriedXY);
+        const double zTop =
+            zOffsetQueriedXYp1 +
+            alphaX * (zOffsetQueriedXp1Yp1 - zOffsetQueriedXYp1);
+        EXPECT_NEAR(
+            newZ, zVal + tFactor * (zBottom + alphaY * (zTop - zBottom)), EPS);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+inline void DeltaLongLatToEastingNorthing(double phi, double dlam, double dphi,
+                                          double a, double b, double &de,
+                                          double &dn) {
+    const double sinphi = sin(phi);
+    const double cosphi = cos(phi);
+    const double a2 = a * a;
+    const double b2 = b * b;
+    const double X = a2 * (cosphi * cosphi) + b2 * (sinphi * sinphi);
+    const double sqrtX = sqrt(X);
+    de = dlam * (a2 * cosphi) / sqrtX;
+    dn = dphi * a2 * b2 / (sqrtX * X);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluator_horizontal_unit_metre) {
+
+    json j(getMinValidContent());
+    j["horizontal_offset_method"] = "addition";
+    j["horizontal_offset_unit"] = "metre";
+    j["vertical_offset_unit"] = "metre";
+    constexpr double tFactor = 0.5;
+    constexpr double gridMinX = 165.8;
+    constexpr double gridMinY = -37.5;
+    constexpr double gridMaxX = 166.2;
+    constexpr double gridMaxY = -37.2;
+    constexpr double gridResX = gridMaxX - gridMinX;
+    constexpr double gridResY = gridMaxY - gridMinY;
+    constexpr int extraPointX = 1;
+    constexpr int extraPointY = 1;
+
+    j["components"] = {
+        {{"displacement_type", "horizontal"},
+         {"uncertainty_type", "none"},
+         {"extent",
+          {{"type", "bbox"},
+           {"parameters",
+            {{"bbox",
+              {gridMinX - extraPointX * gridResX,
+               gridMinY - extraPointY * gridResY, gridMaxX, gridMaxY}}}}}},
+         {"spatial_model",
+          {
+              {"type", "GeoTIFF"},
+              {"interpolation_method", "XXXXXXX"},
+              {"filename", "bla.tif"},
+          }},
+         {"time_function",
+          {{"type", "piecewise"},
+           {"parameters",
+            {{"before_first", "zero"},
+             {"after_last", "zero"},
+             {"model",
+              {{{"epoch", "2010-01-01T00:00:00Z"}, {"scale_factor", tFactor}},
+               {{"epoch", "2020-01-01T00:00:00Z"},
+                {"scale_factor", tFactor}}}}}}}}}};
+
+    struct Grid : public GridPrototype {
+        bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset,
+                                      double &northingOffset) const {
+            ix -= extraPointX;
+            iy -= extraPointY;
+            if (ix == -1)
+                ix = 0;
+            if (iy == -1)
+                iy = 0;
+            if (ix == 0 && iy == 0) {
+                eastingOffset = 0.4f;
+                northingOffset = -0.2f;
+            } else if (ix == 1 && iy == 0) {
+                eastingOffset = 0.5f;
+                northingOffset = -0.25f;
+            } else if (ix == 0 && iy == 1) {
+                eastingOffset = 0.8f;
+                northingOffset = -0.4f;
+            } else if (ix == 1 && iy == 1) {
+                eastingOffset = 1.f;
+                northingOffset = -0.3f;
+            } else {
+                return false;
+            }
+            return true;
+        }
+
+        bool getZOffset(int ix, int iy, double &zOffset) const {
+            ix -= extraPointX;
+            iy -= extraPointY;
+            if (ix == -1)
+                ix = 0;
+            if (iy == -1)
+                iy = 0;
+            if (ix == 0 && iy == 0) {
+                zOffset = 0.84f;
+            } else if (ix == 1 && iy == 0) {
+                zOffset = 0.75f;
+            } else if (ix == 0 && iy == 1) {
+                zOffset = 0.36f;
+            } else if (ix == 1 && iy == 1) {
+                zOffset = 0.f;
+            } else {
+                return false;
+            }
+            return true;
+        }
+
+        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 std::string(); }
+#endif
+    };
+
+    struct GridSet : public GridSetPrototype<Grid> {
+
+        Grid grid{};
+
+        GridSet() {
+            grid.minx = DegToRad(gridMinX - extraPointX * gridResX);
+            grid.miny = DegToRad(gridMinY - extraPointY * gridResY);
+            grid.resx = DegToRad(gridResX);
+            grid.resy = DegToRad(gridResY);
+            grid.width = 2 + extraPointX;
+            grid.height = 2 + extraPointY;
+        }
+
+        const Grid *gridAt(double /*x */, double /* y */) { return &grid; }
+    };
+
+    struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
+        std::unique_ptr<GridSet> open(const std::string &filename) {
+            if (filename != "bla.tif")
+                return nullptr;
+            return std::unique_ptr<GridSet>(new GridSet());
+        }
+
+        bool isGeographicCRS(const std::string & /* crsDef */) { return true; }
+
+#ifdef DEBUG_DEFMODEL
+        void log(const std::string & /* msg */) {}
+#endif
+
+        void geographicToGeocentric(double lam, double phi, double height,
+                                    double a, double /*b*/, double es,
+                                    double &X, double &Y, double &Z) {
+            PJ_CONTEXT *ctx = proj_context_create();
+            PJ *cart = proj_create(
+                ctx,
+                ("+proj=cart +a=" + osgeo::proj::internal::toString(a, 18) +
+                 " +es=" + osgeo::proj::internal::toString(es, 18))
+                    .c_str());
+            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;
+            proj_destroy(cart);
+            proj_context_destroy(ctx);
+        }
+
+        void geocentricToGeographic(double X, double Y, double Z, double a,
+                                    double /*b*/, double es, double &lam,
+                                    double &phi, double &height) {
+            PJ_CONTEXT *ctx = proj_context_create();
+            PJ *cart = proj_create(
+                ctx,
+                ("+proj=cart +a=" + osgeo::proj::internal::toString(a, 18) +
+                 " +es=" + osgeo::proj::internal::toString(es, 18))
+                    .c_str());
+            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;
+            proj_destroy(cart);
+            proj_context_destroy(ctx);
+        }
+    };
+    EvaluatorIface iface;
+
+    constexpr double a = 6378137;
+    constexpr double b = 6356752.314140;
+    constexpr double tValid = 2018;
+    constexpr double zVal = 100;
+
+    const struct {
+        double lon;
+        double lat;
+        double expected_de;
+        double expected_dn;
+        double expected_dz;
+        const char *displacement_type;
+        const char *interpolation_method;
+    } testPoints[] = {
+        {gridMinX - extraPointX * gridResX - 1e-11,
+         gridMinY - extraPointY * gridResY - 1e-11, 0.4, -0.2, 0, "horizontal",
+         "bilinear"},
+        {gridMinX, gridMinY, 0.4, -0.2, 0, "horizontal", "bilinear"},
+        {gridMaxX, gridMinY, 0.5, -0.25, 0, "horizontal", "bilinear"},
+        {gridMinX, gridMaxY, 0.8, -0.4, 0, "horizontal", "bilinear"},
+        {gridMaxX, gridMaxY, 1, -0.3, 0, "horizontal", "bilinear"},
+        {gridMaxX + 1e-11, gridMaxY + 1e-11, 1, -0.3, 0, "horizontal",
+         "bilinear"},
+        {165.9, -37.3, 0.70833334, -0.32083334, 0, "horizontal", "bilinear"},
+        {165.9, -37.3, 0.70833334, -0.32083334, 0.4525, "3d", "bilinear"},
+
+        {gridMinX, gridMinY, 0.4, -0.2, 0, "horizontal", "geocentric_bilinear"},
+        {gridMaxX, gridMinY, 0.5, -0.25, 0, "horizontal",
+         "geocentric_bilinear"},
+        {gridMinX, gridMaxY, 0.8, -0.4, 0, "horizontal", "geocentric_bilinear"},
+        {gridMaxX, gridMaxY, 1, -0.3, 0, "horizontal", "geocentric_bilinear"},
+        {165.9, -37.3, 0.7083692044608846, -0.3209642339711405, 0, "horizontal",
+         "geocentric_bilinear"},
+        {165.9, -37.3, 0.7083692044608846, -0.3209642339711405, 0.4525, "3d",
+         "geocentric_bilinear"},
+    };
+
+    for (const auto &testPoint : testPoints) {
+        j["components"][0]["displacement_type"] = testPoint.displacement_type;
+        j["components"][0]["spatial_model"]["interpolation_method"] =
+            testPoint.interpolation_method;
+        Evaluator<Grid, GridSet, EvaluatorIface> eval(
+            MasterFile::parse(j.dump()), iface, a, b);
+
+        const double lon = testPoint.lon;
+        const double lat = testPoint.lat;
+        double newLon;
+        double newLat;
+        double newZ;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ))
+            << lon << " " << lat << " " << testPoint.displacement_type
+            << testPoint.interpolation_method;
+        EXPECT_NEAR(newZ - zVal, tFactor * testPoint.expected_dz, 1e-8)
+            << lon << " " << lat << " " << testPoint.displacement_type
+            << testPoint.interpolation_method;
+
+        double de;
+        double dn;
+        DeltaLongLatToEastingNorthing(DegToRad(lat), newLon - DegToRad(lon),
+                                      newLat - DegToRad(lat), a, b, de, dn);
+        EXPECT_NEAR(de, tFactor * testPoint.expected_de, 1e-8)
+            << lon << " " << lat << " " << testPoint.displacement_type
+            << testPoint.interpolation_method;
+        EXPECT_NEAR(dn, tFactor * testPoint.expected_dn, 1e-8)
+            << lon << " " << lat << " " << testPoint.displacement_type
+            << testPoint.interpolation_method;
+
+        if (lon == gridMinX && lat == gridMinY) {
+            // Redo the exact same test, to test caching
+            double newLon2;
+            double newLat2;
+            double newZ2;
+            EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                     tValid, newLon2, newLat2, newZ2));
+            EXPECT_EQ(newLon2, newLon);
+            EXPECT_EQ(newLat2, newLat);
+            EXPECT_EQ(newZ2, newZ);
+
+            // Shift in longitude
+            EXPECT_TRUE(eval.forward(iface, DegToRad(lon - gridResX / 2),
+                                     DegToRad(lat), zVal, tValid, newLon2,
+                                     newLat2, newZ2));
+
+            // Redo test at original position
+            EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                     tValid, newLon2, newLat2, newZ2));
+            EXPECT_EQ(newLon2, newLon);
+            EXPECT_EQ(newLat2, newLat);
+            EXPECT_EQ(newZ2, newZ);
+
+            // Shift in latitude
+            EXPECT_TRUE(eval.forward(iface, DegToRad(lon),
+                                     DegToRad(lat - gridResY / 2), zVal, tValid,
+                                     newLon2, newLat2, newZ2));
+
+            // Redo test at original position
+            EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                     tValid, newLon2, newLat2, newZ2));
+            EXPECT_EQ(newLon2, newLon);
+            EXPECT_EQ(newLat2, newLat);
+            EXPECT_EQ(newZ2, newZ);
+        }
+    }
+
+    // Test inverse()
+    {
+        j["horizontal_offset_method"] = "addition";
+        j["components"][0]["displacement_type"] = "3d";
+        j["components"][0]["spatial_model"]["interpolation_method"] =
+            "bilinear";
+        Evaluator<Grid, GridSet, EvaluatorIface> eval(
+            MasterFile::parse(j.dump()), iface, a, b);
+
+        const double lon = 165.9;
+        const double lat = -37.3;
+        double newLon;
+        double newLat;
+        double newZ;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+
+        double invLon;
+        double invLat;
+        double invZ;
+        EXPECT_TRUE(eval.inverse(iface, newLon, newLat, newZ, tValid, invLon,
+                                 invLat, invZ));
+        EXPECT_NEAR(RadToDeg(invLon), lon, 1e-10);
+        EXPECT_NEAR(RadToDeg(invLat), lat, 1e-10);
+        EXPECT_NEAR(invZ, zVal, 1e-4);
+    }
+
+    // Test horizontal_offset_method = geocentric
+    {
+        j["horizontal_offset_method"] = "geocentric";
+        j["components"][0]["displacement_type"] = "3d";
+        j["components"][0]["spatial_model"]["interpolation_method"] =
+            "bilinear";
+        Evaluator<Grid, GridSet, EvaluatorIface> eval(
+            MasterFile::parse(j.dump()), iface, a, b);
+
+        const double lon = gridMinX;
+        const double lat = gridMinY;
+        double newLon;
+        double newLat;
+        double newZ;
+        EXPECT_TRUE(eval.forward(iface, DegToRad(lon), DegToRad(lat), zVal,
+                                 tValid, newLon, newLat, newZ));
+
+        double de;
+        double dn;
+        constexpr double expected_de = 0.40000000948081327;
+        constexpr double expected_dn = -0.19999999810542682;
+        DeltaLongLatToEastingNorthing(DegToRad(lat), newLon - DegToRad(lon),
+                                      newLat - DegToRad(lat), a, b, de, dn);
+        EXPECT_NEAR(de, tFactor * expected_de, 1e-10);
+        EXPECT_NEAR(dn, tFactor * expected_dn, 1e-10);
+        EXPECT_NEAR(newZ - zVal, tFactor * 0.84, 1e-4);
+    }
+}
+
+// ---------------------------------------------------------------------------
+
+TEST(defmodel, evaluator_projected_crs) {
+
+    json j(getMinValidContent());
+    j["horizontal_offset_method"] = "addition";
+    j["horizontal_offset_unit"] = "metre";
+    j["vertical_offset_unit"] = "metre";
+    constexpr double gridMinX = 10000;
+    constexpr double gridMinY = 20000;
+    constexpr double gridMaxX = 30000;
+    constexpr double gridMaxY = 40000;
+    constexpr double gridResX = gridMaxX - gridMinX;
+    constexpr double gridResY = gridMaxY - gridMinY;
+
+    j["extent"]
+     ["parameters"] = {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}};
+    j["components"] = {
+        {{"displacement_type", "horizontal"},
+         {"uncertainty_type", "none"},
+         {"extent",
+          {{"type", "bbox"},
+           {"parameters",
+            {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}}}},
+         {"spatial_model",
+          {
+              {"type", "GeoTIFF"},
+              {"interpolation_method", "bilinear"},
+              {"filename", "bla.tif"},
+          }},
+         {"time_function", {{"type", "constant"}}}}};
+
+    struct Grid : public GridPrototype {
+        bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset,
+                                      double &northingOffset) const {
+            if (ix == 0 && iy == 0) {
+                eastingOffset = 0.4;
+                northingOffset = -0.2;
+            } else if (ix == 1 && iy == 0) {
+                eastingOffset = 0.5;
+                northingOffset = -0.25;
+            } else if (ix == 0 && iy == 1) {
+                eastingOffset = 0.8;
+                northingOffset = -0.4;
+            } else if (ix == 1 && iy == 1) {
+                eastingOffset = 1.;
+                northingOffset = -0.3;
+            } else {
+                return false;
+            }
+            return true;
+        }
+
+#ifdef DEBUG_DEFMODEL
+        std::string name() const { return std::string(); }
+#endif
+    };
+
+    struct GridSet : public GridSetPrototype<Grid> {
+
+        Grid grid{};
+
+        GridSet() {
+            grid.minx = gridMinX;
+            grid.miny = gridMinY;
+            grid.resx = gridResX;
+            grid.resy = gridResY;
+            grid.width = 2;
+            grid.height = 2;
+        }
+
+        const Grid *gridAt(double /*x */, double /* y */) { return &grid; }
+    };
+
+    struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
+        std::unique_ptr<GridSet> open(const std::string &filename) {
+            if (filename != "bla.tif")
+                return nullptr;
+            return std::unique_ptr<GridSet>(new GridSet());
+        }
+
+        bool isGeographicCRS(const std::string & /* crsDef */) { return false; }
+
+#ifdef DEBUG_DEFMODEL
+        void log(const std::string & /* msg */) {}
+#endif
+    };
+    EvaluatorIface iface;
+
+    constexpr double a = 6378137;
+    constexpr double b = 6356752.314140;
+    constexpr double tValid = 2018;
+    constexpr double zVal = 100;
+
+    Evaluator<Grid, GridSet, EvaluatorIface> eval(MasterFile::parse(j.dump()),
+                                                  iface, a, b);
+
+    double newX;
+    double newY;
+    double newZ;
+    EXPECT_TRUE(eval.forward(iface, gridMinX, gridMinY, zVal, tValid, newX,
+                             newY, newZ));
+    EXPECT_NEAR(newX - gridMinX, 0.4, 1e-8);
+    EXPECT_NEAR(newY - gridMinY, -0.2, 1e-8);
+    EXPECT_NEAR(newZ - zVal, 0, 1e-8);
+
+    {
+        json jcopy(j);
+        jcopy["horizontal_offset_unit"] = "degree";
+        EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>(
+                         MasterFile::parse(jcopy.dump()), iface, a, b)),
+                     EvaluatorException);
+    }
+
+    {
+        json jcopy(j);
+        jcopy["horizontal_offset_method"] = "geocentric";
+        EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>(
+                         MasterFile::parse(jcopy.dump()), iface, a, b)),
+                     EvaluatorException);
+    }
+
+    {
+        json jcopy(j);
+        jcopy["components"][0]["spatial_model"]["interpolation_method"] =
+            "geocentric_bilinear";
+        EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>(
+                         MasterFile::parse(jcopy.dump()), iface, a, b)),
+                     EvaluatorException);
+    }
+}
+
+} // namespace
+
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
diff --git a/test/unit/test_grids.cpp b/test/unit/test_grids.cpp
index 1c77e5b6..e9841acd 100644
--- a/test/unit/test_grids.cpp
+++ b/test/unit/test_grids.cpp
@@ -232,17 +232,17 @@ TEST_F(GridTest, GenericShiftGridSet_gtiff_projected) {
         m_ctxt, "tests/test_3d_grid_projected.tif");
     ASSERT_NE(gridSet, nullptr);
     ASSERT_EQ(gridSet->gridAt(-1000, -1000), nullptr);
-    auto grid = gridSet->gridAt(0, 0);
+    auto grid = gridSet->gridAt(1500300.0, 5400300.0);
     ASSERT_NE(grid, nullptr);
     EXPECT_EQ(grid->width(), 2);
     EXPECT_EQ(grid->height(), 2);
     EXPECT_EQ(grid->extentAndRes().isGeographic, false);
-    EXPECT_EQ(grid->extentAndRes().west, -180);
-    EXPECT_EQ(grid->extentAndRes().east, 180);
-    EXPECT_EQ(grid->extentAndRes().south, -90);
-    EXPECT_EQ(grid->extentAndRes().north, 90);
-    EXPECT_EQ(grid->extentAndRes().resX, 360);
-    EXPECT_EQ(grid->extentAndRes().resY, 180);
+    EXPECT_EQ(grid->extentAndRes().west, 1500000.0);
+    EXPECT_EQ(grid->extentAndRes().east, 1501000.0);
+    EXPECT_EQ(grid->extentAndRes().south, 5400000.0);
+    EXPECT_EQ(grid->extentAndRes().north, 5401000.0);
+    EXPECT_EQ(grid->extentAndRes().resX, 1000);
+    EXPECT_EQ(grid->extentAndRes().resY, 1000);
 }
 
 } // namespace