Skip to content

Commit

Permalink
Merge pull request #2449 from rouault/geoidmodel_ballpark
Browse files Browse the repository at this point in the history 
createOperation(): add a ballpark vertical transformation when dealing with GEOIDMODEL[]
  • Loading branch information
@rouault
rouault authored Nov 24, 2020
2 parents 28aa4f1 + d6f6862 commit 5172091
Show file tree
Hide file tree
Showing 4 changed files with 194 additions and 35 deletions.
99 changes: 85 additions & 14 deletions src/iso19111/coordinateoperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11233,6 +11233,12 @@ struct CoordinateOperationFactory::Private {
const crs::GeographicCRS *geogDst,
std::vector<CoordinateOperationNNPtr> &res);

static void createOperationsVertToGeogBallpark(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &context, const crs::VerticalCRS *vertSrc,
const crs::GeographicCRS *geogDst,
std::vector<CoordinateOperationNNPtr> &res);

static void createOperationsBoundToBound(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &context, const crs::BoundCRS *boundSrc,
Expand Down Expand Up @@ -11334,17 +11340,20 @@ struct PrecomputedOpCharacteristics {
bool gridsKnown_ = false;
size_t stepCount_ = 0;
bool isApprox_ = false;
bool hasBallparkVertical_ = false;
bool isNullTransformation_ = false;

PrecomputedOpCharacteristics() = default;
PrecomputedOpCharacteristics(double area, double accuracy,
bool isPROJExportable, bool hasGrids,
bool gridsAvailable, bool gridsKnown,
size_t stepCount, bool isApprox,
bool hasBallparkVertical,
bool isNullTransformation)
: area_(area), accuracy_(accuracy), isPROJExportable_(isPROJExportable),
hasGrids_(hasGrids), gridsAvailable_(gridsAvailable),
gridsKnown_(gridsKnown), stepCount_(stepCount), isApprox_(isApprox),
hasBallparkVertical_(hasBallparkVertical),
isNullTransformation_(isNullTransformation) {}
};

Expand Down Expand Up @@ -11388,6 +11397,15 @@ struct SortFunction {
return false;
}

if (!iterA->second.hasBallparkVertical_ &&
iterB->second.hasBallparkVertical_) {
return true;
}
if (iterA->second.hasBallparkVertical_ &&
!iterB->second.hasBallparkVertical_) {
return false;
}

if (!iterA->second.isNullTransformation_ &&
iterB->second.isNullTransformation_) {
return true;
Expand Down Expand Up @@ -11654,7 +11672,9 @@ struct FilterResults {
? CoordinateOperationContext::SpatialCriterion::
STRICT_CONTAINMENT
: context->getSpatialCriterion();
bool hasFoundOpWithExtent = false;
bool hasOnlyBallpark = true;
bool hasNonBallparkWithoutExtent = false;
bool hasNonBallparkOpWithExtent = false;
const bool allowBallpark = context->getAllowBallparkTransformations();
for (const auto &op : sourceList) {
if (desiredAccuracy != 0) {
Expand All @@ -11669,9 +11689,15 @@ struct FilterResults {
if (areaOfInterest) {
bool emptyIntersection = false;
auto extent = getExtent(op, true, emptyIntersection);
if (!extent)
if (!extent) {
if (!op->hasBallparkTransformation()) {
hasNonBallparkWithoutExtent = true;
}
continue;
hasFoundOpWithExtent = true;
}
if (!op->hasBallparkTransformation()) {
hasNonBallparkOpWithExtent = true;
}
bool extentContains =
extent->contains(NN_NO_CHECK(areaOfInterest));
if (!hasOpThatContainsAreaOfInterestAndNoGrid &&
Expand All @@ -11698,9 +11724,15 @@ struct FilterResults {
BOTH) {
bool emptyIntersection = false;
auto extent = getExtent(op, true, emptyIntersection);
if (!extent)
if (!extent) {
if (!op->hasBallparkTransformation()) {
hasNonBallparkWithoutExtent = true;
}
continue;
hasFoundOpWithExtent = true;
}
if (!op->hasBallparkTransformation()) {
hasNonBallparkOpWithExtent = true;
}
bool extentContainsExtent1 =
!extent1 || extent->contains(NN_NO_CHECK(extent1));
bool extentContainsExtent2 =
Expand Down Expand Up @@ -11730,12 +11762,16 @@ struct FilterResults {
}
}
}
if (!op->hasBallparkTransformation()) {
hasOnlyBallpark = false;
}
res.emplace_back(op);
}

// In case no operation has an extent and no result is found,
// retain all initial operations that match accuracy criterion.
if (res.empty() && !hasFoundOpWithExtent) {
if ((res.empty() && !hasNonBallparkOpWithExtent) ||
(hasOnlyBallpark && hasNonBallparkWithoutExtent)) {
for (const auto &op : sourceList) {
if (desiredAccuracy != 0) {
const double accuracy = getAccuracy(op);
Expand Down Expand Up @@ -11839,6 +11875,8 @@ struct FilterResults {
area, getAccuracy(op), isPROJExportable, hasGrids,
gridsAvailable, gridsKnown, stepCount,
op->hasBallparkTransformation(),
op->nameStr().find("ballpark vertical transformation") !=
std::string::npos,
isNullTransformation(op->nameStr()));
}

Expand Down Expand Up @@ -13647,11 +13685,17 @@ bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
ENTER_FUNCTION();

if (geogSrc && vertDst) {
res = createOperationsGeogToVertFromGeoid(sourceCRS, targetCRS, vertDst,
context);
createOperationsFromDatabase(targetCRS, sourceCRS, context, geodDst,
geodSrc, geogDst, geogSrc, vertDst,
vertSrc, res);
res = applyInverse(res);
} else if (geogDst && vertSrc) {
res = applyInverse(createOperationsGeogToVertFromGeoid(
targetCRS, sourceCRS, vertSrc, context));
if (!res.empty()) {
createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context,
vertSrc, geogDst, res);
}
}

if (!res.empty()) {
Expand Down Expand Up @@ -14762,17 +14806,32 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog(
match.get(),
util::IComparable::Criterion::EQUIVALENT) &&
!match->identifiers().empty()) {
res = createOperations(
auto resTmp = createOperations(
NN_NO_CHECK(
util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
match)),
targetCRS, context);
res.insert(res.end(), resTmp.begin(), resTmp.end());
return;
}
}
}
}

createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context, vertSrc,
geogDst, res);
}

// ---------------------------------------------------------------------------

void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
Private::Context &, const crs::VerticalCRS *vertSrc,
const crs::GeographicCRS *geogDst,
std::vector<CoordinateOperationNNPtr> &res) {

ENTER_FUNCTION();

const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0];
const double convSrc = srcAxis->unit().conversionToSI();
double convDst = 1.0;
Expand All @@ -14791,12 +14850,24 @@ void CoordinateOperationFactory::Private::createOperationsVertToGeog(
((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp));

const double factor = convSrc / convDst;
auto conv = Transformation::createChangeVerticalUnit(
util::PropertyMap().set(
common::IdentifiedObject::NAME_KEY,

const auto &sourceCRSExtent = getExtent(sourceCRS);
const auto &targetCRSExtent = getExtent(targetCRS);
const bool sameExtent =
sourceCRSExtent && targetCRSExtent &&
sourceCRSExtent->_isEquivalentTo(
targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT);

util::PropertyMap map;
map.set(common::IdentifiedObject::NAME_KEY,
buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT),
sourceCRS, targetCRS,
BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT)
.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
sameExtent ? NN_NO_CHECK(sourceCRSExtent)
: metadata::Extent::WORLD);

auto conv = Transformation::createChangeVerticalUnit(
map, sourceCRS, targetCRS,
common::Scale(heightDepthReversal ? -factor : factor), {});
conv->setHasBallparkTransformation(true);
res.push_back(conv);
Expand Down
69 changes: 66 additions & 3 deletions test/cli/testprojinfo_out.dist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -990,7 +990,7 @@ PROJ.4 string:
+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs

Testing RH2000 height to SWEREF99: projinfo -s EPSG:5613 -t EPSG:4977
Candidate operations found: 1
Candidate operations found: 2
-------------------------------------
Operation No. 1:

Expand Down Expand Up @@ -1042,8 +1042,55 @@ COORDINATEOPERATION["Inverse of SWEREF99 to RH2000 height",
BBOX[55.28,10.93,69.07,24.17]],
ID["INVERSE(DERIVED_FROM(PROJ))","EPSG_4977_TO_EPSG_5613"]]

-------------------------------------
Operation No. 2:

unknown id, Transformation from RH2000 height to SWEREF99 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=noop

WKT2:2019 string:
COORDINATEOPERATION["Transformation from RH2000 height to SWEREF99 (ballpark vertical transformation, without ellipsoid height to vertical height correction)",
SOURCECRS[
VERTCRS["RH2000 height",
DYNAMIC[
FRAMEEPOCH[2000]],
VDATUM["Rikets hojdsystem 2000"],
CS[vertical,1],
AXIS["gravity-related height (H)",up,
LENGTHUNIT["metre",1]],
ID["EPSG",5613]]],
TARGETCRS[
GEOGCRS["SWEREF99",
DATUM["SWEREF99",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,3],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["ellipsoidal height (h)",up,
ORDER[3],
LENGTHUNIT["metre",1]],
ID["EPSG",4977]]],
METHOD["Change of Vertical Unit",
ID["EPSG",1069]],
PARAMETER["Unit conversion scalar",1,
SCALEUNIT["unity",1],
ID["EPSG",1051]],
USAGE[
SCOPE["unknown"],
AREA["World"],
BBOX[-90,-180,90,180]]]

Testing NAD83(2011) + NAVD88 height -> NAD83(2011) : projinfo -s EPSG:6349 -t EPSG:6319 --spatial-test intersects -o PROJ
Candidate operations found: 2
Candidate operations found: 3
-------------------------------------
Operation No. 1:

Expand All @@ -1070,8 +1117,16 @@ PROJ string:
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 3:

unknown id, Transformation from NAVD88 height to NAD83(2011) (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=noop

Testing NGF IGN69 height to RGF93: projinfo -s EPSG:5720 -t EPSG:4965 -o PROJ
Candidate operations found: 1
Candidate operations found: 2
-------------------------------------
Operation No. 1:

Expand All @@ -1085,6 +1140,14 @@ PROJ string:
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1

-------------------------------------
Operation No. 2:

unknown id, Transformation from NGF-IGN69 height to RGF93 (ballpark vertical transformation, without ellipsoid height to vertical height correction), unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=noop

Testing EPSG:32631 --3d
PROJ.4 string:
+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +type=crs
Expand Down
49 changes: 37 additions & 12 deletions test/unit/test_c_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4650,10 +4650,21 @@ TEST_F(CApi, proj_create_vertical_crs_ex) {
ObjectKeeper keeper_geog_crs(geog_crs);
ASSERT_NE(geog_crs, nullptr);

auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
nullptr);
ObjectKeeper keeper_P(P);
ASSERT_NE(P, nullptr);
PJ_OPERATION_FACTORY_CONTEXT *ctxt =
proj_create_operation_factory_context(m_ctxt, nullptr);
ASSERT_NE(ctxt, nullptr);
ContextKeeper keeper_ctxt(ctxt);
proj_operation_factory_context_set_grid_availability_use(
m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED);
proj_operation_factory_context_set_spatial_criterion(
m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
PJ_OBJ_LIST *operations =
proj_create_operations(m_ctxt, compound, geog_crs, ctxt);
ASSERT_NE(operations, nullptr);
ObjListKeeper keeper_operations(operations);
EXPECT_GE(proj_list_get_count(operations), 1);
auto P = proj_list_get(m_ctxt, operations, 0);
ObjectKeeper keeper_transform(P);

auto name = proj_get_name(P);
ASSERT_TRUE(name != nullptr);
Expand Down Expand Up @@ -4706,10 +4717,21 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
ObjectKeeper keeper_geog_crs(geog_crs);
ASSERT_NE(geog_crs, nullptr);

auto P = proj_create_crs_to_crs_from_pj(m_ctxt, compound, geog_crs, nullptr,
nullptr);
ObjectKeeper keeper_P(P);
ASSERT_NE(P, nullptr);
PJ_OPERATION_FACTORY_CONTEXT *ctxt =
proj_create_operation_factory_context(m_ctxt, nullptr);
ASSERT_NE(ctxt, nullptr);
ContextKeeper keeper_ctxt(ctxt);
proj_operation_factory_context_set_grid_availability_use(
m_ctxt, ctxt, PROJ_GRID_AVAILABILITY_IGNORED);
proj_operation_factory_context_set_spatial_criterion(
m_ctxt, ctxt, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
PJ_OBJ_LIST *operations =
proj_create_operations(m_ctxt, compound, geog_crs, ctxt);
ASSERT_NE(operations, nullptr);
ObjListKeeper keeper_operations(operations);
EXPECT_GE(proj_list_get_count(operations), 1);
auto P = proj_list_get(m_ctxt, operations, 0);
ObjectKeeper keeper_transform(P);

auto name = proj_get_name(P);
ASSERT_TRUE(name != nullptr);
Expand Down Expand Up @@ -4739,10 +4761,13 @@ TEST_F(CApi, proj_create_vertical_crs_ex_with_geog_crs) {
ObjectKeeper keeper_compound_from_projjson(compound_from_projjson);
ASSERT_NE(compound_from_projjson, nullptr);

auto P2 = proj_create_crs_to_crs_from_pj(m_ctxt, compound_from_projjson,
geog_crs, nullptr, nullptr);
ObjectKeeper keeper_P2(P2);
ASSERT_NE(P2, nullptr);
PJ_OBJ_LIST *operations2 =
proj_create_operations(m_ctxt, compound_from_projjson, geog_crs, ctxt);
ASSERT_NE(operations2, nullptr);
ObjListKeeper keeper_operations2(operations2);
EXPECT_GE(proj_list_get_count(operations2), 1);
auto P2 = proj_list_get(m_ctxt, operations2, 0);
ObjectKeeper keeper_transform2(P2);

auto name_bis = proj_get_name(P2);
ASSERT_TRUE(name_bis != nullptr);
Expand Down
Loading

0 comments on commit 5172091

Please sign in to comment.