diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 1 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 1 | ||||
| -rw-r--r-- | src/pj_list.h | 1 | ||||
| -rw-r--r-- | src/proj_constants.h | 6 | ||||
| -rw-r--r-- | src/projections/col_urban.cpp | 76 |
5 files changed, 85 insertions, 0 deletions
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; +} |