diff options
Diffstat (limited to 'test')
| -rw-r--r-- | test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/gie/Makefile.am | 9 | ||||
| -rw-r--r-- | test/gie/builtins.gie | 124 | ||||
| -rw-r--r-- | test/gie/defmodel.gie | 124 | ||||
| -rw-r--r-- | test/unit/CMakeLists.txt | 11 | ||||
| -rw-r--r-- | test/unit/Makefile.am | 11 | ||||
| -rw-r--r-- | test/unit/test_defmodel.cpp | 1515 | ||||
| -rw-r--r-- | test/unit/test_grids.cpp | 31 | ||||
| -rw-r--r-- | test/unit/test_io.cpp | 168 | ||||
| -rw-r--r-- | test/unit/test_operation.cpp | 372 |
10 files changed, 2356 insertions, 10 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/builtins.gie b/test/gie/builtins.gie index 6e5b326f..76e09ca9 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -2158,6 +2158,7 @@ expect 135 50 =============================================================================== # Interrupted Goode Homolosine # PCyl, Sph. +# (Each of the 12 sub-projections tested separately) =============================================================================== ------------------------------------------------------------------------------- @@ -2166,12 +2167,53 @@ operation +proj=igh +a=6400000 tolerance 0.1 mm accept 2 1 expect 223878.497456271 111701.072127637 +roundtrip 1 accept 2 -1 expect 223708.371313058 -111701.072127637 +roundtrip 1 accept -2 1 expect -222857.740596992 111701.072127637 +roundtrip 1 accept -2 -1 expect -223027.866740205 -111701.072127637 +roundtrip 1 + +accept -100.0 22.0 +expect -11170107.212763708 2457423.5868080168 +roundtrip 1 +accept -30.0 22.0 +expect -2863013.673043605 2457423.586808016 +roundtrip 1 +accept -100.0 67.0 +expect -11170107.212763708 7205942.523056464 +roundtrip 1 +accept -30.0 67.0 +expect 17045.719482862 7205942.523056464 +roundtrip 1 +accept -160.0 -22.0 +expect -17872171.540421933 -2457423.586808016 +roundtrip 1 +accept -60.0 -22.0 +expect -6702064.327658225 -2457423.586808016 +roundtrip 1 +accept 20.0 -22.0 +expect 2234021.442552742 -2457423.586808016 +roundtrip 1 +accept 140.0 -22.0 +expect 15638150.097869191 -2457423.586808016 +roundtrip 1 +accept -160.0 -67.0 +expect -17872171.540421933 -7205942.523056464 +roundtrip 1 +accept -60.0 -67.0 +expect -6702064.327658225 -7205942.523056464 +roundtrip 1 +accept 20.0 -67.0 +expect 2234021.442552742 -7205942.523056464 +roundtrip 1 +accept 140.0 -67.0 +expect 15638150.097869191 -7205942.523056464 +roundtrip 1 direction inverse accept 200 100 @@ -2183,6 +2225,84 @@ expect -0.001790497 0.000895247 accept -200 -100 expect -0.001790496 -0.000895247 +=============================================================================== +# Interrupted Goode Homolosine Ocean View +# PCyl, Sph. +# (Each of the 12 sub-projections tested separately) +=============================================================================== + +------------------------------------------------------------------------------- +operation +proj=igh_o +a=6400000 +------------------------------------------------------------------------------- +tolerance 0.1 mm +accept 2 1 +expect 223197.992883418 111701.072127637 +roundtrip 1 +accept 2 -1 +expect 223708.371313058 -111701.072127637 +roundtrip 1 +accept -2 1 +expect -223538.245169845 111701.072127637 +roundtrip 1 +accept -2 -1 +expect -223027.866740205 -111701.072127637 +roundtrip 1 + +accept -140.0 22.0 +expect -15638150.097869192 2457423.586808016 +roundtrip 1 +accept 170.0 70.0 +expect 16560870.317293623 7463176.386461447 +roundtrip 1 +accept -10.0 22.0 +expect -1117010.721276371 2457423.586808016 +roundtrip 1 +accept 130.0 22.0 +expect 14521139.376592822 2457423.586808016 +roundtrip 1 +accept -170.0 70.0 +expect -17167948.303394791 7463176.386461447 +roundtrip 1 +accept -140.0 67.0 +expect -15638150.097869191 7205942.523056464 +roundtrip 1 +accept -10.0 67.0 +expect -1117010.721276371 7205942.523056464 +roundtrip 1 +accept 130.0 67.0 +expect 14521139.376592822 7205942.523056464 +roundtrip 1 +accept -110.0 -22.0 +expect -12287117.934040081 -2457423.586808016 +roundtrip 1 +accept 20.0 -22.0 +expect 2234021.442552742 -2457423.586808016 +roundtrip 1 +accept 150.0 -22.0 +expect 16755160.819145568 -2457423.586808016 +roundtrip 1 +accept -110.0 -67.0 +expect -12287117.934040081 -7205942.523056464 +roundtrip 1 +accept 20.0 -67.0 +expect 2234021.442552742 -7205942.523056464 +roundtrip 1 +accept 95.0 -67.0 +expect 13699006.578494834 -7205942.523056464 +roundtrip 1 +accept 150.0 -67.0 +expect 16755160.819145564 -7205942.523056464 +roundtrip 1 + +direction inverse +accept 200 100 +expect 0.001790494 0.000895247 +accept 200 -100 +expect 0.001790491 -0.000895247 +accept -200 100 +expect -0.001790492 0.000895247 +accept -200 -100 +expect -0.001790496 -0.000895247 =============================================================================== # International Map of the World Polyconic @@ -3411,12 +3531,16 @@ operation +proj=moll +a=6400000 tolerance 0.1 mm accept 2 1 expect 201113.698641813 124066.283433860 +roundtrip 1 accept 2 -1 expect 201113.698641813 -124066.283433860 +roundtrip 1 accept -2 1 expect -201113.698641813 124066.283433860 +roundtrip 1 accept -2 -1 expect -201113.698641813 -124066.283433860 +roundtrip 1 direction inverse accept 200 100 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 5240949e..e9841acd 100644 --- a/test/unit/test_grids.cpp +++ b/test/unit/test_grids.cpp @@ -66,7 +66,7 @@ TEST_F(GridTest, VerticalShiftGridSet_null) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 3); EXPECT_EQ(grid->height(), 3); - EXPECT_EQ(grid->extentAndRes().westLon, -M_PI); + EXPECT_EQ(grid->extentAndRes().west, -M_PI); EXPECT_TRUE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out = -1.0f; @@ -103,7 +103,7 @@ TEST_F(GridTest, HorizontalShiftGridSet_null) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 3); EXPECT_EQ(grid->height(), 3); - EXPECT_EQ(grid->extentAndRes().westLon, -M_PI); + EXPECT_EQ(grid->extentAndRes().west, -M_PI); EXPECT_TRUE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out1 = -1.0f; @@ -131,7 +131,7 @@ TEST_F(GridTest, HorizontalShiftGridSet_gtiff) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 4); EXPECT_EQ(grid->height(), 4); - EXPECT_EQ(grid->extentAndRes().westLon, 4.0 / 180 * M_PI); + EXPECT_EQ(grid->extentAndRes().west, 4.0 / 180 * M_PI); EXPECT_FALSE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out1 = -1.0f; @@ -154,7 +154,7 @@ TEST_F(GridTest, GenericShiftGridSet_null) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 3); EXPECT_EQ(grid->height(), 3); - EXPECT_EQ(grid->extentAndRes().westLon, -M_PI); + EXPECT_EQ(grid->extentAndRes().west, -M_PI); EXPECT_TRUE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out = -1.0f; @@ -180,7 +180,8 @@ TEST_F(GridTest, GenericShiftGridSet_gtiff) { ASSERT_NE(grid, nullptr); EXPECT_EQ(grid->width(), 5); EXPECT_EQ(grid->height(), 5); - EXPECT_EQ(grid->extentAndRes().westLon, 21.0 / 180 * M_PI); + EXPECT_EQ(grid->extentAndRes().isGeographic, true); + EXPECT_EQ(grid->extentAndRes().west, 21.0 / 180 * M_PI); EXPECT_FALSE(grid->isNullGrid()); EXPECT_FALSE(grid->hasChanged()); float out = -1.0f; @@ -224,4 +225,24 @@ TEST_F(GridTest, EXPECT_EQ(grid->height(), 8); } +// --------------------------------------------------------------------------- + +TEST_F(GridTest, GenericShiftGridSet_gtiff_projected) { + auto gridSet = NS_PROJ::GenericShiftGridSet::open( + m_ctxt, "tests/test_3d_grid_projected.tif"); + ASSERT_NE(gridSet, nullptr); + ASSERT_EQ(gridSet->gridAt(-1000, -1000), nullptr); + 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, 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 diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 7ad9cf4a..17fddaaf 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -2484,6 +2484,174 @@ TEST(wkt_parse, COMPD_CS_non_conformant_horizontal_plus_horizontal_as_in_LAS) { // --------------------------------------------------------------------------- +TEST(wkt_parse, + COMPD_CS_non_conformant_horizontal_TOWGS84_plus_horizontal_as_in_LAS) { + + const auto wkt = "COMPD_CS[\"WGS 84 + WGS 84\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]]]"; + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj); + ASSERT_TRUE(crs != nullptr); + auto baseCRS = nn_dynamic_pointer_cast<GeographicCRS>(crs->baseCRS()); + ASSERT_TRUE(baseCRS != nullptr); + EXPECT_EQ(baseCRS->nameStr(), "WGS 84"); + EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U); + + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, + COMPD_CS_horizontal_bound_geog_plus_vertical_ellipsoidal_height) { + // See https://github.com/OSGeo/PROJ/issues/2228 + const char *wkt = + "COMPD_CS[\"NAD83 + Ellipsoid (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (Meters)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast<BoundCRS>(obj); + ASSERT_TRUE(crs != nullptr); + auto baseCRS = nn_dynamic_pointer_cast<GeographicCRS>(crs->baseCRS()); + ASSERT_TRUE(baseCRS != nullptr); + EXPECT_EQ(baseCRS->nameStr(), "NAD83"); + EXPECT_EQ(baseCRS->coordinateSystem()->axisList().size(), 3U); + + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, + COMPD_CS_horizontal_projected_plus_vertical_ellipsoidal_height) { + // Variant of above + const char *wkt = + "COMPD_CS[\"WGS 84 / UTM zone 31N + Ellipsoid (Meters)\",\n" + " PROJCS[\"WGS 84 / UTM zone 31N\",\n" + " GEOGCS[\"WGS 84\",\n" + " DATUM[\"WGS_1984\",\n" + " SPHEROID[\"WGS 84\",6378137,298.257223563,\n" + " AUTHORITY[\"EPSG\",\"7030\"]],\n" + " AUTHORITY[\"EPSG\",\"6326\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4326\"]],\n" + " PROJECTION[\"Transverse_Mercator\"],\n" + " PARAMETER[\"latitude_of_origin\",0],\n" + " PARAMETER[\"central_meridian\",3],\n" + " PARAMETER[\"scale_factor\",0.9996],\n" + " PARAMETER[\"false_easting\",500000],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Easting\",EAST],\n" + " AXIS[\"Northing\",NORTH],\n" + " AUTHORITY[\"EPSG\",\"32631\"]],\n" + " VERT_CS[\"Ellipsoid (Meters)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "WGS 84 / UTM zone 31N"); + EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U); + + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); +} + +// --------------------------------------------------------------------------- + +TEST(wkt_parse, + COMPD_CS_horizontal_geog_plus_vertical_ellipsoidal_height_non_metre) { + // See https://github.com/OSGeo/PROJ/issues/2232 + const char *wkt = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto crs = nn_dynamic_pointer_cast<GeographicCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), "NAD83 (Ellipsoid (US Feet))"); + EXPECT_EQ(crs->coordinateSystem()->axisList().size(), 3U); + EXPECT_NEAR(crs->coordinateSystem()->axisList()[2]->unit().conversionToSI(), + 0.304800609601219, 1e-15); + + EXPECT_EQ( + crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL, dbContext) + .get()), + wkt); +} + +// --------------------------------------------------------------------------- + TEST(wkt_parse, COORDINATEOPERATION) { std::string src_wkt; diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 92b5f923..51bbf50b 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6379,6 +6379,66 @@ TEST(operation, geogCRS_to_boundCRS_of_geogCRS) { // --------------------------------------------------------------------------- +TEST(operation, boundCRS_to_geogCRS_same_datum_context) { + auto boundCRS = BoundCRS::createFromTOWGS84( + GeographicCRS::EPSG_4269, std::vector<double>{1, 2, 3, 4, 5, 6, 7}); + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + boundCRS, GeographicCRS::EPSG_4269, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=noop"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, boundCRS_to_geogCRS_hubCRS_and_targetCRS_same_but_baseCRS_not) { + const char *wkt = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + + auto dbContext = DatabaseContext::create(); + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto boundCRS = nn_dynamic_pointer_cast<BoundCRS>(obj); + ASSERT_TRUE(boundCRS != nullptr); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(boundCRS), GeographicCRS::EPSG_4979, ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=unitconvert +z_in=us-ft +z_out=m"); +} + +// --------------------------------------------------------------------------- + TEST(operation, boundCRS_to_boundCRS) { auto utm31 = ProjectedCRS::create( PropertyMap(), GeographicCRS::EPSG_4807, @@ -7195,6 +7255,237 @@ TEST(operation, // --------------------------------------------------------------------------- +TEST(operation, + compoundCRS_with_boundVerticalCRS_from_grids_to_geogCRS_with_ftus_ctxt) { + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 + NAVD88 height - Geoid12B (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"@foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast<CompoundCRS>(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + const char *wktDst = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto objDst = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst); + auto dstCRS = nn_dynamic_pointer_cast<GeographicCRS>(objDst); + ASSERT_TRUE(dstCRS != nullptr); + + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +z_in=m " + "+xy_out=deg +z_out=us-ft " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + +TEST( + operation, + compoundCRS_with_boundGeogCRS_boundVerticalCRS_from_grids_to_boundGeogCRS_with_ftus_ctxt) { + + // Variant of above but with TOWGS84 in source & target CRS + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 + NAVD88 height - Geoid12B (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"@foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast<CompoundCRS>(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + const char *wktDst = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto objDst = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst); + auto dstCRS = nn_dynamic_pointer_cast<BoundCRS>(objDst); + ASSERT_TRUE(dstCRS != nullptr); + + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_EQ(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=rad +z_in=m " + "+xy_out=deg +z_out=us-ft"); +} + +// --------------------------------------------------------------------------- + +TEST( + operation, + compoundCRS_with_boundVerticalCRS_from_grids_to_boundGeogCRS_with_ftus_ctxt) { + + // Variant of above but with TOWGS84 in target CRS only + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 + NAVD88 height - Geoid12B (Meters)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"@foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast<CompoundCRS>(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + const char *wktDst = + "COMPD_CS[\"NAD83 + Ellipsoid (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"Ellipsoid (US Feet)\",\n" + " VERT_DATUM[\"Ellipsoid\",2002],\n" + " UNIT[\"US survey foot\",0.304800609601219,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"; + auto objDst = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktDst); + auto dstCRS = nn_dynamic_pointer_cast<BoundCRS>(objDst); + ASSERT_TRUE(dstCRS != nullptr); + + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + NN_NO_CHECK(srcCRS), NN_NO_CHECK(dstCRS), ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step +proj=vgridshift +grids=@foo.gtx +multiplier=1 " + "+step +proj=unitconvert +xy_in=rad +z_in=m " + "+xy_out=deg +z_out=us-ft " + "+step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST(operation, geocent_to_compoundCRS) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=geocent +datum=WGS84 +units=m +type=crs"); @@ -7545,7 +7836,7 @@ TEST( CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); auto list = CoordinateOperationFactory::create()->createOperations( NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst), ctxt); - ASSERT_GE(list.size(), 1U); + ASSERT_EQ(list.size(), 1U); EXPECT_EQ(list[0]->nameStr(), "Inverse of unnamed + " "Transformation from NAD83 to WGS84 + " @@ -7566,6 +7857,85 @@ TEST( // --------------------------------------------------------------------------- +TEST(operation, compoundCRS_to_compoundCRS_issue_2232) { + auto objSrc = WKTParser().createFromWKT( + "COMPD_CS[\"NAD83 / Alabama West + NAVD88 height - Geoid12B " + "(Meters)\",\n" + " PROJCS[\"NAD83 / Alabama West\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " PROJECTION[\"Transverse_Mercator\"],\n" + " PARAMETER[\"latitude_of_origin\",30],\n" + " PARAMETER[\"central_meridian\",-87.5],\n" + " PARAMETER[\"scale_factor\",0.999933333],\n" + " PARAMETER[\"false_easting\",600000],\n" + " PARAMETER[\"false_northing\",0],\n" + " UNIT[\"metre\",1,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"X\",EAST],\n" + " AXIS[\"Y\",NORTH],\n" + " AUTHORITY[\"EPSG\",\"26930\"]],\n" + " VERT_CS[\"NAVD88 height - Geoid12B (Meters)\",\n" + " VERT_DATUM[\"North American Vertical Datum 1988\",2005,\n" + " EXTENSION[\"PROJ4_GRIDS\",\"foo.gtx\"],\n" + " AUTHORITY[\"EPSG\",\"5103\"]],\n" + " UNIT[\"metre\",1.0,\n" + " AUTHORITY[\"EPSG\",\"9001\"]],\n" + " AXIS[\"Gravity-related height\",UP],\n" + " AUTHORITY[\"EPSG\",\"5703\"]]]"); + auto src = nn_dynamic_pointer_cast<CRS>(objSrc); + ASSERT_TRUE(src != nullptr); + + auto objDst = WKTParser().createFromWKT( + "COMPD_CS[\"NAD83 + some CRS (US Feet)\",\n" + " GEOGCS[\"NAD83\",\n" + " DATUM[\"North_American_Datum_1983\",\n" + " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n" + " AUTHORITY[\"EPSG\",\"7019\"]],\n" + " TOWGS84[0,0,0,0,0,0,0],\n" + " AUTHORITY[\"EPSG\",\"6269\"]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " AUTHORITY[\"EPSG\",\"8901\"]],\n" + " UNIT[\"degree\",0.0174532925199433,\n" + " AUTHORITY[\"EPSG\",\"9122\"]],\n" + " AUTHORITY[\"EPSG\",\"4269\"]],\n" + " VERT_CS[\"some CRS (US Feet)\",\n" + " VERT_DATUM[\"some datum\",2005],\n" + " UNIT[\"US survey foot\",0.3048006096012192,\n" + " AUTHORITY[\"EPSG\",\"9003\"]],\n" + " AXIS[\"Up\",UP]]]"); + auto dst = nn_dynamic_pointer_cast<CRS>(objDst); + ASSERT_TRUE(dst != nullptr); + + auto dbContext = DatabaseContext::create(); + auto authFactory = AuthorityFactory::create(dbContext, "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setGridAvailabilityUse( + CoordinateOperationContext::GridAvailabilityUse:: + IGNORE_GRID_AVAILABILITY); + ctxt->setSpatialCriterion( + CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION); + + auto list = CoordinateOperationFactory::create()->createOperations( + NN_CHECK_ASSERT(src), NN_CHECK_ASSERT(dst), ctxt); + EXPECT_GE(list.size(), 1U); + + auto list2 = CoordinateOperationFactory::create()->createOperations( + NN_CHECK_ASSERT(dst), NN_CHECK_ASSERT(src), ctxt); + EXPECT_EQ(list2.size(), list.size()); +} + +// --------------------------------------------------------------------------- + TEST(operation, compoundCRS_to_compoundCRS_context) { auto authFactory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); |