diff options
| -rw-r--r-- | src/iso19111/coordinateoperation.cpp | 87 | ||||
| -rw-r--r-- | test/unit/test_operation.cpp | 95 |
2 files changed, 149 insertions, 33 deletions
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index ef262cc8..6593d3e9 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -608,6 +608,7 @@ struct CoordinateOperation::Private { interpolationCRS_(other.interpolationCRS_), sourceCoordinateEpoch_(other.sourceCoordinateEpoch_), targetCoordinateEpoch_(other.targetCoordinateEpoch_), + hasBallparkTransformation_(other.hasBallparkTransformation_), strongRef_(other.strongRef_ ? internal::make_unique<CRSStrongRef>( *(other.strongRef_)) : nullptr) {} @@ -11038,7 +11039,7 @@ struct CoordinateOperationFactory::Private { const crs::CRSNNPtr &targetCRS, Context &context); private: - static constexpr bool allowEmptyIntersection = true; + static constexpr bool disallowEmptyIntersection = true; static void buildCRSIds(const crs::CRSNNPtr &crs, Private::Context &context, @@ -12553,7 +12554,8 @@ createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc, static CoordinateOperationNNPtr createHorizVerticalPROJBased( const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS, const operation::CoordinateOperationNNPtr &horizTransform, - const operation::CoordinateOperationNNPtr &verticalTransform) { + const operation::CoordinateOperationNNPtr &verticalTransform, + bool checkExtent) { auto geogDst = util::nn_dynamic_pointer_cast<crs::GeographicCRS>(targetCRS); assert(geogDst); @@ -12579,10 +12581,16 @@ static CoordinateOperationNNPtr createHorizVerticalPROJBased( verticalTransform->coordinateOperationAccuracies(), verticalTransform->hasBallparkTransformation()); } else { - bool dummy = false; + bool emptyIntersection = false; auto ops = std::vector<CoordinateOperationNNPtr>{horizTransform, verticalTransform}; - auto extent = getExtent(ops, true, dummy); + auto extent = getExtent(ops, true, emptyIntersection); + if (checkExtent && emptyIntersection) { + std::string msg( + "empty intersection of area of validity of concatenated " + "operations"); + throw InvalidOperationEmptyIntersection(msg); + } auto properties = util::PropertyMap(); properties.set(common::IdentifiedObject::NAME_KEY, computeConcatenatedName(ops)); @@ -12844,7 +12852,7 @@ CoordinateOperationFactory::Private::createOperationsGeogToGeog( } auto op = ConcatenatedOperation::createComputeMetadata( - steps, !allowEmptyIntersection); + steps, disallowEmptyIntersection); op->setHasBallparkTransformation(!sameDatum); res.emplace_back(op); return res; @@ -13113,8 +13121,11 @@ void CoordinateOperationFactory::Private::createOperationsWithDatumPivot( } logTrace("transformation " + debugStr); #endif - res.emplace_back(ConcatenatedOperation::createComputeMetadata( - subOps, !allowEmptyIntersection)); + try { + res.emplace_back(ConcatenatedOperation::createComputeMetadata( + subOps, disallowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } } }; @@ -13700,7 +13711,7 @@ CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid( context, targetOp, vertDst, tmp); assert(!tmp.empty()); auto ret = ConcatenatedOperation::createComputeMetadata( - {op, tmp.front()}, !allowEmptyIntersection); + {op, tmp.front()}, disallowEmptyIntersection); return ret; }; @@ -13858,7 +13869,7 @@ std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private:: res.emplace_back( ConcatenatedOperation::createComputeMetadata( {opFirst, opsSecond.front()}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } } } @@ -14043,7 +14054,7 @@ void CoordinateOperationFactory::Private::createOperationsGeodToGeod( auto opSecond = Conversion::createGeographicGeocentric(interm_crs, targetCRS); res.emplace_back(ConcatenatedOperation::createComputeMetadata( - {opFirst, opSecond}, !allowEmptyIntersection)); + {opFirst, opSecond}, disallowEmptyIntersection)); } else { // Apply previous case in reverse way std::vector<CoordinateOperationNNPtr> resTmp; @@ -14109,7 +14120,7 @@ void CoordinateOperationFactory::Private::createOperationsDerivedTo( for (const auto &opSecond : opsSecond) { try { res.emplace_back(ConcatenatedOperation::createComputeMetadata( - {opFirst, opSecond}, !allowEmptyIntersection)); + {opFirst, opSecond}, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14168,7 +14179,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( ConcatenatedOperation::createComputeMetadata( {NN_NO_CHECK(opIntermediate), boundSrc->transformation()}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } else { @@ -14191,7 +14202,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( subops.emplace_back(boundSrc->transformation()); res.emplace_back( ConcatenatedOperation::createComputeMetadata( - subops, !allowEmptyIntersection)); + subops, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14231,7 +14242,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( subops.emplace_back(opLast); res.emplace_back( ConcatenatedOperation::createComputeMetadata( - subops, !allowEmptyIntersection)); + subops, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14282,7 +14293,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( try { res.emplace_back( ConcatenatedOperation::createComputeMetadata( - {opFirst, transf}, !allowEmptyIntersection)); + {opFirst, transf}, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14308,7 +14319,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( res.emplace_back( ConcatenatedOperation::createComputeMetadata( {opFirst, boundSrc->transformation()}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14329,15 +14340,19 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( for (const auto &opFirst : opsFirst) { for (const auto &opLast : opsLast) { // Exclude artificial transformations from the hub - // to the target CRS - if (!opLast->hasBallparkTransformation()) { + // to the target CRS, if it is the only one. + if (opsLast.size() > 1 || + !opLast->hasBallparkTransformation()) { try { res.emplace_back( ConcatenatedOperation::createComputeMetadata( {opFirst, opLast}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } + } else { + // std::cerr << "excluded " << opLast->nameStr() << + // std::endl; } } } @@ -14395,7 +14410,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( try { res.emplace_back( ConcatenatedOperation::createComputeMetadata( - {op, conv}, !allowEmptyIntersection)); + {op, conv}, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14417,10 +14432,13 @@ void CoordinateOperationFactory::Private::createOperationsBoundToGeog( ConcatenatedOperation:: createComputeMetadata( {opFirst, opLast}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } catch ( const InvalidOperationEmptyIntersection &) { } + } else { + // std::cerr << "excluded " << opLast->nameStr() << + // std::endl; } } } @@ -14625,7 +14643,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToBound( } res.emplace_back( ConcatenatedOperation::createComputeMetadata( - ops, !allowEmptyIntersection)); + ops, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14672,7 +14690,7 @@ void CoordinateOperationFactory::Private::createOperationsBoundToBound( try { res.emplace_back( ConcatenatedOperation::createComputeMetadata( - {opFirst, opLast}, !allowEmptyIntersection)); + {opFirst, opLast}, disallowEmptyIntersection)); } catch (const InvalidOperationEmptyIntersection &) { } } @@ -14844,7 +14862,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( res.emplace_back( ConcatenatedOperation::createComputeMetadata( {opsFirst.front(), opLast}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } catch (const std::exception &) { } } @@ -15082,7 +15100,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToGeog( try { auto op = createHorizVerticalPROJBased( sourceCRS, targetCRS, horizTransform, - verticalTransform); + verticalTransform, disallowEmptyIntersection); res.emplace_back(op); } catch (const std::exception &) { } @@ -15117,9 +15135,13 @@ void CoordinateOperationFactory::Private::createOperationsToGeod( for (const auto &op : sourceToGeog3DOps) { auto newOp = op->shallowClone(); setCRSs(newOp.get(), sourceCRS, intermGeog3DCRS); - res.emplace_back(ConcatenatedOperation::createComputeMetadata( - {newOp, geog3DToTargetOps.front()}, - !allowEmptyIntersection)); + try { + res.emplace_back( + ConcatenatedOperation::createComputeMetadata( + {newOp, geog3DToTargetOps.front()}, + disallowEmptyIntersection)); + } catch (const InvalidOperationEmptyIntersection &) { + } } } } @@ -15186,11 +15208,14 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( createOperations(intermGeogSrc, intermGeogDst, context); for (const auto &op1 : opsSrcToGeog) { if (op1->hasBallparkTransformation()) { + // std::cerr << "excluded " << op1->nameStr() << std::endl; continue; } for (const auto &op2 : opsGeogSrcToGeogDst) { for (const auto &op3 : opsGeogToTarget) { if (op3->hasBallparkTransformation()) { + // std::cerr << "excluded " << op3->nameStr() << + // std::endl; continue; } try { @@ -15204,7 +15229,7 @@ void CoordinateOperationFactory::Private::createOperationsCompoundToCompound( CoordinateOperationNNPtr>{op1, op2, op3}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } catch (const std::exception &) { } } @@ -15340,14 +15365,14 @@ void CoordinateOperationFactory::Private::createOperationsBoundToCompound( createComputeMetadata( {opSrc, intermOps.front(), opDst}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } } else { res.emplace_back( ConcatenatedOperation:: createComputeMetadata( {opSrc, opDst}, - !allowEmptyIntersection)); + disallowEmptyIntersection)); } } } diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 23f3489b..c7c0a262 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -6801,7 +6801,7 @@ TEST(operation, nadgrids_with_pm) { dst = authFactory->createCoordinateReferenceSystem("4258"); list = CoordinateOperationFactory::create()->createOperations( NN_NO_CHECK(src), dst, ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_GE(list.size(), 1U); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " "+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 " @@ -6819,7 +6819,7 @@ TEST(operation, nadgrids_with_pm) { ASSERT_TRUE(crsFromWkt); list = CoordinateOperationFactory::create()->createOperations( NN_NO_CHECK(crsFromWkt), dst, ctxt); - ASSERT_EQ(list.size(), 1U); + ASSERT_GE(list.size(), 1U); EXPECT_EQ(list[0]->exportToPROJString(PROJStringFormatter::create().get()), "+proj=pipeline " "+step +inv +proj=tmerc +lat_0=39.6666666666667 +lon_0=1 " @@ -7621,6 +7621,97 @@ TEST( // --------------------------------------------------------------------------- +TEST(operation, + compoundCRS_with_boundGeogCRS_and_geoid_to_geodCRS_NAD2011_ctxt) { + + auto dbContext = DatabaseContext::create(); + + const char *wktSrc = + "COMPD_CS[\"NAD83 / California zone 5 (ftUS) + " + "NAVD88 height - Geoid12B (ftUS)\"," + " PROJCS[\"NAD83 / California zone 5 (ftUS)\"," + " GEOGCS[\"NAD83\"," + " DATUM[\"North_American_Datum_1983\"," + " SPHEROID[\"GRS 1980\",6378137,298.257222101," + " AUTHORITY[\"EPSG\",\"7019\"]]," + " TOWGS84[0,0,0,0,0,0,0]," + " AUTHORITY[\"EPSG\",\"6269\"]]," + " PRIMEM[\"Greenwich\",0," + " AUTHORITY[\"EPSG\",\"8901\"]]," + " UNIT[\"degree\",0.0174532925199433," + " AUTHORITY[\"EPSG\",\"9122\"]]," + " AUTHORITY[\"EPSG\",\"4269\"]]," + " PROJECTION[\"Lambert_Conformal_Conic_2SP\"]," + " PARAMETER[\"standard_parallel_1\",35.46666666666667]," + " PARAMETER[\"standard_parallel_2\",34.03333333333333]," + " PARAMETER[\"latitude_of_origin\",33.5]," + " PARAMETER[\"central_meridian\",-118]," + " PARAMETER[\"false_easting\",6561666.667]," + " PARAMETER[\"false_northing\",1640416.667]," + " UNIT[\"US survey foot\",0.3048006096012192," + " AUTHORITY[\"EPSG\",\"9003\"]]," + " AXIS[\"X\",EAST]," + " AXIS[\"Y\",NORTH]," + " AUTHORITY[\"EPSG\",\"2229\"]]," + "VERT_CS[\"NAVD88 height - Geoid12B (ftUS)\"," + " VERT_DATUM[\"North American Vertical Datum 1988\",2005," + " AUTHORITY[\"EPSG\",\"5103\"]]," + " UNIT[\"US survey foot\",0.3048006096012192," + " AUTHORITY[\"EPSG\",\"9003\"]]," + " AXIS[\"Gravity-related height\",UP]," + " AUTHORITY[\"EPSG\",\"6360\"]]]"; + auto objSrc = + WKTParser().attachDatabaseContext(dbContext).createFromWKT(wktSrc); + auto srcCRS = nn_dynamic_pointer_cast<CompoundCRS>(objSrc); + ASSERT_TRUE(srcCRS != nullptr); + + auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG"); + // NAD83(2011) geocentric + auto dstCRS = authFactoryEPSG->createCoordinateReferenceSystem("6317"); + + auto authFactory = AuthorityFactory::create(dbContext, std::string()); + 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), dstCRS, ctxt); + bool found = false; + for (const auto &op : list) { + if (op->nameStr() == + "Inverse of unnamed + " + "Transformation from NAD83 to WGS84 + " + "Ballpark geographic offset from WGS 84 to NAD83(2011) + " + "Transformation from NAVD88 height (ftUS) to NAVD88 height + " + "Inverse of NAD83(2011) to NAVD88 height (1) + " + "Conversion from NAD83(2011) (geog3D) to NAD83(2011) " + "(geocentric)") { + found = true; + EXPECT_EQ( + op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=unitconvert +xy_in=us-ft +xy_out=m " + "+step +inv +proj=lcc +lat_0=33.5 +lon_0=-118 " + "+lat_1=35.4666666666667 +lat_2=34.0333333333333 " + "+x_0=2000000.0001016 +y_0=500000.0001016 +ellps=GRS80 " + "+step +proj=unitconvert +z_in=us-ft +z_out=m " + "+step +proj=vgridshift +grids=us_noaa_g2012bu0.tif " + "+multiplier=1 " + "+step +proj=cart +ellps=GRS80"); + } + } + EXPECT_TRUE(found); + if (!found) { + for (const auto &op : list) { + std::cerr << op->nameStr() << std::endl; + } + } +} + +// --------------------------------------------------------------------------- + TEST(operation, geocent_to_compoundCRS) { auto objSrc = PROJStringParser().createFromPROJString( "+proj=geocent +datum=WGS84 +units=m +type=crs"); |
