From 572433ec4d7ef08d99ef163bbf42bdf69d870985 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 12 Oct 2019 00:15:16 +0200 Subject: createOperations(): allow transforming from a compoundCRS of a bound verticalCRS to a 2D CRS --- src/iso19111/coordinateoperation.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) (limited to 'src') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 3ba5a5a2..915ce88f 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -12619,9 +12619,7 @@ CoordinateOperationFactory::Private::createOperations( auto vertCRSOfBaseOfBoundSrc = dynamic_cast(boundSrc->baseCRS().get()); - if (vertCRSOfBaseOfBoundSrc && hubSrcGeog && - hubSrcGeog->coordinateSystem()->axisList().size() == 3 && - geogDst->coordinateSystem()->axisList().size() == 3) { + if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) { auto opsFirst = createOperations(sourceCRS, hubSrc, context); if (context.skipHorizontalTransformation) { if (!opsFirst.empty()) -- cgit v1.2.3 From c1ef192087eaf7d302065ae9196ca31abde02900 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 3 Oct 2019 12:05:44 +0000 Subject: ob_tran.cpp: add comment to link maths with Snyder's formulas --- src/projections/ob_tran.cpp | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'src') diff --git a/src/projections/ob_tran.cpp b/src/projections/ob_tran.cpp index 4ae8dbe7..badc6dd8 100644 --- a/src/projections/ob_tran.cpp +++ b/src/projections/ob_tran.cpp @@ -31,8 +31,10 @@ static PJ_XY o_forward(PJ_LP lp, PJ *P) { /* spheroid */ coslam = cos(lp.lam); sinphi = sin(lp.phi); cosphi = cos(lp.phi); + /* Formula (5-8b) of Snyder's "Map projections: a working manual" */ lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), Q->sphip * cosphi * coslam + Q->cphip * sinphi) + Q->lamp); + /* Formula (5-7) */ lp.phi = aasin(P->ctx,Q->sphip * sinphi - Q->cphip * cosphi * coslam); return Q->link->fwd(lp, Q->link); @@ -62,7 +64,9 @@ static PJ_LP o_inverse(PJ_XY xy, PJ *P) { /* spheroid */ coslam = cos(lp.lam -= Q->lamp); sinphi = sin(lp.phi); cosphi = cos(lp.phi); + /* Formula (5-9) */ lp.phi = aasin(P->ctx,Q->sphip * sinphi + Q->cphip * cosphi * coslam); + /* Formula (5-10b) */ lp.lam = aatan2(cosphi * sin(lp.lam), Q->sphip * cosphi * coslam - Q->cphip * sinphi); } -- cgit v1.2.3 From f94f1ddeaa83bab4f6f3f9d5720a777e569c781b Mon Sep 17 00:00:00 2001 From: "backporting[bot]" <42222807+backporting[bot]@users.noreply.github.com> Date: Fri, 18 Oct 2019 22:40:53 +0200 Subject: [Backport 6.2] Fix segfaults in case of out-of-memory situations (#1684) (fixes #1678) --- src/datum_set.cpp | 12 ++++++++++++ src/ell_set.cpp | 8 ++++++++ 2 files changed, 20 insertions(+) (limited to 'src') diff --git a/src/datum_set.cpp b/src/datum_set.cpp index c1cb4cb9..873d7be5 100644 --- a/src/datum_set.cpp +++ b/src/datum_set.cpp @@ -85,10 +85,22 @@ int pj_datum_set(projCtx ctx, paralist *pl, PJ *projdef) entry[ sizeof(entry) - 1 ] = '\0'; curr = curr->next = pj_mkparam(entry); + if (nullptr == curr) + { + pj_ctx_set_errno(ctx, ENOMEM); + return 1; + } } if( pj_datums[i].defn && strlen(pj_datums[i].defn) > 0 ) + { curr = curr->next = pj_mkparam(pj_datums[i].defn); + if (nullptr == curr) + { + pj_ctx_set_errno(ctx, ENOMEM); + return 1; + } + } (void)curr; /* make clang static analyzer happy */ } diff --git a/src/ell_set.cpp b/src/ell_set.cpp index d2930ca4..bb46b3a4 100644 --- a/src/ell_set.cpp +++ b/src/ell_set.cpp @@ -1,5 +1,6 @@ /* set ellipsoid parameters a and es */ +#include #include #include #include @@ -156,7 +157,14 @@ static int ellps_ellps (PJ *P) { err = proj_errno_reset (P); paralist* new_params = pj_mkparam (ellps->major); + if (nullptr == new_params) + return proj_errno_set (P, ENOMEM); new_params->next = pj_mkparam (ellps->ell); + if (nullptr == new_params->next) + { + pj_dealloc(new_params); + return proj_errno_set (P, ENOMEM); + } paralist* old_params = P->params; P->params = new_params; -- cgit v1.2.3 From 00bac075a1b626f87262ea5c64a805a104765cba Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 18 Oct 2019 19:54:55 +0000 Subject: createOperations(): fix double vertical unit conversion from CompoundCRS to other CRS when the horizontal part of the projected CRS uses non-metre unit Fix issue reported on https://lists.osgeo.org/pipermail/proj/2019-October/008939.html --- src/iso19111/crs.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index eea7cce8..a13662cd 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -3239,13 +3239,16 @@ void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter, if (!formatter->getCRSExport()) { formatter->addStep("unitconvert"); formatter->addParam("xy_in", "m"); - formatter->addParam("z_in", "m"); + if (!formatter->omitZUnitConversion()) + formatter->addParam("z_in", "m"); if (projUnit.empty()) { formatter->addParam("xy_out", toSI); - formatter->addParam("z_out", toSI); + if (!formatter->omitZUnitConversion()) + formatter->addParam("z_out", toSI); } else { formatter->addParam("xy_out", projUnit); - formatter->addParam("z_out", projUnit); + if (!formatter->omitZUnitConversion()) + formatter->addParam("z_out", projUnit); } } else { if (projUnit.empty()) { -- cgit v1.2.3 From f96f8ef4ffc0002219799cc5b0b2364c94a4ed49 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 25 Oct 2019 09:04:16 +0000 Subject: importFromWkt(): fix axis orientation for non-standard ESRI WKT (fixes #1690) --- src/iso19111/io.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index a9607247..2a6d8a5e 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3185,7 +3185,16 @@ ConversionNNPtr WKTParser::Private::buildProjectionFromESRI( } } - const auto *wkt2_mapping = getMapping(esriMapping->wkt2_name); + const char *projectionMethodWkt2Name = esriMapping->wkt2_name; + if (ci_equal(esriProjectionName, "Krovak")) { + const std::string projCRSName = + stripQuotes(projCRSNode->GP()->children()[0]); + if (projCRSName.find("_East_North") != std::string::npos) { + projectionMethodWkt2Name = EPSG_NAME_METHOD_KROVAK_NORTH_ORIENTED; + } + } + + const auto *wkt2_mapping = getMapping(projectionMethodWkt2Name); if (ci_equal(esriProjectionName, "Stereographic")) { try { if (std::fabs(io::asDouble( -- cgit v1.2.3 From c206d14f681bf9aaf603c445b82a84bbc1978f1e Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 28 Oct 2019 17:53:05 +0100 Subject: Update ABI version numbers for 6.2.1 release --- src/Makefile.am | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index fe7f2572..0a868412 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -42,7 +42,7 @@ geodtest_LDADD = libproj.la lib_LTLIBRARIES = libproj.la -libproj_la_LDFLAGS = -no-undefined -version-info 17:0:2 +libproj_la_LDFLAGS = -no-undefined -version-info 17:1:2 libproj_la_LIBADD = @SQLITE3_LIBS@ libproj_la_SOURCES = \ -- cgit v1.2.3 From 45b36ca27685c13c62c340bbdd3dd449f7e94af9 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 28 Oct 2019 18:50:18 +0100 Subject: Update version numbers for 6.2.1 --- src/proj.h | 2 +- src/proj_api.h | 2 +- src/release.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/proj.h b/src/proj.h index 7a47f043..73c1165d 100644 --- a/src/proj.h +++ b/src/proj.h @@ -154,7 +154,7 @@ extern "C" { /* The version numbers should be updated with every release! **/ #define PROJ_VERSION_MAJOR 6 #define PROJ_VERSION_MINOR 2 -#define PROJ_VERSION_PATCH 0 +#define PROJ_VERSION_PATCH 1 extern char const PROJ_DLL pj_release[]; /* global release id string */ diff --git a/src/proj_api.h b/src/proj_api.h index b933fc6b..23eecfa7 100644 --- a/src/proj_api.h +++ b/src/proj_api.h @@ -38,7 +38,7 @@ #endif #ifndef PJ_VERSION -#define PJ_VERSION 620 +#define PJ_VERSION 621 #endif #ifdef PROJ_RENAME_SYMBOLS diff --git a/src/release.cpp b/src/release.cpp index 13ed45f0..7f62da3f 100644 --- a/src/release.cpp +++ b/src/release.cpp @@ -11,7 +11,7 @@ char const pj_release[] = STR(PROJ_VERSION_MAJOR)"." STR(PROJ_VERSION_MINOR)"." STR(PROJ_VERSION_PATCH)", " - "September 1st, 2019"; + "November 1st, 2019"; const char *pj_get_release() { return pj_release; -- cgit v1.2.3 From 970946b5668b31c1f2498db4b1dbc83871d52a91 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 30 Oct 2019 15:56:41 +0000 Subject: createFromWkt(): be tolerant to missing scale_factor parameter (fixes #1700) This is invalid WKT, but GDAL 2.4 used to accept it and make a reasonable use of it... Currently we default it to 0 which is non sensical. Better use 1 as GDAL 2.4 did, and emit a warning. Other fix: proj_create_from_wkt() was documented to operate by default in non-strict validation mode, but it was actually in strict mode. So do as documented. --- src/iso19111/c_api.cpp | 20 ++++++++++--- src/iso19111/coordinateoperation.cpp | 33 +++++++++++++--------- src/iso19111/io.cpp | 55 +++++++++++++++++++++++++++++++++++- 3 files changed, 90 insertions(+), 18 deletions(-) (limited to 'src') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 1e5c106a..c99eb12c 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -498,7 +498,7 @@ template static PROJ_STRING_LIST to_string_list(T &&set) { * proj_string_list_destroy(). * @param out_grammar_errors Pointer to a PROJ_STRING_LIST object, or NULL. * If provided, *out_grammar_errors will contain a list of errors regarding the - * WKT grammaer. It must be freed with proj_string_list_destroy(). + * WKT grammar. It must be freed with proj_string_list_destroy(). * @return Object that must be unreferenced with proj_destroy(), or NULL in * case of error. */ @@ -522,6 +522,7 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt, if (dbContext) { parser.attachDatabaseContext(NN_NO_CHECK(dbContext)); } + parser.setStrict(false); for (auto iter = options; iter && iter[0]; ++iter) { const char *value; if ((value = getOptionValue(*iter, "STRICT="))) { @@ -536,10 +537,19 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt, auto obj = nn_dynamic_pointer_cast( parser.createFromWKT(wkt)); + std::vector warningsFromParsing; if (out_grammar_errors) { - auto warnings = parser.warningList(); - if (!warnings.empty()) { - *out_grammar_errors = to_string_list(warnings); + auto rawWarnings = parser.warningList(); + std::vector grammarWarnings; + for (const auto &msg : rawWarnings) { + if (msg.find("Default it to") != std::string::npos) { + warningsFromParsing.push_back(msg); + } else { + grammarWarnings.push_back(msg); + } + } + if (!grammarWarnings.empty()) { + *out_grammar_errors = to_string_list(grammarWarnings); } } @@ -548,6 +558,8 @@ PJ *proj_create_from_wkt(PJ_CONTEXT *ctx, const char *wkt, if (derivedCRS) { auto warnings = derivedCRS->derivingConversionRef()->validateParameters(); + warnings.insert(warnings.end(), warningsFromParsing.begin(), + warningsFromParsing.end()); if (!warnings.empty()) { *out_warnings = to_string_list(warnings); } diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 915ce88f..c4e3973e 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -5927,24 +5927,31 @@ void Conversion::_exportToPROJString( if (!param->proj_name) { continue; } - auto value = + const auto value = parameterValueMeasure(param->wkt2_name, param->epsg_code); + double valueConverted = 0; + if (value == nullMeasure) { + // Deal with missing values. In an ideal world, this would + // not happen + if (param->epsg_code == + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) { + valueConverted = 1.0; + } + } else if (param->unit_type == + common::UnitOfMeasure::Type::ANGULAR) { + valueConverted = + value.convertToUnit(common::UnitOfMeasure::DEGREE); + } else { + valueConverted = value.getSIValue(); + } + if (mapping->epsg_code == EPSG_CODE_METHOD_LAMBERT_CONIC_CONFORMAL_1SP && strcmp(param->proj_name, "lat_1") == 0) { - formatter->addParam( - param->proj_name, - value.convertToUnit(common::UnitOfMeasure::DEGREE)); - formatter->addParam( - "lat_0", - value.convertToUnit(common::UnitOfMeasure::DEGREE)); - } else if (param->unit_type == - common::UnitOfMeasure::Type::ANGULAR) { - formatter->addParam( - param->proj_name, - value.convertToUnit(common::UnitOfMeasure::DEGREE)); + formatter->addParam(param->proj_name, valueConverted); + formatter->addParam("lat_0", valueConverted); } else { - formatter->addParam(param->proj_name, value.getSIValue()); + formatter->addParam(param->proj_name, valueConverted); } } diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 2a6d8a5e..c0628270 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -3476,6 +3476,14 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( } propertiesMethod.set(IdentifiedObject::NAME_KEY, projectionName); + std::vector foundParameters; + if (mapping) { + size_t countParams = 0; + while (mapping->params[countParams] != nullptr) { + ++countParams; + } + foundParameters.resize(countParams); + } for (const auto &childNode : projCRSNode->GP()->children()) { if (ci_equal(childNode->GP()->value(), WKTConstants::PARAMETER)) { const auto &childNodeChildren = childNode->GP()->children(); @@ -3496,11 +3504,18 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( continue; } } - const auto *paramMapping = + auto *paramMapping = mapping ? getMappingFromWKT1(mapping, parameterName) : nullptr; if (mapping && mapping->epsg_code == EPSG_CODE_METHOD_MERCATOR_VARIANT_B && ci_equal(parameterName, "latitude_of_origin")) { + for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) { + if (mapping->params[idx]->epsg_code == + EPSG_CODE_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN) { + foundParameters[idx] = true; + break; + } + } parameterName = EPSG_NAME_PARAMETER_LATITUDE_OF_NATURAL_ORIGIN; propertiesParameter.set( Identifier::CODE_KEY, @@ -3508,6 +3523,12 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( propertiesParameter.set(Identifier::CODESPACE_KEY, Identifier::EPSG); } else if (paramMapping) { + for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) { + if (mapping->params[idx] == paramMapping) { + foundParameters[idx] = true; + break; + } + } parameterName = paramMapping->wkt2_name; if (paramMapping->epsg_code != 0) { propertiesParameter.set(Identifier::CODE_KEY, @@ -3531,6 +3552,38 @@ ConversionNNPtr WKTParser::Private::buildProjectionStandard( } } + // Add back important parameters that should normally be present, but + // are sometimes missing. Currently we only deal with Scale factor at + // natural origin. This is to avoid a default value of 0 to slip in later. + // But such WKT should be considered invalid. + if (mapping) { + for (size_t idx = 0; mapping->params[idx] != nullptr; ++idx) { + if (!foundParameters[idx] && + mapping->params[idx]->epsg_code == + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN) { + + emitRecoverableWarning( + "The WKT string lacks a value " + "for " EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN + ". Default it to 1."); + + PropertyMap propertiesParameter; + propertiesParameter.set( + Identifier::CODE_KEY, + EPSG_CODE_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN); + propertiesParameter.set(Identifier::CODESPACE_KEY, + Identifier::EPSG); + propertiesParameter.set( + IdentifiedObject::NAME_KEY, + EPSG_NAME_PARAMETER_SCALE_FACTOR_AT_NATURAL_ORIGIN); + parameters.push_back( + OperationParameter::create(propertiesParameter)); + values.push_back(ParameterValue::create( + Measure(1.0, UnitOfMeasure::SCALE_UNITY))); + } + } + } + return Conversion::create( PropertyMap().set(IdentifiedObject::NAME_KEY, "unnamed"), propertiesMethod, parameters, values) -- cgit v1.2.3 From 97a468bbbea4a18126f2b5e1516183f30fa0ac6b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Nov 2019 14:33:55 +0000 Subject: createOperations(): in some circumstances we wrongly promoted a Helmert geog2D transformation to a geog3D Fixes for example EPSG:4979 to EPSG:2189, as raised in https://github.com/OSGeo/gdal/issues/1972#issuecomment-548814354 --- src/iso19111/coordinateoperation.cpp | 51 +++++++++++++++++++++++++++++------- 1 file changed, 41 insertions(+), 10 deletions(-) (limited to 'src') diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index c4e3973e..0eaba883 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -6513,6 +6513,24 @@ static void getTransformationType(const crs::CRSNNPtr &sourceCRSIn, isGeog2D = nSrcAxisCount == 2 && nTargetAxisCount == 2; isGeog3D = !isGeog2D && nSrcAxisCount >= 2 && nTargetAxisCount >= 2; } + +// --------------------------------------------------------------------------- + +static int +useOperationMethodEPSGCodeIfPresent(const util::PropertyMap &properties, + int nDefaultOperationMethodEPSGCode) { + const auto *operationMethodEPSGCode = + properties.get("OPERATION_METHOD_EPSG_CODE"); + if (operationMethodEPSGCode) { + const auto boxedValue = dynamic_cast( + (*operationMethodEPSGCode).get()); + if (boxedValue && + boxedValue->type() == util::BoxedValue::Type::INTEGER) { + return boxedValue->integerValue(); + } + } + return nDefaultOperationMethodEPSGCode; +} //! @endcond // --------------------------------------------------------------------------- @@ -6541,12 +6559,13 @@ TransformationNNPtr Transformation::createGeocentricTranslations( isGeog3D); return create( properties, sourceCRSIn, targetCRSIn, nullptr, - createMethodMapNameEPSGCode( + createMethodMapNameEPSGCode(useOperationMethodEPSGCodeIfPresent( + properties, isGeocentric ? EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOCENTRIC : isGeog2D ? EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_2D - : EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_3D), + : EPSG_CODE_METHOD_GEOCENTRIC_TRANSLATION_GEOGRAPHIC_3D)), VectorOfParameters{ createOpParamNameEPSGCode(EPSG_CODE_PARAMETER_X_AXIS_TRANSLATION), createOpParamNameEPSGCode(EPSG_CODE_PARAMETER_Y_AXIS_TRANSLATION), @@ -6599,11 +6618,12 @@ TransformationNNPtr Transformation::createPositionVector( isGeog3D); return createSevenParamsTransform( properties, - createMethodMapNameEPSGCode( + createMethodMapNameEPSGCode(useOperationMethodEPSGCodeIfPresent( + properties, isGeocentric ? EPSG_CODE_METHOD_POSITION_VECTOR_GEOCENTRIC : isGeog2D ? EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_2D - : EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_3D), + : EPSG_CODE_METHOD_POSITION_VECTOR_GEOGRAPHIC_3D)), sourceCRSIn, targetCRSIn, translationXMetre, translationYMetre, translationZMetre, rotationXArcSecond, rotationYArcSecond, rotationZArcSecond, scaleDifferencePPM, accuracies); @@ -6648,11 +6668,12 @@ TransformationNNPtr Transformation::createCoordinateFrameRotation( isGeog3D); return createSevenParamsTransform( properties, - createMethodMapNameEPSGCode( + createMethodMapNameEPSGCode(useOperationMethodEPSGCodeIfPresent( + properties, isGeocentric ? EPSG_CODE_METHOD_COORDINATE_FRAME_GEOCENTRIC : isGeog2D ? EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_2D - : EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_3D), + : EPSG_CODE_METHOD_COORDINATE_FRAME_GEOGRAPHIC_3D)), sourceCRSIn, targetCRSIn, translationXMetre, translationYMetre, translationZMetre, rotationXArcSecond, rotationYArcSecond, rotationZArcSecond, scaleDifferencePPM, accuracies); @@ -6791,12 +6812,13 @@ TransformationNNPtr Transformation::createTimeDependentPositionVector( isGeog3D); return createFifteenParamsTransform( properties, - createMethodMapNameEPSGCode( + createMethodMapNameEPSGCode(useOperationMethodEPSGCodeIfPresent( + properties, isGeocentric ? EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOCENTRIC : isGeog2D ? EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_2D - : EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_3D), + : EPSG_CODE_METHOD_TIME_DEPENDENT_POSITION_VECTOR_GEOGRAPHIC_3D)), sourceCRSIn, targetCRSIn, translationXMetre, translationYMetre, translationZMetre, rotationXArcSecond, rotationYArcSecond, rotationZArcSecond, scaleDifferencePPM, rateTranslationX, @@ -6868,12 +6890,13 @@ TransformationNNPtr Transformation::createTimeDependentCoordinateFrameRotation( isGeog3D); return createFifteenParamsTransform( properties, - createMethodMapNameEPSGCode( + createMethodMapNameEPSGCode(useOperationMethodEPSGCodeIfPresent( + properties, isGeocentric ? EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOCENTRIC : isGeog2D ? EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_2D - : EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_3D), + : EPSG_CODE_METHOD_TIME_DEPENDENT_COORDINATE_FRAME_GEOGRAPHIC_3D)), sourceCRSIn, targetCRSIn, translationXMetre, translationYMetre, translationZMetre, rotationXArcSecond, rotationYArcSecond, rotationZArcSecond, scaleDifferencePPM, rateTranslationX, @@ -7455,6 +7478,14 @@ createPropertiesForInverse(const CoordinateOperation *op, bool derivedFrom, addModifiedIdentifier(map, op, true, derivedFrom); + const auto so = dynamic_cast(op); + if (so) { + const int soMethodEPSGCode = so->method()->getEPSGCode(); + if (soMethodEPSGCode > 0) { + map.set("OPERATION_METHOD_EPSG_CODE", soMethodEPSGCode); + } + } + return map; } -- cgit v1.2.3 From 8fdd512ea8d9e0e788d93c5f356bc3d1728ba2aa Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 9 Nov 2019 12:22:48 +0000 Subject: Doc: document oddity related to identification of CRS from ESRI WKT Or more generally formulations that don't have an explicit axis order. Refs https://github.com/pyproj4/pyproj/issues/475 projinfo 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]' returns EPSG:4326 with 100% confidence. But its axis order is not the same as EPSG:4326. I've pondered about this, like decreasing the confidence of the match, but this would have downstream effects on GDAL (shapefiles with the above content in a .prj would no longer be identified as EPSG:4326). So for now, document that oddity. --- src/iso19111/c_api.cpp | 24 ++++++++++++++++------ src/iso19111/crs.cpp | 55 ++++++++++++++++++++++++++++++++++++-------------- 2 files changed, 58 insertions(+), 21 deletions(-) (limited to 'src') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index c99eb12c..fb3b7f40 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -2101,14 +2101,26 @@ PJ *proj_get_target_crs(PJ_CONTEXT *ctx, const PJ *obj) { * The method returns a list of matching reference CRS, and the percentage * (0-100) of confidence in the match. The list is sorted by decreasing * confidence. - * - * 100% means that the name of the reference entry + *
    + *
  • 100% means that the name of the reference entry * perfectly matches the CRS name, and both are equivalent. In which case a * single result is returned. - * 90% means that CRS are equivalent, but the names are not exactly the same. - * 70% means that CRS are equivalent), but the names do not match at all. - * 25% means that the CRS are not equivalent, but there is some similarity in - * the names. + * Note: in the case of a GeographicCRS whose axis + * order is implicit in the input definition (for example ESRI WKT), then axis + * order is ignored for the purpose of identification. That is the CRS built + * from + * GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]], + * PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] + * will be identified to EPSG:4326, but will not pass a + * isEquivalentTo(EPSG_4326, util::IComparable::Criterion::EQUIVALENT) test, + * but rather isEquivalentTo(EPSG_4326, + * util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS) + *
    • + *
  • 90% means that CRS are equivalent, but the names are not exactly the same.
    • + *
  • 70% means that CRS are equivalent), but the names do not match at all.
    • + *
  • 25% means that the CRS are not equivalent, but there is some similarity in + * the names.
    • + *
