From d8fe9964bcbc6d1eeaddf6f5d47ca6d444dc8744 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 22 Nov 2020 15:08:41 +0100 Subject: Add +proj=topocentric geocentric->topocentric conversion (fixes #500) --- .../ECEF_ENU_Longitude_Latitude_relationships.png | Bin 0 -> 30731 bytes .../ECEF_ENU_Longitude_Latitude_relationships.svg | 348 +++++++++++++++++++++ docs/source/operations/conversions/index.rst | 1 + docs/source/operations/conversions/topocentric.rst | 115 +++++++ .../internal/coordinateoperation_constants.hpp | 41 +++ src/Makefile.am | 1 + src/conversions/topocentric.cpp | 165 ++++++++++ src/iso19111/coordinateoperation.cpp | 55 +++- src/lib_proj.cmake | 1 + src/pj_list.h | 1 + src/proj_constants.h | 20 ++ test/gie/builtins.gie | 42 +++ test/unit/test_operation.cpp | 122 ++++++++ 13 files changed, 904 insertions(+), 8 deletions(-) create mode 100644 docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.png create mode 100644 docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.svg create mode 100644 docs/source/operations/conversions/topocentric.rst create mode 100644 src/conversions/topocentric.cpp diff --git a/docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.png b/docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.png new file mode 100644 index 00000000..224ff738 Binary files /dev/null and b/docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.png differ diff --git a/docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.svg b/docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.svg new file mode 100644 index 00000000..c97ca181 --- /dev/null +++ b/docs/source/operations/conversions/images/ECEF_ENU_Longitude_Latitude_relationships.svg @@ -0,0 +1,348 @@ + + + +Z + + + + + + +Y + + + + + + +X + + + + + + +North + + + + + + +East + + + + + + +Up + + + + + + +ecef + + +ecef + + +ecef + + + +φ + +λ + + \ No newline at end of file diff --git a/docs/source/operations/conversions/index.rst b/docs/source/operations/conversions/index.rst index b0bd0d37..b335c2ac 100644 --- a/docs/source/operations/conversions/index.rst +++ b/docs/source/operations/conversions/index.rst @@ -19,4 +19,5 @@ conversions. pop push set + topocentric unitconvert diff --git a/docs/source/operations/conversions/topocentric.rst b/docs/source/operations/conversions/topocentric.rst new file mode 100644 index 00000000..828ef4e0 --- /dev/null +++ b/docs/source/operations/conversions/topocentric.rst @@ -0,0 +1,115 @@ +.. _topocentric: + +================================================================================ +Geocentric to topocentric conversion +================================================================================ + +.. versionadded:: 8.0.0 + +Convert geocentric coordinates to topocentric coordinates (in the forward path). + ++---------------------+--------------------------------------------------------+ +| **Alias** | topocentric | ++---------------------+--------------------------------------------------------+ +| **Domain** | 3D | ++---------------------+--------------------------------------------------------+ +| **Input type** | Geocentric cartesian coordinates | ++---------------------+--------------------------------------------------------+ +| **Output type** | Topocentric cartesian coordinates | ++---------------------+--------------------------------------------------------+ + +This operation converts geocentric coordinate values (X, Y, Z) to topocentric +(E/East, N/North, U/Up) values. This is also sometimes known as the ECEF (Earth +Centered Earth Fixed) to ENU conversion. + +Topocentric coordinates are expressed in a frame whose East and North axis form +a local tangent plane to the Earth's ellipsoidal surface fixed to a specific +location (the topocentric origin), and the Up axis points upwards along the +normal to that plane. + +.. image:: ./images/ECEF_ENU_Longitude_Latitude_relationships.png + :align: center + :scale: 100% + :alt: ENU coordinate frame + +.. + Source : https://en.wikipedia.org/wiki/Local_tangent_plane_coordinates#/media/File:ECEF_ENU_Longitude_Latitude_relationships.svg + Public domain +.. + +The topocentric origin is a required parameter of the conversion, and can be +expressed either as geocentric coordinates (``X_0``, ``Y_0`` and ``Z_0``) or +as geographic coordinates (``lat_0``, ``lon_0``, ``h_0``). + +When conversion between geographic and topocentric coordinates is desired, the +topocentric conversion must be preceded by the :ref:`cart` conversion to +perform the initial geographic to geocentric coordinates conversion. + +The formulas used come from the "Geocentric/topocentric conversions" paragraph +of :cite:`IOGP2018` + +Usage +################################################################################ + +Convert geocentric coordinates to topocentric coordinates, with the topocentric +origin specified in geocentric coordinates:: + + echo 3771793.968 140253.342 5124304.349 2020 | \ + cct -d 3 +proj=topocentric +ellps=WGS84 +X_0=3652755.3058 +Y_0=319574.6799 +Z_0=5201547.3536 + + -189013.869 -128642.040 -4220.171 2020.0000 + +Convert geographic coordinates to topocentric coordinates, with the topocentric +origin specified in geographic coordinates:: + + echo 2.12955 53.80939444 73 2020 | cct -d 3 +proj=pipeline \ + +step +proj=cart +ellps=WGS84 \ + +step +proj=topocentric +ellps=WGS84 +lon_0=5 +lat_0=55 +h_0=200 + + -189013.869 -128642.040 -4220.171 2020.0000 + + +Parameters +################################################################################ + +.. include:: ../options/ellps.rst + +Topocentric origin described as geocentric coordinates +------------------------------------------------------ + +.. note:: + + The below options are mutually exclusive with the ones to express the origin as geographic coordinates. + +.. option:: +X_0= + + Geocentric X value of the topocentric origin (in metre) + +.. option:: +Y_0= + + Geocentric Y value of the topocentric origin (in metre) + +.. option:: +Z_0= + + Geocentric Z value of the topocentric origin (in metre) + +Topocentric origin described as geographic coordinates +------------------------------------------------------ + +.. note:: + + The below options are mutually exclusive with the ones to express the origin as geocentric coordinates. + +.. option:: +lat_0= + + Latitude of topocentric origin (in degree) + +.. option:: +lon_0= + + Longitude of topocentric origin (in degree) + +.. option:: +h_0= + + Ellipsoidal height of topocentric origin (in metre) + + *Defaults to 0.0.* diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp index 592f6263..0ed3a027 100644 --- a/include/proj/internal/coordinateoperation_constants.hpp +++ b/include/proj/internal/coordinateoperation_constants.hpp @@ -532,6 +532,34 @@ static const ParamMapping *const paramsColombiaUrban[] = { ¶mProjectionPlaneOriginHeight, nullptr}; +static const ParamMapping paramGeocentricXTopocentricOrigin = { + EPSG_NAME_PARAMETER_GEOCENTRIC_X_TOPOCENTRIC_ORIGIN, + EPSG_CODE_PARAMETER_GEOCENTRIC_X_TOPOCENTRIC_ORIGIN, nullptr, + common::UnitOfMeasure::Type::LINEAR, "X_0"}; + +static const ParamMapping paramGeocentricYTopocentricOrigin = { + EPSG_NAME_PARAMETER_GEOCENTRIC_Y_TOPOCENTRIC_ORIGIN, + EPSG_CODE_PARAMETER_GEOCENTRIC_Y_TOPOCENTRIC_ORIGIN, nullptr, + common::UnitOfMeasure::Type::LINEAR, "Y_0"}; + +static const ParamMapping paramGeocentricZTopocentricOrigin = { + EPSG_NAME_PARAMETER_GEOCENTRIC_Z_TOPOCENTRIC_ORIGIN, + EPSG_CODE_PARAMETER_GEOCENTRIC_Z_TOPOCENTRIC_ORIGIN, nullptr, + common::UnitOfMeasure::Type::LINEAR, "Z_0"}; + +static const ParamMapping *const paramsGeocentricTopocentric[] = { + ¶mGeocentricXTopocentricOrigin, ¶mGeocentricYTopocentricOrigin, + ¶mGeocentricZTopocentricOrigin, nullptr}; + +static const ParamMapping paramHeightTopoOriginWithH0 = { + EPSG_NAME_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN, + EPSG_CODE_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN, nullptr, + common::UnitOfMeasure::Type::LINEAR, "h_0"}; + +static const ParamMapping *const paramsGeographicTopocentric[] = { + ¶mLatTopoOrigin, ¶mLonTopoOrigin, ¶mHeightTopoOriginWithH0, + nullptr}; + static const MethodMapping projectionMethodMappings[] = { {EPSG_NAME_METHOD_TRANSVERSE_MERCATOR, EPSG_CODE_METHOD_TRANSVERSE_MERCATOR, "Transverse_Mercator", "tmerc", nullptr, paramsNatOriginScaleK}, @@ -836,6 +864,14 @@ static const MethodMapping projectionMethodMappings[] = { {EPSG_NAME_METHOD_COLOMBIA_URBAN, EPSG_CODE_METHOD_COLOMBIA_URBAN, nullptr, "col_urban", nullptr, paramsColombiaUrban}, + + {EPSG_NAME_METHOD_GEOCENTRIC_TOPOCENTRIC, + EPSG_CODE_METHOD_GEOCENTRIC_TOPOCENTRIC, nullptr, "topocentric", nullptr, + paramsGeocentricTopocentric}, + + {EPSG_NAME_METHOD_GEOGRAPHIC_TOPOCENTRIC, + EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC, nullptr, nullptr, nullptr, + paramsGeographicTopocentric}, }; #define METHOD_NAME_CODE(method) \ @@ -878,6 +914,8 @@ static const struct MethodNameCode { METHOD_NAME_CODE(AXIS_ORDER_REVERSAL_2D), METHOD_NAME_CODE(AXIS_ORDER_REVERSAL_3D), METHOD_NAME_CODE(GEOGRAPHIC_GEOCENTRIC), + METHOD_NAME_CODE(GEOCENTRIC_TOPOCENTRIC), + METHOD_NAME_CODE(GEOGRAPHIC_TOPOCENTRIC), // Transformations METHOD_NAME_CODE(LONGITUDE_ROTATION), METHOD_NAME_CODE(AFFINE_PARAMETRIC_TRANSFORMATION), @@ -943,6 +981,9 @@ static const struct ParamNameCode { PARAM_NAME_CODE(LONGITUDE_OF_ORIGIN), PARAM_NAME_CODE(ELLIPSOID_SCALE_FACTOR), PARAM_NAME_CODE(PROJECTION_PLANE_ORIGIN_HEIGHT), + PARAM_NAME_CODE(GEOCENTRIC_X_TOPOCENTRIC_ORIGIN), + PARAM_NAME_CODE(GEOCENTRIC_Y_TOPOCENTRIC_ORIGIN), + PARAM_NAME_CODE(GEOCENTRIC_Z_TOPOCENTRIC_ORIGIN), // Parameters of transformations PARAM_NAME_CODE(SEMI_MAJOR_AXIS_DIFFERENCE), PARAM_NAME_CODE(FLATTENING_DIFFERENCE), diff --git a/src/Makefile.am b/src/Makefile.am index 5b36c8bd..4b904d9f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -187,6 +187,7 @@ libproj_la_SOURCES = \ conversions/geoc.cpp \ conversions/geocent.cpp \ conversions/noop.cpp \ + conversions/topocentric.cpp \ conversions/set.cpp \ conversions/unitconvert.cpp \ \ diff --git a/src/conversions/topocentric.cpp b/src/conversions/topocentric.cpp new file mode 100644 index 00000000..214f8221 --- /dev/null +++ b/src/conversions/topocentric.cpp @@ -0,0 +1,165 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Convert between geocentric coordinates and topocentric (ENU) coordinates + * + ****************************************************************************** + * Copyright (c) 2020, Even Rouault + * + * 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. + *****************************************************************************/ + +#define PJ_LIB__ + +#include "proj_internal.h" +#include +#include + +PROJ_HEAD(topocentric, "Geocentric/Topocentric conversion"); + +// Notations and formulas taken from IOGP Publication 373-7-2 - +// Geomatics Guidance Note number 7, part 2 - October 2020 + +namespace { // anonymous namespace +struct pj_opaque { + double X0; + double Y0; + double Z0; + double sinphi0; + double cosphi0; + double sinlam0; + double coslam0; +}; +} // anonymous namespace + +// Convert from geocentric to topocentric +static PJ_COORD topocentric_fwd(PJ_COORD in, PJ * P) +{ + struct pj_opaque *Q = static_cast(P->opaque); + PJ_COORD out; + const double dX = in.xyz.x - Q->X0; + const double dY = in.xyz.y - Q->Y0; + const double dZ = in.xyz.z - Q->Z0; + out.xyz.x = -dX * Q->sinlam0 + dY * Q->coslam0; + out.xyz.y = -dX * Q->sinphi0 * Q->coslam0 - dY * Q->sinphi0 * Q->sinlam0 + dZ * Q->cosphi0; + out.xyz.z = dX * Q->cosphi0 * Q->coslam0 + dY * Q->cosphi0 * Q->sinlam0 + dZ * Q->sinphi0; + return out; +} + +// Convert from topocentric to geocentric +static PJ_COORD topocentric_inv(PJ_COORD in, PJ * P) +{ + struct pj_opaque *Q = static_cast(P->opaque); + PJ_COORD out; + out.xyz.x = Q->X0 - in.xyz.x * Q->sinlam0 - in.xyz.y * Q->sinphi0 * Q->coslam0 + in.xyz.z * Q->cosphi0 * Q->coslam0; + out.xyz.y = Q->Y0 + in.xyz.x * Q->coslam0 - in.xyz.y * Q->sinphi0 * Q->sinlam0 + in.xyz.z * Q->cosphi0 * Q->sinlam0; + out.xyz.z = Q->Z0 + in.xyz.y * Q->cosphi0 + in.xyz.z * Q->sinphi0; + return out; +} + + +/*********************************************************************/ +PJ *CONVERSION(topocentric,1) { +/*********************************************************************/ + struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = static_cast(Q); + + // The topocentric origin can be specified either in geocentric coordinates + // (X_0,Y_0,Z_0) or as geographic coordinates (lon_0,lat_0,h_0) + // Checks: + // - X_0 or lon_0 must be specified + // - If X_0 is specified, the Y_0 and Z_0 must also be + // - If lon_0 is specified, then lat_0 must also be + // - If any of X_0, Y_0, Z_0 is specified, then any of lon_0,lat_0,h_0 must + // not be, and vice versa. + const auto hasX0 = pj_param_exists(P->params, "X_0"); + const auto hasY0 = pj_param_exists(P->params, "Y_0"); + const auto hasZ0 = pj_param_exists(P->params, "Z_0"); + const auto hasLon0 = pj_param_exists(P->params, "lon_0"); + const auto hasLat0 = pj_param_exists(P->params, "lat_0"); + const auto hash0 = pj_param_exists(P->params, "h_0"); + if( !hasX0 && !hasLon0 ) + { + return pj_default_destructor(P, PJD_ERR_MISSING_ARGS); + } + if ( (hasX0 || hasY0 || hasZ0) && + (hasLon0 || hasLat0 || hash0) ) + { + return pj_default_destructor(P, PJD_ERR_MUTUALLY_EXCLUSIVE_ARGS); + } + if( hasX0 && (!hasY0 || !hasZ0) ) + { + return pj_default_destructor(P, PJD_ERR_MISSING_ARGS); + } + if( hasLon0 && !hasLat0 ) // allow missing h_0 + { + return pj_default_destructor(P, PJD_ERR_MISSING_ARGS); + } + + // Pass a dummy ellipsoid definition that will be overridden just afterwards + PJ* cart = proj_create(P->ctx, "+proj=cart +a=1"); + if (cart == nullptr) + return pj_default_destructor(P, ENOMEM); + /* inherit ellipsoid definition from P to cart */ + pj_inherit_ellipsoid_def (P, cart); + + if( hasX0 ) + { + Q->X0 = pj_param(P->ctx, P->params, "dX_0").f; + Q->Y0 = pj_param(P->ctx, P->params, "dY_0").f; + Q->Z0 = pj_param(P->ctx, P->params, "dZ_0").f; + + // Compute lam0, phi0 from X0,Y0,Z0 + PJ_XYZ xyz; + xyz.x = Q->X0; + xyz.y = Q->Y0; + xyz.z = Q->Z0; + const auto lpz = pj_inv3d(xyz, cart); + Q->sinphi0 = sin(lpz.phi); + Q->cosphi0 = cos(lpz.phi); + Q->sinlam0 = sin(lpz.lam); + Q->coslam0 = cos(lpz.lam); + } + else + { + // Compute X0,Y0,Z0 from lam0, phi0, h0 + PJ_LPZ lpz; + lpz.lam = P->lam0; + lpz.phi = P->phi0; + lpz.z = pj_param(P->ctx, P->params, "dh_0").f; + const auto xyz = pj_fwd3d(lpz, cart); + Q->X0 = xyz.x; + Q->Y0 = xyz.y; + Q->Z0 = xyz.z; + + Q->sinphi0 = sin(P->phi0); + Q->cosphi0 = cos(P->phi0); + Q->sinlam0 = sin(P->lam0); + Q->coslam0 = cos(P->lam0); + } + + proj_destroy(cart); + + P->fwd4d = topocentric_fwd; + P->inv4d = topocentric_inv; + P->left = PJ_IO_UNITS_CARTESIAN; + P->right = PJ_IO_UNITS_CARTESIAN; + return P; +} diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 1d3b7a90..dfd349b5 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -6034,28 +6034,39 @@ void Conversion::_exportToPROJString( !isHeightDepthReversal; bool applyTargetCRSModifiers = applySourceCRSModifiers; + if (formatter->getCRSExport()) { + if (methodEPSGCode == EPSG_CODE_METHOD_GEOCENTRIC_TOPOCENTRIC || + methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) { + throw io::FormattingException("Transformation cannot be exported " + "as a PROJ.4 string (but can be part " + "of a PROJ pipeline)"); + } + } + auto l_sourceCRS = sourceCRS(); + crs::GeographicCRSPtr srcGeogCRS; if (!formatter->getCRSExport() && l_sourceCRS && applySourceCRSModifiers) { - crs::CRS *horiz = l_sourceCRS.get(); - const auto compound = dynamic_cast(horiz); + crs::CRSPtr horiz = l_sourceCRS; + const auto compound = + dynamic_cast(l_sourceCRS.get()); if (compound) { const auto &components = compound->componentReferenceSystems(); if (!components.empty()) { - horiz = components.front().get(); + horiz = components.front().as_nullable(); } } - auto geogCRS = dynamic_cast(horiz); - if (geogCRS) { + srcGeogCRS = std::dynamic_pointer_cast(horiz); + if (srcGeogCRS) { formatter->setOmitProjLongLatIfPossible(true); formatter->startInversion(); - geogCRS->_exportToPROJString(formatter); + srcGeogCRS->_exportToPROJString(formatter); formatter->stopInversion(); formatter->setOmitProjLongLatIfPossible(false); } - auto projCRS = dynamic_cast(horiz); + auto projCRS = dynamic_cast(horiz.get()); if (projCRS) { formatter->startInversion(); formatter->pushOmitZUnitConversion(); @@ -6325,6 +6336,30 @@ void Conversion::_exportToPROJString( } bConversionDone = true; bEllipsoidParametersDone = true; + } else if (methodEPSGCode == EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC) { + if (!srcGeogCRS) { + throw io::FormattingException( + "Export of Geographic/Topocentric conversion to a PROJ string " + "requires an input geographic CRS"); + } + + formatter->addStep("cart"); + srcGeogCRS->ellipsoid()->_exportToPROJString(formatter); + + formatter->addStep("topocentric"); + const auto latOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_LATITUDE_TOPOGRAPHIC_ORIGIN, + common::UnitOfMeasure::DEGREE); + const auto lonOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_LONGITUDE_TOPOGRAPHIC_ORIGIN, + common::UnitOfMeasure::DEGREE); + const auto heightOrigin = parameterValueNumeric( + EPSG_CODE_PARAMETER_ELLIPSOIDAL_HEIGHT_TOPOCENTRIC_ORIGIN, + common::UnitOfMeasure::METRE); + formatter->addParam("lat_0", latOrigin); + formatter->addParam("lon_0", lonOrigin); + formatter->addParam("h_0", heightOrigin); + bConversionDone = true; } auto l_targetCRS = targetCRS(); @@ -6449,7 +6484,9 @@ void Conversion::_exportToPROJString( } if (!bEllipsoidParametersDone) { - auto targetGeogCRS = horiz->extractGeographicCRS(); + auto targetGeodCRS = horiz->extractGeodeticCRS(); + auto targetGeogCRS = + std::dynamic_pointer_cast(targetGeodCRS); if (targetGeogCRS) { if (formatter->getCRSExport()) { targetGeogCRS->addDatumInfoToPROJString(formatter); @@ -6458,6 +6495,8 @@ void Conversion::_exportToPROJString( targetGeogCRS->primeMeridian()->_exportToPROJString( formatter); } + } else if (targetGeodCRS) { + targetGeodCRS->ellipsoid()->_exportToPROJString(formatter); } } diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 67bc1f4e..9f482d19 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -173,6 +173,7 @@ set(SRC_LIBPROJ_CONVERSIONS conversions/geoc.cpp conversions/geocent.cpp conversions/noop.cpp + conversions/topocentric.cpp conversions/set.cpp conversions/unitconvert.cpp ) diff --git a/src/pj_list.h b/src/pj_list.h index bcdc189e..d00e780f 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -155,6 +155,7 @@ PROJ_HEAD(tinshift, "Triangulation based transformation") PROJ_HEAD(tissot, "Tissot Conic") PROJ_HEAD(tmerc, "Transverse Mercator") PROJ_HEAD(tobmerc, "Tobler-Mercator") +PROJ_HEAD(topocentric, "Geocentric/Topocentric conversion") PROJ_HEAD(tpeqd, "Two Point Equidistant") PROJ_HEAD(tpers, "Tilted perspective") PROJ_HEAD(unitconvert, "Unit conversion") diff --git a/src/proj_constants.h b/src/proj_constants.h index a3da2c10..ce3b2157 100644 --- a/src/proj_constants.h +++ b/src/proj_constants.h @@ -651,4 +651,24 @@ #define EPSG_CODE_METHOD_HEIGHT_DEPTH_REVERSAL 1068 #define EPSG_NAME_METHOD_HEIGHT_DEPTH_REVERSAL "Height Depth Reversal" +/* ------------------------------------------------------------------------ */ + +#define EPSG_NAME_METHOD_GEOCENTRIC_TOPOCENTRIC "Geocentric/topocentric conversions" +#define EPSG_CODE_METHOD_GEOCENTRIC_TOPOCENTRIC 9836 + +#define EPSG_NAME_PARAMETER_GEOCENTRIC_X_TOPOCENTRIC_ORIGIN "Geocentric X of topocentric origin" +#define EPSG_CODE_PARAMETER_GEOCENTRIC_X_TOPOCENTRIC_ORIGIN 8837 + +#define EPSG_NAME_PARAMETER_GEOCENTRIC_Y_TOPOCENTRIC_ORIGIN "Geocentric Y of topocentric origin" +#define EPSG_CODE_PARAMETER_GEOCENTRIC_Y_TOPOCENTRIC_ORIGIN 8838 + +#define EPSG_NAME_PARAMETER_GEOCENTRIC_Z_TOPOCENTRIC_ORIGIN "Geocentric Z of topocentric origin" +#define EPSG_CODE_PARAMETER_GEOCENTRIC_Z_TOPOCENTRIC_ORIGIN 8839 + +/* ------------------------------------------------------------------------ */ + +#define EPSG_NAME_METHOD_GEOGRAPHIC_TOPOCENTRIC "Geographic/topocentric conversions" +#define EPSG_CODE_METHOD_GEOGRAPHIC_TOPOCENTRIC 9837 + + #endif /* PROJ_CONSTANTS_INCLUDED */ diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index add5d925..72984931 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -6681,4 +6681,46 @@ accept -74.25 4.8 expect 80859.033 122543.174 roundtrip 1 +=============================================================================== +# Geocentric/topocentric conversion +=============================================================================== + +# Test parameter and point from IOGP Publication 373-7-2 - Geomatics Guidance Note number 7, part 2 - October 2020 +operation +proj=topocentric +ellps=WGS84 +X_0=3652755.3058 +Y_0=319574.6799 +Z_0=5201547.3536 +tolerance 1 mm +accept 3771793.968 140253.342 5124304.349 +expect -189013.869 -128642.040 -4220.171 +roundtrip 1 + +=============================================================================== +# Geographic/topocentric conversion +=============================================================================== + +# Test parameter and point from IOGP Publication 373-7-2 - Geomatics Guidance Note number 7, part 2 - October 2020 +operation +proj=pipeline +step +proj=cart +ellps=WGS84 +step +proj=topocentric +ellps=WGS84 +lon_0=5 +lat_0=55 +h_0=200 +tolerance 1 mm +accept 2.12955 53.80939444444444 73 +expect -189013.869 -128642.040 -4220.171 +roundtrip 1 + +=============================================================================== +# Error cases of topocentric +=============================================================================== + +# missing X_0,Y_0,Z_0 or lon_0,lat_0 +operation +proj=topocentric +ellps=WGS84 +expect failure errno missing_args + +# missing Z_0 +operation +proj=topocentric +ellps=WGS84 +X_0=0 +Y_0=0 +expect failure errno missing_args + +# missing lat_0 +operation +proj=topocentric +ellps=WGS84 +lon_0=0 +expect failure errno missing_args + +# X_0 and lon_0 are mutually exclusive +operation +proj=topocentric +ellps=WGS84 +X_0=0 +lon_0=0 +expect failure errno mutually_exclusive_args + diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index b6e555b1..269b1bbd 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -10098,6 +10098,128 @@ TEST(operation, derivedGeographicCRS_with_to_wgs84_to_geographicCRS) { // --------------------------------------------------------------------------- +TEST(operation, geographic_topocentric) { + auto wkt = + "PROJCRS[\"EPSG topocentric example A\",\n" + " BASEGEOGCRS[\"WGS 84\",\n" + " ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n" + " MEMBER[\"World Geodetic System 1984 (Transit)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G730)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G873)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1150)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1674)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1762)\"],\n" + " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ENSEMBLEACCURACY[2.0]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " ID[\"EPSG\",4979]],\n" + " CONVERSION[\"EPSG topocentric example A\",\n" + " METHOD[\"Geographic/topocentric conversions\",\n" + " ID[\"EPSG\",9837]],\n" + " PARAMETER[\"Latitude of topocentric origin\",55,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8834]],\n" + " PARAMETER[\"Longitude of topocentric origin\",5,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433],\n" + " ID[\"EPSG\",8835]],\n" + " PARAMETER[\"Ellipsoidal height of topocentric origin\",0,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8836]]],\n" + " CS[Cartesian,3],\n" + " AXIS[\"topocentric East (U)\",east,\n" + " ORDER[1],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " AXIS[\"topocentric North (V)\",north,\n" + " ORDER[2],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " AXIS[\"topocentric height (W)\",up,\n" + " ORDER[3],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " USAGE[\n" + " SCOPE[\"Example only (fictitious).\"],\n" + " AREA[\"Description of the extent of the CRS.\"],\n" + " BBOX[-90,-180,90,180]],\n" + " ID[\"EPSG\",5819]]"; + auto obj = WKTParser().createFromWKT(wkt); + auto dst = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(dst != nullptr); + auto op = CoordinateOperationFactory::create()->createOperation( + GeographicCRS::EPSG_4979, NN_CHECK_ASSERT(dst)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + "+step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m " + "+step +proj=cart +ellps=WGS84 " + "+step +proj=topocentric +lat_0=55 +lon_0=5 +h_0=0 +ellps=WGS84"); +} + +// --------------------------------------------------------------------------- + +TEST(operation, geocentric_topocentric) { + auto wkt = + "PROJCRS[\"EPSG topocentric example B\",\n" + " BASEGEODCRS[\"WGS 84\",\n" + " ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n" + " MEMBER[\"World Geodetic System 1984 (Transit)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G730)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G873)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1150)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1674)\"],\n" + " MEMBER[\"World Geodetic System 1984 (G1762)\"],\n" + " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n" + " LENGTHUNIT[\"metre\",1]],\n" + " ENSEMBLEACCURACY[2.0]],\n" + " PRIMEM[\"Greenwich\",0,\n" + " ANGLEUNIT[\"degree\",0.0174532925199433]],\n" + " ID[\"EPSG\",4978]],\n" + " CONVERSION[\"EPSG topocentric example B\",\n" + " METHOD[\"Geocentric/topocentric conversions\",\n" + " ID[\"EPSG\",9836]],\n" + " PARAMETER[\"Geocentric X of topocentric origin\",3771793.97,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8837]],\n" + " PARAMETER[\"Geocentric Y of topocentric origin\",140253.34,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8838]],\n" + " PARAMETER[\"Geocentric Z of topocentric origin\",5124304.35,\n" + " LENGTHUNIT[\"metre\",1],\n" + " ID[\"EPSG\",8839]]],\n" + " CS[Cartesian,3],\n" + " AXIS[\"topocentric East (U)\",east,\n" + " ORDER[1],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " AXIS[\"topocentric North (V)\",north,\n" + " ORDER[2],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " AXIS[\"topocentric height (W)\",up,\n" + " ORDER[3],\n" + " LENGTHUNIT[\"metre\",1]],\n" + " USAGE[\n" + " SCOPE[\"Example only (fictitious).\"],\n" + " AREA[\"Description of the extent of the CRS.\"],\n" + " BBOX[-90,-180,90,180]],\n" + " ID[\"EPSG\",5820]]"; + auto dbContext = DatabaseContext::create(); + // Need a database so that EPSG:4978 is resolved + auto obj = WKTParser().attachDatabaseContext(dbContext).createFromWKT(wkt); + auto dst = nn_dynamic_pointer_cast(obj); + ASSERT_TRUE(dst != nullptr); + auto f(NS_PROJ::io::WKTFormatter::create( + NS_PROJ::io::WKTFormatter::Convention::WKT2_2019)); + std::cerr << GeodeticCRS::EPSG_4978->exportToWKT(f.get()) << std::endl; + auto op = CoordinateOperationFactory::create()->createOperation( + GeodeticCRS::EPSG_4978, NN_CHECK_ASSERT(dst)); + ASSERT_TRUE(op != nullptr); + EXPECT_EQ(op->exportToPROJString(PROJStringFormatter::create().get()), + "+proj=topocentric +X_0=3771793.97 +Y_0=140253.34 " + "+Z_0=5124304.35 +ellps=WGS84"); +} + +// --------------------------------------------------------------------------- + TEST(operation, mercator_variant_A_to_variant_B) { auto projCRS = ProjectedCRS::create( PropertyMap(), GeographicCRS::EPSG_4326, -- cgit v1.2.3