diff options
| -rw-r--r-- | data/sql/grid_alternatives.sql | 16 | ||||
| -rw-r--r-- | scripts/reference_exported_symbols.txt | 1 | ||||
| -rw-r--r-- | src/apply_gridshift.cpp | 6 | ||||
| -rw-r--r-- | src/apps/proj.cpp | 8 | ||||
| -rw-r--r-- | src/gridcatalog.cpp | 4 | ||||
| -rw-r--r-- | src/iso19111/c_api.cpp | 28 | ||||
| -rw-r--r-- | src/iso19111/crs.cpp | 49 | ||||
| -rw-r--r-- | src/iso19111/io.cpp | 10 | ||||
| -rw-r--r-- | src/nad_cvt.cpp | 43 | ||||
| -rw-r--r-- | src/proj.h | 3 | ||||
| -rw-r--r-- | src/proj_internal.h | 4 | ||||
| -rw-r--r-- | src/transform.cpp | 16 | ||||
| -rw-r--r-- | test/cli/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | test/cli/Makefile.am | 7 | ||||
| -rw-r--r-- | test/cli/ntv2_out.dist | 3 | ||||
| -rwxr-xr-x | test/cli/testntv2 | 8 | ||||
| -rwxr-xr-x | test/cli/testproj | 55 | ||||
| -rw-r--r-- | test/cli/testproj_out.dist | 1 | ||||
| -rw-r--r-- | test/unit/pj_transform_test.cpp | 34 | ||||
| -rw-r--r-- | test/unit/test_c_api.cpp | 30 | ||||
| -rw-r--r-- | test/unit/test_crs.cpp | 40 | ||||
| -rw-r--r-- | test/unit/test_io.cpp | 35 |
22 files changed, 376 insertions, 26 deletions
diff --git a/data/sql/grid_alternatives.sql b/data/sql/grid_alternatives.sql index 6a893c5e..070a0a03 100644 --- a/data/sql/grid_alternatives.sql +++ b/data/sql/grid_alternatives.sql @@ -1106,6 +1106,22 @@ INSERT INTO grid_alternatives(original_grid_name, 'proj-datumgrid-oceania', NULL, NULL, NULL, NULL); +INSERT INTO grid_alternatives(original_grid_name, + proj_grid_name, + proj_grid_format, + proj_method, + inverse_direction, + package_name, + url, direct_download, open_license, directory) + SELECT grid_name, + 'AUSGeoid98.gtx', + 'GTX', + 'vgridshift', + 1, + 'proj-datumgrid-oceania', + NULL, NULL, NULL, NULL FROM grid_transformation WHERE + grid_name LIKE '%DAT.htm' AND name LIKE 'GDA94 to AHD height%'; + -- Netherlands / RDNAP (non-free grids) INSERT INTO grid_alternatives(original_grid_name, diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt index e20f29b3..3f4cd51b 100644 --- a/scripts/reference_exported_symbols.txt +++ b/scripts/reference_exported_symbols.txt @@ -801,6 +801,7 @@ proj_context_use_proj4_init_rules proj_convert_conversion_to_other_method proj_coord proj_coord_error() +proj_coordoperation_create_inverse proj_coordoperation_get_accuracy proj_coordoperation_get_grid_used proj_coordoperation_get_grid_used_count diff --git a/src/apply_gridshift.cpp b/src/apply_gridshift.cpp index e6c7b2b1..39e7c3b5 100644 --- a/src/apply_gridshift.cpp +++ b/src/apply_gridshift.cpp @@ -115,7 +115,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse, /* Determine which grid is the correct given an input coordinate. */ /************************************************************************/ -static struct CTABLE* find_ctable(projCtx ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables) { +struct CTABLE* find_ctable(projCtx ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables) { int itable; /* keep trying till we find a table that works */ @@ -210,7 +210,7 @@ int pj_apply_gridshift_3( projCtx ctx, PJ_GRIDINFO **gridlist, int gridlist_coun ct = find_ctable(ctx, input, gridlist_count, gridlist); if( ct != nullptr ) { - output = nad_cvt( input, inverse, ct ); + output = nad_cvt( ctx, input, inverse, ct, gridlist_count, gridlist); if ( output.lam != HUGE_VAL && debug_count++ < 20 ) pj_log( ctx, PJ_LOG_DEBUG_MINOR, "pj_apply_gridshift(): used %s", ct->id ); @@ -356,7 +356,7 @@ PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction) { } inverse = direction == PJ_FWD ? 0 : 1; - out = nad_cvt(lp, inverse, ct); + out = nad_cvt(P->ctx, lp, inverse, ct, P->gridlist_count, P->gridlist); if (out.lam == HUGE_VAL || out.phi == HUGE_VAL) pj_ctx_set_errno(P->ctx, PJD_ERR_GRID_AREA); diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index 7fe08023..09c8a81d 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -490,6 +490,14 @@ int main(int argc, char **argv) { exit(0); } + // Ugly hack. See https://github.com/OSGeo/PROJ/issues/1782 + if( Proj->right == PJ_IO_UNITS_WHATEVER && Proj->descr && + strncmp(Proj->descr, "General Oblique Transformation", + strlen("General Oblique Transformation")) == 0 ) + { + Proj->right = PJ_IO_UNITS_PROJECTED; + } + if (inverse) { if (!Proj->inv) emess(3,"inverse projection not available"); diff --git a/src/gridcatalog.cpp b/src/gridcatalog.cpp index 15d81dd7..9b94fef8 100644 --- a/src/gridcatalog.cpp +++ b/src/gridcatalog.cpp @@ -164,7 +164,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse, return PJD_ERR_FAILED_TO_LOAD_GRID; } - output_after = nad_cvt( input, inverse, gi->ct ); + output_after = nad_cvt( defn->ctx, input, inverse, gi->ct, 0, nullptr ); if( output_after.lam == HUGE_VAL ) { if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) @@ -213,7 +213,7 @@ int pj_gc_apply_gridshift( PJ *defn, int inverse, return PJD_ERR_FAILED_TO_LOAD_GRID; } - output_before = nad_cvt( input, inverse, gi->ct ); + output_before = nad_cvt( defn->ctx, input, inverse, gi->ct, 0, nullptr ); if( output_before.lam == HUGE_VAL ) { if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 9db9e5b3..5e2ac522 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -7769,6 +7769,34 @@ PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) { // --------------------------------------------------------------------------- +/** \brief Returns a PJ* coordinate operation object which represents the + * inverse operation of the specified coordinate operation. + * + * @param ctx PROJ context, or NULL for default context + * @param obj Object of type CoordinateOperation (must not be NULL) + * @return a new PJ* object to free with proj_destroy() in case of success, or + * nullptr in case of error + * @since 6.3 + */ +PJ *proj_coordoperation_create_inverse(PJ_CONTEXT *ctx, const PJ *obj) { + + SANITIZE_CTX(ctx); + auto co = dynamic_cast<const CoordinateOperation *>(obj->iso_obj.get()); + if (!co) { + proj_log_error(ctx, __FUNCTION__, + "Object is not a CoordinateOperation"); + return nullptr; + } + try { + return pj_obj_create(ctx, co->inverse()); + } catch (const std::exception &e) { + proj_log_debug(ctx, __FUNCTION__, e.what()); + return nullptr; + } +} + +// --------------------------------------------------------------------------- + /** \brief Returns the number of steps of a concatenated operation. * * The input object must be a concatenated operation. diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 49cc050f..46aeee20 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -49,6 +49,7 @@ #include <algorithm> #include <cassert> +#include <cmath> #include <cstring> #include <iostream> #include <memory> @@ -595,12 +596,44 @@ CRSNNPtr CRS::alterId(const std::string &authName, //! @cond Doxygen_Suppress -static bool isAxisListNorthEast( +static bool mustAxisOrderBeSwitchedForVisualizationInternal( const std::vector<cs::CoordinateSystemAxisNNPtr> &axisList) { const auto &dir0 = axisList[0]->direction(); const auto &dir1 = axisList[1]->direction(); - return (&dir0 == &cs::AxisDirection::NORTH && - &dir1 == &cs::AxisDirection::EAST); + if (&dir0 == &cs::AxisDirection::NORTH && + &dir1 == &cs::AxisDirection::EAST) { + return true; + } + + // Address EPSG:32661 "WGS 84 / UPS North (N,E)" + if (&dir0 == &cs::AxisDirection::SOUTH && + &dir1 == &cs::AxisDirection::SOUTH) { + const auto &meridian0 = axisList[0]->meridian(); + const auto &meridian1 = axisList[1]->meridian(); + return meridian0 != nullptr && meridian1 != nullptr && + std::abs(meridian0->longitude().convertToUnit( + common::UnitOfMeasure::DEGREE) - + 180.0) < 1e-10 && + std::abs(meridian1->longitude().convertToUnit( + common::UnitOfMeasure::DEGREE) - + 90.0) < 1e-10; + } + + // Address EPSG:32761 "WGS 84 / UPS South (N,E)" + if (&dir0 == &cs::AxisDirection::NORTH && + &dir1 == &cs::AxisDirection::NORTH) { + const auto &meridian0 = axisList[0]->meridian(); + const auto &meridian1 = axisList[1]->meridian(); + return meridian0 != nullptr && meridian1 != nullptr && + std::abs(meridian0->longitude().convertToUnit( + common::UnitOfMeasure::DEGREE) - + 0.0) < 1e-10 && + std::abs(meridian1->longitude().convertToUnit( + common::UnitOfMeasure::DEGREE) - + 90.0) < 1e-10; + } + + return false; } // --------------------------------------------------------------------------- @@ -616,12 +649,14 @@ bool CRS::mustAxisOrderBeSwitchedForVisualization() const { const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this); if (geogCRS) { - return isAxisListNorthEast(geogCRS->coordinateSystem()->axisList()); + return mustAxisOrderBeSwitchedForVisualizationInternal( + geogCRS->coordinateSystem()->axisList()); } const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this); if (projCRS) { - return isAxisListNorthEast(projCRS->coordinateSystem()->axisList()); + return mustAxisOrderBeSwitchedForVisualizationInternal( + projCRS->coordinateSystem()->axisList()); } return false; @@ -655,7 +690,7 @@ CRSNNPtr CRS::normalizeForVisualization() const { const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this); if (geogCRS) { const auto &axisList = geogCRS->coordinateSystem()->axisList(); - if (isAxisListNorthEast(axisList)) { + if (mustAxisOrderBeSwitchedForVisualizationInternal(axisList)) { auto cs = axisList.size() == 2 ? cs::EllipsoidalCS::create(util::PropertyMap(), axisList[1], axisList[0]) @@ -670,7 +705,7 @@ CRSNNPtr CRS::normalizeForVisualization() const { const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this); if (projCRS) { const auto &axisList = projCRS->coordinateSystem()->axisList(); - if (isAxisListNorthEast(axisList)) { + if (mustAxisOrderBeSwitchedForVisualizationInternal(axisList)) { auto cs = axisList.size() == 2 ? cs::CartesianCS::create(util::PropertyMap(), axisList[1], diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index e0e6152a..c704e1c1 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -2683,7 +2683,9 @@ WKTParser::Private::buildGeodeticCRS(const WKTNodeNNPtr &node) { auto cs = buildCS(csNode, node, angularUnit); auto ellipsoidalCS = nn_dynamic_pointer_cast<EllipsoidalCS>(cs); if (ellipsoidalCS) { - assert(!ci_equal(nodeName, WKTConstants::GEOCCS)); + if (ci_equal(nodeName, WKTConstants::GEOCCS)) { + throw ParsingException("ellipsoidal CS not expected in GEOCCS"); + } try { auto crs = GeographicCRS::create(props, datum, datumEnsemble, NN_NO_CHECK(ellipsoidalCS)); @@ -7154,6 +7156,12 @@ void PROJStringFormatter::stopInversion() { // the current end of steps for (auto iter = startIter; iter != d->steps_.end(); ++iter) { iter->inverted = !iter->inverted; + for (auto ¶mValue : iter->paramValues) { + if (paramValue.key == "omit_fwd") + paramValue.key = "omit_inv"; + else if (paramValue.key == "omit_inv") + paramValue.key = "omit_fwd"; + } } // And reverse the order of steps in that range as well. std::reverse(startIter, d->steps_.end()); diff --git a/src/nad_cvt.cpp b/src/nad_cvt.cpp index 79441d0a..e8b8e9b7 100644 --- a/src/nad_cvt.cpp +++ b/src/nad_cvt.cpp @@ -7,10 +7,12 @@ #include "proj_internal.h" #include <math.h> +#include <limits> + #define MAX_ITERATIONS 10 #define TOL 1e-12 -PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) { +PJ_LP nad_cvt(projCtx ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables) { PJ_LP t, tb,del, dif; int i = MAX_ITERATIONS; const double toltol = TOL*TOL; @@ -40,18 +42,37 @@ PJ_LP nad_cvt(PJ_LP in, int inverse, struct CTABLE *ct) { do { del = nad_intr(t, ct); - /* This case used to return failure, but I have - changed it to return the first order approximation - of the inverse shift. This avoids cases where the - grid shift *into* this grid came from another grid. - While we aren't returning optimally correct results - I feel a close result in this case is better than - no result. NFW - To demonstrate use -112.5839956 49.4914451 against - the NTv2 grid shift file from Canada. */ - if (del.lam == HUGE_VAL) + /* In case we are (or have switched) on the null grid, exit now */ + if( del.lam == 0 && del.phi == 0 ) break; + /* We can possibly go outside of the initial guessed grid, so try */ + /* to fetch a new grid into which iterate... */ + if (del.lam == HUGE_VAL) + { + if( grid_count == 0 ) + break; + PJ_LP lp; + lp.lam = t.lam + ct->ll.lam; + lp.phi = t.phi + ct->ll.phi; + auto newCt = find_ctable(ctx, lp, grid_count, tables); + if( newCt == nullptr || newCt == ct ) + break; + pj_log(ctx, PJ_LOG_DEBUG_MINOR, "Switching from grid %s to grid %s", + ct->id, + newCt->id); + ct = newCt; + t.lam = lp.lam - ct->ll.lam; + t.phi = lp.phi - ct->ll.phi; + tb = in; + tb.lam -= ct->ll.lam; + tb.phi -= ct->ll.phi; + tb.lam = adjlon (tb.lam - M_PI) + M_PI; + dif.lam = std::numeric_limits<double>::max(); + dif.phi = std::numeric_limits<double>::max(); + continue; + } + dif.lam = t.lam - del.lam - tb.lam; dif.phi = t.phi + del.phi - tb.phi; t.lam -= dif.lam; @@ -1091,6 +1091,9 @@ int PROJ_DLL proj_coordoperation_get_towgs84_values(PJ_CONTEXT *ctx, int value_count, int emit_error_if_incompatible); +PJ PROJ_DLL *proj_coordoperation_create_inverse(PJ_CONTEXT *ctx, const PJ *obj); + + int PROJ_DLL proj_concatoperation_get_step_count(PJ_CONTEXT *ctx, const PJ *concatoperation); diff --git a/src/proj_internal.h b/src/proj_internal.h index 3ca927a3..3e219682 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -863,8 +863,10 @@ int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *); int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *); /* nadcon related protos */ +struct CTABLE* find_ctable(projCtx_t *ctx, PJ_LP input, int grid_count, PJ_GRIDINFO **tables); + PJ_LP nad_intr(PJ_LP, struct CTABLE *); -PJ_LP nad_cvt(PJ_LP, int, struct CTABLE *); +PJ_LP nad_cvt(projCtx_t *ctx, PJ_LP in, int inverse, struct CTABLE *ct, int grid_count, PJ_GRIDINFO **tables); struct CTABLE *nad_init(projCtx_t *ctx, char *); struct CTABLE *nad_ctable_init( projCtx_t *ctx, struct projFileAPI_t* fid ); int nad_ctable_load( projCtx_t *ctx, struct CTABLE *, struct projFileAPI_t* fid ); diff --git a/src/transform.cpp b/src/transform.cpp index d111d835..781c0061 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -231,6 +231,14 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double * return 0; } + // Ugly hack. See https://github.com/OSGeo/PROJ/issues/1782 + if( P->right == PJ_IO_UNITS_WHATEVER && P->descr && + strncmp(P->descr, "General Oblique Transformation", + strlen("General Oblique Transformation")) == 0 ) + { + P->right = PJ_IO_UNITS_PROJECTED; + } + for( i = 0; i <n; i++ ) { PJ_XY projected_loc; @@ -342,6 +350,14 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double * return 0; } + // Ugly hack. See https://github.com/OSGeo/PROJ/issues/1782 + if( P->right == PJ_IO_UNITS_WHATEVER && P->descr && + strncmp(P->descr, "General Oblique Transformation", + strlen("General Oblique Transformation")) == 0 ) + { + P->right = PJ_IO_UNITS_PROJECTED; + } + /* Fallback to the original PROJ.4 API 2d inversion - inv */ for( i = 0; i < n; i++ ) { PJ_XY projected_loc; diff --git a/test/cli/CMakeLists.txt b/test/cli/CMakeLists.txt index 4e1ab75a..feef5bf0 100644 --- a/test/cli/CMakeLists.txt +++ b/test/cli/CMakeLists.txt @@ -7,6 +7,7 @@ set(PROJINFO_BIN "projinfo") set(CCT_BIN "cct") proj_add_test_script_sh("test27" PROJ_BIN) proj_add_test_script_sh("test83" PROJ_BIN) +proj_add_test_script_sh("testproj" PROJ_BIN) proj_add_test_script_sh("testvarious" CS2CS_BIN) proj_add_test_script_sh("testdatumfile" CS2CS_BIN "connu") proj_add_test_script_sh("testIGNF" CS2CS_BIN "ntf_r93.gsb") diff --git a/test/cli/Makefile.am b/test/cli/Makefile.am index 758352c6..fb7169e4 100644 --- a/test/cli/Makefile.am +++ b/test/cli/Makefile.am @@ -10,6 +10,7 @@ CCTEXE = $(EXEPATH)/cct # PROJ.4 test scripts TEST27 = $(THIS_DIR)/test27 TEST83 = $(THIS_DIR)/test83 +TESTPROJ = $(THIS_DIR)/testproj TESTNTV2 = $(THIS_DIR)/testntv2 TESTVARIOUS = $(THIS_DIR)/testvarious TESTFLAKY = $(THIS_DIR)/testflaky @@ -24,6 +25,7 @@ EXTRA_DIST = pj_out27.dist pj_out83.dist td_out.dist \ testIGNF proj_outIGNF.dist \ testprojinfo testprojinfo_out.dist \ testcct testcct_out.dist \ + testproj testproj_out.dist \ CMakeLists.txt testprojinfo-check: @@ -35,6 +37,9 @@ test27-check: test83-check: $(TEST83) $(PROJEXE) +testproj-check: + $(TESTPROJ) $(PROJEXE) + testvarious-check: PROJ_LIB=$(PROJ_LIB) $(TESTVARIOUS) $(CS2CSEXE) @@ -56,4 +61,4 @@ testntv2-check: testcct-check: PROJ_LIB=$(PROJ_LIB) $(TESTCCT) $(CCTEXE) -check-local: testprojinfo-check test27-check test83-check testvarious-check testdatumfile-check testign-check testntv2-check testcct-check +check-local: testprojinfo-check test27-check test83-check testproj-check testvarious-check testdatumfile-check testign-check testntv2-check testcct-check diff --git a/test/cli/ntv2_out.dist b/test/cli/ntv2_out.dist index 9a97f9cf..650a69d8 100644 --- a/test/cli/ntv2_out.dist +++ b/test/cli/ntv2_out.dist @@ -9,3 +9,6 @@ Try with NTv2 and NTv1 together ... falls back to NTv1 99d00'00.000"W 65d00'00.000"N 0.0 99d0'1.5885"W 65d0'1.3482"N 0.000 111d00'00.000"W 46d00'00.000"N 0.0 111d0'3.1897"W 45d59'59.7489"N 0.000 111d00'00.000"W 47d30'00.000"N 0.0 111d0'2.7989"W 47d29'59.9896"N 0.000 +############################################################## +Switching between NTv2 subgrids +-112.5839956 49.4914451 0 -112.58307487 49.49145197 0.00000000 diff --git a/test/cli/testntv2 b/test/cli/testntv2 index 73371dbe..2a31304e 100755 --- a/test/cli/testntv2 +++ b/test/cli/testntv2 @@ -52,6 +52,14 @@ $EXE +proj=latlong +ellps=clrk66 +nadgrids=ntv2_0.gsb,ntv1_can.dat,conus \ 111d00'00.000"W 46d00'00.000"N 0.0 111d00'00.000"W 47d30'00.000"N 0.0 EOF + +echo "##############################################################" >> ${OUT} +echo Switching between NTv2 subgrids >> ${OUT} +# Initial guess is in ALraymnd, going to parent CAwest afterwards +$EXE +proj=latlong +datum=NAD83 +to +proj=latlong +ellps=clrk66 +nadgrids=ntv2_0.gsb -E -d 8 >>${OUT} <<EOF +-112.5839956 49.4914451 0 +EOF + # ############################################################################## # Done! diff --git a/test/cli/testproj b/test/cli/testproj new file mode 100755 index 00000000..8686224e --- /dev/null +++ b/test/cli/testproj @@ -0,0 +1,55 @@ +: +# Script to test proj exe +# +TEST_CLI_DIR=`dirname $0` +EXE=$1 + +usage() +{ + echo "Usage: ${0} <path to 'proj' program>" + echo + exit 1 +} + +if test -z "${EXE}"; then + EXE=../../src/cs2cs +fi + +if test ! -x ${EXE}; then + echo "*** ERROR: Can not find '${EXE}' program!" + exit 1 +fi + +if test -z "${PROJ_LIB}"; then + export PROJ_LIB="`dirname $0`/../../data" +fi + +echo "============================================" +echo "Running ${0} using ${EXE}:" +echo "============================================" + +OUT=testproj_out +# +echo "doing tests into file ${OUT}, please wait" +# +$EXE +ellps=WGS84 +proj=ob_tran +o_proj=latlon +o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 +to_meter=0.0174532925199433 +no_defs -E -f '%.3f' >${OUT} <<EOF +2 49 +EOF + +# +# do 'diff' with distribution results +echo "diff ${OUT} with testproj_out.dist" +diff -u -b ${OUT} ${TEST_CLI_DIR}/testproj_out.dist +if [ $? -ne 0 ] ; then + echo "" + echo "PROBLEMS HAVE OCCURRED" + echo "test file ${OUT} saved" + echo + exit 100 +else + echo "TEST OK" + echo "test file ${OUT} removed" + echo + /bin/rm -f ${OUT} + exit 0 +fi diff --git a/test/cli/testproj_out.dist b/test/cli/testproj_out.dist new file mode 100644 index 00000000..fad0a014 --- /dev/null +++ b/test/cli/testproj_out.dist @@ -0,0 +1 @@ +2 49 2.000 49.000 diff --git a/test/unit/pj_transform_test.cpp b/test/unit/pj_transform_test.cpp index b3a061b4..1f4473c1 100644 --- a/test/unit/pj_transform_test.cpp +++ b/test/unit/pj_transform_test.cpp @@ -614,4 +614,38 @@ TEST(proj_api_h, pj_set_finder) { pj_set_finder(nullptr); } +// --------------------------------------------------------------------------- + +TEST(pj_transform_test, ob_tran_to_meter_as_dest) { + auto src = pj_init_plus( + "+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs"); + auto dst = pj_init_plus("+ellps=WGS84 +proj=ob_tran +o_proj=latlon " + "+o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 " + "+to_meter=0.0174532925199433 +no_defs"); + double x = 2 * DEG_TO_RAD; + double y = 49 * DEG_TO_RAD; + EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); + EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; + EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; + pj_free(src); + pj_free(dst); +} + +// --------------------------------------------------------------------------- + +TEST(pj_transform_test, ob_tran_to_meter_as_srouce) { + auto src = pj_init_plus("+ellps=WGS84 +proj=ob_tran +o_proj=latlon " + "+o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 " + "+to_meter=0.0174532925199433 +no_defs"); + auto dst = pj_init_plus( + "+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs"); + double x = 2 * DEG_TO_RAD; + double y = 49 * DEG_TO_RAD; + EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); + EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; + EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; + pj_free(src); + pj_free(dst); +} + } // namespace diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index f87f6589..22e2ac11 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -3486,6 +3486,36 @@ TEST_F(CApi, proj_normalize_for_visualization_on_crs) { // --------------------------------------------------------------------------- +TEST_F(CApi, proj_coordoperation_create_inverse) { + + auto P = proj_create( + m_ctxt, "+proj=pipeline +step +proj=axisswap +order=2,1 +step " + "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push " + "+v_3 +step +proj=cart +ellps=evrst30 +step +proj=helmert " + "+x=293 +y=836 +z=318 +rx=0.5 +ry=1.6 +rz=-2.8 +s=2.1 " + "+convention=position_vector +step +inv +proj=cart " + "+ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert " + "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"); + ObjectKeeper keeper_P(P); + ASSERT_NE(P, nullptr); + auto Pinversed = proj_coordoperation_create_inverse(m_ctxt, P); + ObjectKeeper keeper_Pinversed(Pinversed); + ASSERT_NE(Pinversed, nullptr); + + auto projstr = proj_as_proj_string(m_ctxt, Pinversed, PJ_PROJ_5, nullptr); + ASSERT_NE(projstr, nullptr); + EXPECT_EQ(std::string(projstr), + "+proj=pipeline +step +proj=axisswap +order=2,1 +step " + "+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 " + "+step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=293 " + "+y=836 +z=318 +rx=0.5 +ry=1.6 +rz=-2.8 +s=2.1 " + "+convention=position_vector +step +inv +proj=cart " + "+ellps=evrst30 +step +proj=pop +v_3 +step +proj=unitconvert " + "+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"); +} + +// --------------------------------------------------------------------------- + TEST_F(CApi, proj_get_remarks) { auto co = proj_create_from_database(m_ctxt, "EPSG", "8048", PJ_CATEGORY_COORDINATE_OPERATION, false, diff --git a/test/unit/test_crs.cpp b/test/unit/test_crs.cpp index e896853f..33d67e0a 100644 --- a/test/unit/test_crs.cpp +++ b/test/unit/test_crs.cpp @@ -943,6 +943,16 @@ TEST(crs, EPSG_32661_projected_north_pole_north_east) { ASSERT_TRUE(op != nullptr); EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), proj_string); + + auto opNormalized = op->normalizeForVisualization(); + auto proj_string_normalized = + "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step " + "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 " + "+ellps=WGS84"; + EXPECT_EQ( + opNormalized->exportToPROJString(PROJStringFormatter::create().get()), + proj_string_normalized); } // --------------------------------------------------------------------------- @@ -964,6 +974,16 @@ TEST(crs, EPSG_5041_projected_north_pole_east_north) { ASSERT_TRUE(op != nullptr); EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), proj_string); + + auto opNormalized = op->normalizeForVisualization(); + auto proj_string_normalized = + "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step " + "+proj=stere +lat_0=90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 " + "+ellps=WGS84"; + EXPECT_EQ( + opNormalized->exportToPROJString(PROJStringFormatter::create().get()), + proj_string_normalized); } // --------------------------------------------------------------------------- @@ -985,6 +1005,16 @@ TEST(crs, EPSG_32761_projected_south_pole_north_east) { ASSERT_TRUE(op != nullptr); EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), proj_string); + + auto opNormalized = op->normalizeForVisualization(); + auto proj_string_normalized = + "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step " + "+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 " + "+ellps=WGS84"; + EXPECT_EQ( + opNormalized->exportToPROJString(PROJStringFormatter::create().get()), + proj_string_normalized); } // --------------------------------------------------------------------------- @@ -1006,6 +1036,16 @@ TEST(crs, EPSG_5042_projected_south_pole_east_north) { ASSERT_TRUE(op != nullptr); EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), proj_string); + + auto opNormalized = op->normalizeForVisualization(); + auto proj_string_normalized = + "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad " + "+step " + "+proj=stere +lat_0=-90 +lon_0=0 +k=0.994 +x_0=2000000 +y_0=2000000 " + "+ellps=WGS84"; + EXPECT_EQ( + opNormalized->exportToPROJString(PROJStringFormatter::create().get()), + proj_string_normalized); } // --------------------------------------------------------------------------- diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 38dfc2b4..15ab8706 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -5577,6 +5577,17 @@ TEST(wkt_parse, invalid_GEOCCS) { "NORTH],AXIS[\"longitude\",EAST]]"), ParsingException); + // ellipsoidal CS is invalid in a GEOCCS + EXPECT_THROW(WKTParser().createFromWKT( + "GEOCCS[\"WGS 84\",DATUM[\"World Geodetic System 1984\"," + "ELLIPSOID[\"WGS 84\",6378274,298.257223564," + "LENGTHUNIT[\"metre\",1]]]," + "CS[ellipsoidal,2],AXIS[\"geodetic latitude (Lat)\",north," + "ANGLEUNIT[\"degree\",0.0174532925199433]]," + "AXIS[\"geodetic longitude (Lon)\",east," + "ANGLEUNIT[\"degree\",0.0174532925199433]]]"), + ParsingException); + // 3 axis required EXPECT_THROW(WKTParser().createFromWKT( "GEOCCS[\"x\",DATUM[\"x\",SPHEROID[\"x\",1,0.5]],PRIMEM[" @@ -6856,6 +6867,30 @@ TEST(io, projstringformatter_optim_hgridshift_vgridshift_hgridshift_inv) { "+step +proj=pop +v_1 +v_2"); } + // Test omit_fwd->omit_inv when inversing the pipeline + { + auto fmt = PROJStringFormatter::create(); + fmt->startInversion(); + fmt->ingestPROJString("+proj=hgridshift +grids=foo +omit_fwd"); + fmt->stopInversion(); + + EXPECT_EQ(fmt->toString(), + "+proj=pipeline " + "+step +inv +proj=hgridshift +grids=foo +omit_inv"); + } + + // Test omit_inv->omit_fwd when inversing the pipeline + { + auto fmt = PROJStringFormatter::create(); + fmt->startInversion(); + fmt->ingestPROJString("+proj=hgridshift +grids=foo +omit_inv"); + fmt->stopInversion(); + + EXPECT_EQ(fmt->toString(), + "+proj=pipeline " + "+step +inv +proj=hgridshift +grids=foo +omit_fwd"); + } + // Variant with first hgridshift inverted, and second forward { auto fmt = PROJStringFormatter::create(); |