* Other confidence values may be returned by some specialized implementations. * * This is implemented for GeodeticCRS, ProjectedCRS, VerticalCRS and diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index a13662cd..46cde43a 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -697,14 +697,26 @@ CRSNNPtr CRS::normalizeForVisualization() const { * The method returns a list of matching reference CRS, and the percentage * (0-100) of confidence in the match. The list is sorted by decreasing * confidence. - * - * 100% means that the name of the reference entry + *
    + *
  • 100% means that the name of the reference entry * perfectly matches the CRS name, and both are equivalent. In which case a * single result is returned. - * 90% means that CRS are equivalent, but the names are not exactly the same. - * 70% means that CRS are equivalent), but the names do not match at all. - * 25% means that the CRS are not equivalent, but there is some similarity in - * the names. + * Note: in the case of a GeographicCRS whose axis + * order is implicit in the input definition (for example ESRI WKT), then axis + * order is ignored for the purpose of identification. That is the CRS built + * from + * GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]], + * PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] + * will be identified to EPSG:4326, but will not pass a + * isEquivalentTo(EPSG_4326, util::IComparable::Criterion::EQUIVALENT) test, + * but rather isEquivalentTo(EPSG_4326, + * util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS) + *
    • + *
  • 90% means that CRS are equivalent, but the names are not exactly the same.
    • + *
  • 70% means that CRS are equivalent), but the names do not match at all.
    • + *
  • 25% means that the CRS are not equivalent, but there is some similarity in + * the names.
    • + *
