Merge pull request #2395 from rouault/colombia_urban

Add +proj=col_urban projection, implementing a EPSG projection method used by a number of projected CRS in Colombia (fixes #589)

11 files changed, 176 insertions, 1 deletions

diff --git a/docs/source/operations/projections/col_urban.rst b/docs/source/operations/projections/col_urban.rst
new file mode 100644
index 00000000..df2873dd
--- /dev/null
+++ b/docs/source/operations/projections/col_urban.rst
@@ -0,0 +1,49 @@
+.. _col_urban:
+
+********************************************************************************
+Colombia Urban
+********************************************************************************
+
+.. versionadded:: 7.2
+
++---------------------+----------------------------------------------------------+
+| **Classification**  | Miscellaneous                                            |
++---------------------+----------------------------------------------------------+
+| **Available forms** | Forward and inverse ellipsoidal projection               |
++---------------------+----------------------------------------------------------+
+| **Alias**           | col_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=col_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=<value>
+
+    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 169187cc..48420971 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
+   col_urban
    comill
    crast
    denoy
diff --git a/include/proj/internal/coordinateoperation_constants.hpp b/include/proj/internal/coordinateoperation_constants.hpp
index 7038ba9f..592f6263 100644
--- a/include/proj/internal/coordinateoperation_constants.hpp
+++ b/include/proj/internal/coordinateoperation_constants.hpp
@@ -519,6 +519,19 @@ static const ParamMapping *const paramsVerticalPerspective[] = {
     &paramFalseNorthing, // 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[] = {
+    &paramLatitudeNatOrigin,
+    &paramLongitudeNatOrigin,
+    &paramFalseEasting,
+    &paramFalseNorthing,
+    &paramProjectionPlaneOriginHeight,
+    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,
+     "col_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 68aeadf2..33cc95a6 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/col_urban.cpp \
 	\
 	conversions/axisswap.cpp \
 	conversions/cart.cpp \
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 1cd7f421..67bc1f4e 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/col_urban.cpp
 )
 
 set(SRC_LIBPROJ_CONVERSIONS
diff --git a/src/pj_list.h b/src/pj_list.h
index 5b91af9b..bcdc189e 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(col_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 f64bf496..a3da2c10 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/col_urban.cpp b/src/projections/col_urban.cpp
new file mode 100644
index 00000000..5bc8407f
--- /dev/null
+++ b/src/projections/col_urban.cpp
@@ -0,0 +1,76 @@
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj.h"
+#include "proj_internal.h"
+
+PROJ_HEAD(col_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 col_urban_forward (PJ_LP lp, PJ *P) {
+    PJ_XY xy;
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(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 col_urban_inverse (PJ_XY xy, PJ *P) {
+    PJ_LP lp;
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(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(col_urban) {
+    struct pj_opaque *Q = static_cast<struct pj_opaque*>(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 = col_urban_forward;
+    P->inv = col_urban_inverse;
+
+    return P;
+}
diff --git a/test/cli/testvarious b/test/cli/testvarious
index 6a77729e..292ee316 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} <<EOF
+4.8 -74.25
+EOF
+
 # Done!
 # do 'diff' with distribution results
 echo "diff ${OUT} with ${OUT}.dist"
diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist
index 22a26380..70b2ab6e 100644
--- a/test/cli/tv_out.dist
+++ b/test/cli/tv_out.dist
@@ -482,3 +482,6 @@ Check +proj=longlat +over +datum=WGS84 +to proj=merc +a=6378137 +b=6378137 +lat_
 ##############################################################
 Test EPSG:xxxx EPSG:yyyy filename
 2 49	2dN	49dE 0.000
+##############################################################
+Test Colombia Urban
+4.8 -74.25	122543.174	80859.033 0.000
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index cbba686c..cfce5041 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -6587,4 +6587,18 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+
+===============================================================================
+# Colombia Urbian
+# Test point from IOGP Publication 373-7-2 - Geomatics Guidance Note number 7, part 2 - March 2020
+===============================================================================
+
+-------------------------------------------------------------------------------
+operation +proj=col_urban +lat_0=4.68048611111111 +lon_0=-74.1465916666667 +x_0=92334.879 +y_0=109320.965 +h_0=2550 +ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 1 mm
+accept    -74.25    4.8
+expect    80859.033 122543.174
+roundtrip 1
+
 </gie-strict>