From 076e6a82b7aa84ebc6d092e5f808d1daf3c00393 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 25 Oct 2020 18:06:57 +0100 Subject: [PATCH] Add +proj=colombia_urban projection, implementing a EPSG projection method used by a number of projected CRS in Colombia (fixes #589) --- .../operations/projections/colombia_urban.rst | 49 ++++++++++++ docs/source/operations/projections/index.rst | 1 + .../coordinateoperation_constants.hpp | 19 ++++- src/Makefile.am | 1 + src/lib_proj.cmake | 1 + src/pj_list.h | 1 + src/proj_constants.h | 6 ++ src/projections/colombia_urban.cpp | 76 +++++++++++++++++++ test/cli/testvarious | 6 ++ test/cli/tv_out.dist | 3 + test/gie/builtins.gie | 14 ++++ 11 files changed, 176 insertions(+), 1 deletion(-) create mode 100644 docs/source/operations/projections/colombia_urban.rst create mode 100644 src/projections/colombia_urban.cpp diff --git a/docs/source/operations/projections/colombia_urban.rst b/docs/source/operations/projections/colombia_urban.rst new file mode 100644 index 0000000000..169e741625 --- /dev/null +++ b/docs/source/operations/projections/colombia_urban.rst @@ -0,0 +1,49 @@ +.. _colombia_urban: + +******************************************************************************** +Colombia Urban +******************************************************************************** + +.. versionadded:: 7.2 + ++---------------------+----------------------------------------------------------+ +| **Classification** | Miscellaneous | ++---------------------+----------------------------------------------------------+ +| **Available forms** | Forward and inverse ellipsoidal projection | ++---------------------+----------------------------------------------------------+ +| **Alias** | colombia_urban | ++---------------------+----------------------------------------------------------+ +| **Domain** | 2D | ++---------------------+----------------------------------------------------------+ +| **Input type** | Geodetic coordinates | ++---------------------+----------------------------------------------------------+ +| **Output type** | Projected coordinates | ++---------------------+----------------------------------------------------------+ + +From :cite:`IOGP2018`: + + The capital cites of each department in Colombia use an urban projection for large + scale topographic mapping of the urban areas. It is based on a plane through + the origin at an average height for the area, eliminating the need for corrections + to engineering survey observations. + + proj-string: ``+proj=colombia_urban`` + +Parameters +################################################################################ + +.. include:: ../options/lon_0.rst + +.. include:: ../options/lat_0.rst + +.. include:: ../options/ellps.rst + +.. include:: ../options/x_0.rst + +.. include:: ../options/y_0.rst + +.. option:: +h_0= + + Projection plane origin height (in metre) + + *Defaults to 0.0.* diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst index 169187cc67..0ce7d111c9 100644 --- a/docs/source/operations/projections/index.rst +++ b/docs/source/operations/projections/index.rst @@ -34,6 +34,7 @@ Projections map the spherical 3D space to a flat 2D space. cea chamb collg + colombia_urban comill crast denoy diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp index 7038ba9fe5..bf703fe0e0 100644 --- a/include/proj/internal/coordinateoperation_constants.hpp +++ b/include/proj/internal/coordinateoperation_constants.hpp @@ -519,6 +519,19 @@ static const ParamMapping *const paramsVerticalPerspective[] = { ¶mFalseNorthing, // PROJ addition nullptr}; +static const ParamMapping paramProjectionPlaneOriginHeight = { + EPSG_NAME_PARAMETER_PROJECTION_PLANE_ORIGIN_HEIGHT, + EPSG_CODE_PARAMETER_PROJECTION_PLANE_ORIGIN_HEIGHT, nullptr, + common::UnitOfMeasure::Type::LINEAR, "h_0"}; + +static const ParamMapping *const paramsColombiaUrban[] = { + ¶mLatitudeNatOrigin, + ¶mLongitudeNatOrigin, + ¶mFalseEasting, + ¶mFalseNorthing, + ¶mProjectionPlaneOriginHeight, + nullptr}; + static const MethodMapping projectionMethodMappings[] = { {EPSG_NAME_METHOD_TRANSVERSE_MERCATOR, EPSG_CODE_METHOD_TRANSVERSE_MERCATOR, "Transverse_Mercator", "tmerc", nullptr, paramsNatOriginScaleK}, @@ -820,6 +833,9 @@ static const MethodMapping projectionMethodMappings[] = { {EPSG_NAME_METHOD_VERTICAL_PERSPECTIVE, EPSG_CODE_METHOD_VERTICAL_PERSPECTIVE, nullptr, "nsper", nullptr, paramsVerticalPerspective}, + + {EPSG_NAME_METHOD_COLOMBIA_URBAN, EPSG_CODE_METHOD_COLOMBIA_URBAN, nullptr, + "colombia_urban", nullptr, paramsColombiaUrban}, }; #define METHOD_NAME_CODE(method) \ @@ -855,7 +871,7 @@ static const struct MethodNameCode { METHOD_NAME_CODE(POLAR_STEREOGRAPHIC_VARIANT_A), METHOD_NAME_CODE(POLAR_STEREOGRAPHIC_VARIANT_B), METHOD_NAME_CODE(EQUAL_EARTH), METHOD_NAME_CODE(LABORDE_OBLIQUE_MERCATOR), - METHOD_NAME_CODE(VERTICAL_PERSPECTIVE), + METHOD_NAME_CODE(VERTICAL_PERSPECTIVE), METHOD_NAME_CODE(COLOMBIA_URBAN), // Other conversions METHOD_NAME_CODE(CHANGE_VERTICAL_UNIT), METHOD_NAME_CODE(HEIGHT_DEPTH_REVERSAL), @@ -926,6 +942,7 @@ static const struct ParamNameCode { PARAM_NAME_CODE(LATITUDE_STD_PARALLEL), PARAM_NAME_CODE(LONGITUDE_OF_ORIGIN), PARAM_NAME_CODE(ELLIPSOID_SCALE_FACTOR), + PARAM_NAME_CODE(PROJECTION_PLANE_ORIGIN_HEIGHT), // 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 68aeadf2db..92f4a41a96 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -180,6 +180,7 @@ libproj_la_SOURCES = \ projections/natearth2.cpp \ projections/calcofi.cpp \ projections/eqearth.cpp \ + projections/colombia_urban.cpp \ \ conversions/axisswap.cpp \ conversions/cart.cpp \ diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 1cd7f421d9..b607573861 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -164,6 +164,7 @@ set(SRC_LIBPROJ_PROJECTIONS projections/natearth2.cpp projections/calcofi.cpp projections/eqearth.cpp + projections/colombia_urban.cpp ) set(SRC_LIBPROJ_CONVERSIONS diff --git a/src/pj_list.h b/src/pj_list.h index 5b91af9b7e..ec69f877f2 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -30,6 +30,7 @@ PROJ_HEAD(ccon, "Central Conic") PROJ_HEAD(cea, "Equal Area Cylindrical") PROJ_HEAD(chamb, "Chamberlin Trimetric") PROJ_HEAD(collg, "Collignon") +PROJ_HEAD(colombia_urban, "Colombia Urban") PROJ_HEAD(comill, "Compact Miller") PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)") PROJ_HEAD(defmodel, "Deformation model") diff --git a/src/proj_constants.h b/src/proj_constants.h index f64bf4965d..a3da2c108b 100644 --- a/src/proj_constants.h +++ b/src/proj_constants.h @@ -234,6 +234,9 @@ #define PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION "Pole rotation (GRIB convention)" +#define EPSG_CODE_METHOD_COLOMBIA_URBAN 1052 +#define EPSG_NAME_METHOD_COLOMBIA_URBAN "Colombia Urban" + /* ------------------------------------------------------------------------ */ /* Projection parameters */ @@ -335,6 +338,9 @@ #define EPSG_NAME_PARAMETER_VIEWPOINT_HEIGHT "Viewpoint height" #define EPSG_CODE_PARAMETER_VIEWPOINT_HEIGHT 8840 +#define EPSG_NAME_PARAMETER_PROJECTION_PLANE_ORIGIN_HEIGHT "Projection plane origin height" +#define EPSG_CODE_PARAMETER_PROJECTION_PLANE_ORIGIN_HEIGHT 1039 + /* ------------------------------------------------------------------------ */ /* Other conversions and transformations */ diff --git a/src/projections/colombia_urban.cpp b/src/projections/colombia_urban.cpp new file mode 100644 index 0000000000..4aeb198c59 --- /dev/null +++ b/src/projections/colombia_urban.cpp @@ -0,0 +1,76 @@ +#define PJ_LIB__ + +#include +#include + +#include "proj.h" +#include "proj_internal.h" + +PROJ_HEAD(colombia_urban, "Colombia Urban") + "\n\tMisc\n\th_0="; + +// Notations and formulas taken from IOGP Publication 373-7-2 - +// Geomatics Guidance Note number 7, part 2 - March 2020 + +namespace { // anonymous namespace + +struct pj_opaque { + double h0; // height of projection origin, divided by semi-major axis (a) + double rho0; // adimensional value, contrary to Guidance note 7.2 + double A; + double B; // adimensional value, contrary to Guidance note 7.2 + double C; + double D; // adimensional value, contrary to Guidance note 7.2 +}; +} // anonymous namespace + +static PJ_XY colombia_urban_forward (PJ_LP lp, PJ *P) { + PJ_XY xy; + struct pj_opaque *Q = static_cast(P->opaque); + + const double cosphi = cos(lp.phi); + const double sinphi = sin(lp.phi); + const double nu = 1. / sqrt(1 - P->es * sinphi * sinphi); + const double lam_nu_cosphi = lp.lam * nu * cosphi; + xy.x = Q->A * lam_nu_cosphi; + const double sinphi_m = sin(0.5 * (lp.phi + P->phi0)); + const double rho_m = (1 - P->es) / pow(1 - P->es * sinphi_m * sinphi_m, 1.5); + const double G = 1 + Q->h0 / rho_m; + xy.y = G * Q->rho0 * ((lp.phi - P->phi0) + Q->B * lam_nu_cosphi * lam_nu_cosphi); + + return xy; +} + +static PJ_LP colombia_urban_inverse (PJ_XY xy, PJ *P) { + PJ_LP lp; + struct pj_opaque *Q = static_cast(P->opaque); + + lp.phi = P->phi0 + xy.y / Q->D - Q->B * (xy.x / Q->C) * (xy.x / Q->C); + const double sinphi = sin(lp.phi); + const double nu = 1. / sqrt(1 - P->es * sinphi * sinphi); + lp.lam = xy.x / (Q->C * nu * cos(lp.phi)); + + return lp; +} + +PJ *PROJECTION(colombia_urban) { + struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + const double h0_unscaled = pj_param(P->ctx, P->params, "dh_0").f; + Q->h0 = h0_unscaled / P->a; + const double sinphi0 = sin(P->phi0); + const double nu0 = 1. / sqrt(1 - P->es * sinphi0 * sinphi0); + Q->A = 1 + Q->h0 / nu0; + Q->rho0 = (1 - P->es) / pow(1 - P->es * sinphi0 * sinphi0, 1.5); + Q->B = tan(P->phi0) / (2 * Q->rho0 * nu0); + Q->C = 1 + Q->h0; + Q->D = Q->rho0 * (1 + Q->h0 / (1 - P->es)); + + P->fwd = colombia_urban_forward; + P->inv = colombia_urban_inverse; + + return P; +} diff --git a/test/cli/testvarious b/test/cli/testvarious index 6a77729ec8..292ee31665 100755 --- a/test/cli/testvarious +++ b/test/cli/testvarious @@ -1003,6 +1003,12 @@ echo "2 49" > tmp.txt $EXE EPSG:4326 EPSG:4326 tmp.txt -E >> ${OUT} rm tmp.txt +echo "##############################################################" >> ${OUT} +echo "Test Colombia Urban" >> ${OUT} +$EXE -f %.3f EPSG:4686 EPSG:6247 -E >> ${OUT} <