* Other confidence values may be returned by some specialized implementations. * * This is implemented for GeodeticCRS, ProjectedCRS, VerticalCRS and @@ -1554,17 +1566,30 @@ static bool hasCodeCompatibleOfAuthorityFactory( * authorityFactory is not null. * * The method returns a list of matching reference CRS, and the percentage - * (0-100) of confidence in the match. - * 100% means that the name of the reference entry + * (0-100) of confidence in the match: + *
    + *
  • 100% means that the name of the reference entry * perfectly matches the CRS name, and both are equivalent. In which case a * single result is returned. - * 90% means that CRS are equivalent, but the names are not exactly the same. - * 70% means that CRS are equivalent (equivalent datum and coordinate system), - * but the names do not match at all. - * 60% means that ellipsoid, prime meridian and coordinate systems are - * equivalent, but the CRS and datum names do not match. - * 25% means that the CRS are not equivalent, but there is some similarity in - * the names. + * Note: in the case of a GeographicCRS whose axis + * order is implicit in the input definition (for example ESRI WKT), then axis + * order is ignored for the purpose of identification. That is the CRS built + * from + * GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]], + * PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] + * will be identified to EPSG:4326, but will not pass a + * isEquivalentTo(EPSG_4326, util::IComparable::Criterion::EQUIVALENT) test, + * but rather isEquivalentTo(EPSG_4326, + * util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS) + *
    • + *
  • 90% means that CRS are equivalent, but the names are not exactly the same. + *
  • 70% means that CRS are equivalent (equivalent datum and coordinate system), + * but the names do not match at all.
    • + *
  • 60% means that ellipsoid, prime meridian and coordinate systems are + * equivalent, but the CRS and datum names do not match.
    • + *
  • 25% means that the CRS are not equivalent, but there is some similarity in + * the names.
    • + *
* * @param authorityFactory Authority factory (or null, but degraded * functionality) -- cgit v1.2.3 From f970636daa814be26586bba1cfaa9fbb132da2df Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 12 Nov 2019 12:33:47 +0000 Subject: Fix proj_assign_context()/pj_set_ctx() with pipelines and alternative coord operations Fixes https://github.com/OSGeo/gdal/issues/1989 pj_set_ctx() only changes the context to the main object. It should also recurse down to the steps of the pipeline and the alternative coordinate operations hold in alternativeCoordinateOperations In the GDAL use case with multithreaded reprojection, and objects being transferred between thread, this would cause a failed coordinate transformation to affect an unrelated transformation of another thread... --- src/ctx.cpp | 8 ++++++++ src/pipeline.cpp | 7 ++++++- src/proj_internal.h | 2 ++ 3 files changed, 16 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/ctx.cpp b/src/ctx.cpp index 622a814d..bcb6e1cc 100644 --- a/src/ctx.cpp +++ b/src/ctx.cpp @@ -60,6 +60,14 @@ void pj_set_ctx( projPJ pj, projCtx ctx ) if (pj==nullptr) return; pj->ctx = ctx; + if( pj->is_pipeline ) + { + pj_pipeline_assign_context_to_steps(pj, ctx); + } + for( const auto &alt: pj->alternativeCoordinateOperations ) + { + pj_set_ctx(alt.pj, ctx); + } } /************************************************************************/ diff --git a/src/pipeline.cpp b/src/pipeline.cpp index afa3b19a..76ffc98e 100644 --- a/src/pipeline.cpp +++ b/src/pipeline.cpp @@ -105,6 +105,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20 #include "geodesic.h" #include "proj.h" #include "proj_internal.h" +#include "proj_experimental.h" PROJ_HEAD(pipeline, "Transformation pipeline manager"); PROJ_HEAD(pop, "Retrieve coordinate value from pipeline stack"); @@ -136,7 +137,11 @@ static PJ_LPZ pipeline_reverse_3d (PJ_XYZ xyz, PJ *P); static PJ_XY pipeline_forward (PJ_LP lp, PJ *P); static PJ_LP pipeline_reverse (PJ_XY xy, PJ *P); - +void pj_pipeline_assign_context_to_steps( PJ* P, PJ_CONTEXT* ctx ) +{ + for (int i = 1; i <= static_cast(P->opaque)->steps; i++) + proj_assign_context(static_cast(P->opaque)->pipeline[i], ctx); +} static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P) { diff --git a/src/proj_internal.h b/src/proj_internal.h index 4a126e98..77c4ba17 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -926,6 +926,8 @@ std::string pj_double_quote_string_param_if_needed(const std::string& str); PJ *pj_create_internal (PJ_CONTEXT *ctx, const char *definition); PJ *pj_create_argv_internal (PJ_CONTEXT *ctx, int argc, char **argv); +void pj_pipeline_assign_context_to_steps( PJ* P, PJ_CONTEXT* ctx ); + /* classic public API */ #include "proj_api.h" -- cgit v1.2.3 From bce4b158ab5f7d146de8e8fc98df4612dc8c2c9e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 25 Nov 2019 04:54:17 +0100 Subject: normalizeForVisualization() and other methods applying on a ProjectedCRS: do not mess the derivingConversion object of the original object (fixes #1736) normalizeForVisualization(), promoteTo3D(), demoteTo2D(), alterGeodeticCRS(), alterCSLinearUnit() and alterParametersLinearUnit() all used the object returned by derivingConversionRef() to create a new ProjectedCRS. While doing so, this caused the derivingConversion of the original object to have its targetCRS set to the object returned by normalizeForVisualization() and similar. If that object died, then the weak pointer would be reset, and the original ProjectedCRS() has now its derivingConversionRef()->targetCRS() nullptr So bottom line is use derivingConversion() for anything that is not pure reading !!! This is confirmed to be the fix for the QGIS scenario in https://github.com/qgis/QGIS/issues/30569#issuecomment-540060919 In QGIS use case, the issue arised when using a projected CRS with a non-GIS friendly axis (that is where normalizeForVisualization() created a new projectedCRS) --- src/iso19111/crs.cpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) (limited to 'src') diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp index 46cde43a..e044c0c7 100644 --- a/src/iso19111/crs.cpp +++ b/src/iso19111/crs.cpp @@ -237,7 +237,7 @@ CRSNNPtr CRS::alterGeodeticCRS(const GeodeticCRSNNPtr &newGeodCRS) const { auto projCRS = dynamic_cast(this); if (projCRS) { return ProjectedCRS::create(createPropertyMap(this), newGeodCRS, - projCRS->derivingConversionRef(), + projCRS->derivingConversion(), projCRS->coordinateSystem()); } @@ -264,7 +264,7 @@ CRSNNPtr CRS::alterCSLinearUnit(const common::UnitOfMeasure &unit) const { if (projCRS) { return ProjectedCRS::create( createPropertyMap(this), projCRS->baseCRS(), - projCRS->derivingConversionRef(), + projCRS->derivingConversion(), projCRS->coordinateSystem()->alterUnit(unit)); } } @@ -675,9 +675,8 @@ CRSNNPtr CRS::normalizeForVisualization() const { axisList[0]) : cs::CartesianCS::create(util::PropertyMap(), axisList[1], axisList[0], axisList[2]); - return util::nn_static_pointer_cast( - ProjectedCRS::create(props, projCRS->baseCRS(), - projCRS->derivingConversionRef(), cs)); + return util::nn_static_pointer_cast(ProjectedCRS::create( + props, projCRS->baseCRS(), projCRS->derivingConversion(), cs)); } } @@ -841,7 +840,7 @@ CRSNNPtr CRS::promoteTo3D(const std::string &newName, !newName.empty() ? newName : nameStr()), NN_NO_CHECK( util::nn_dynamic_pointer_cast(base3DCRS)), - projCRS->derivingConversionRef(), cs)); + projCRS->derivingConversion(), cs)); } } @@ -3243,10 +3242,10 @@ bool ProjectedCRS::_isEquivalentTo( ProjectedCRSNNPtr ProjectedCRS::alterParametersLinearUnit(const common::UnitOfMeasure &unit, bool convertToNewUnit) const { - return create(createPropertyMap(this), baseCRS(), - derivingConversionRef()->alterParametersLinearUnit( - unit, convertToNewUnit), - coordinateSystem()); + return create( + createPropertyMap(this), baseCRS(), + derivingConversion()->alterParametersLinearUnit(unit, convertToNewUnit), + coordinateSystem()); } //! @endcond -- cgit v1.2.3