Central conic projection (gnomonic) implementation (as 'proj=ccon') (#662)

Central conic projection implemented as 'ccon'.

11 files changed, 273 insertions, 2 deletions

diff --git a/docs/plot/plotdefs.json b/docs/plot/plotdefs.json
index 7df9d390..b61466f0 100644
--- a/docs/plot/plotdefs.json
+++ b/docs/plot/plotdefs.json
@@ -143,6 +143,17 @@
         "type": "poly"
     },
     {
+        "filename": "ccon.png",
+        "latmax": 70,
+        "latmin": 34,
+        "lonmax": -30,
+        "lonmin": 68,
+        "name": "ccon",
+        "projstring": "+proj=ccon +lat_1=52 +lon_0=19",
+        "res": "med",
+        "type": "poly"
+    },
+    {
         "filename": "cea.png",
         "latmax": 90,
         "latmin": -90,
diff --git a/docs/source/projections/ccon.rst b/docs/source/projections/ccon.rst
new file mode 100644
index 00000000..6197f061
--- /dev/null
+++ b/docs/source/projections/ccon.rst
@@ -0,0 +1,118 @@
+.. _ccon:
+
+********************************************************************************
+Central Conic
+********************************************************************************
+
+.. image:: ./images/ccon.png
+   :scale: 50%
+   :alt:   Central Conic 
+
+This is central (centrographic) projection on cone tangent at ``lat_0`` latitude,
+identical with ``conic()`` projection from ``mapproj`` R package.
+
+Usage
+########
+
+This simple projection is rarely used, as it is not equidistant, equal-area, nor
+conformal. 
+
+An example of usage (and the main reason to implement this projection in proj4)
+is the ATPOL geobotanical grid of Poland, developed in Institute of Botany, 
+Jagiellonian University, Krakow, Poland in 1970s [Zajac1978]_. The grid was 
+originally handwritten on paper maps and further copied by hand. The projection
+(together with strange Earth radius) was chosen by its creators as the compromise
+fit to existing maps during first software development in DOS era. Many years later
+it is still de facto standard grid in Polish geobotanical research.
+
+The ATPOL coordinates can be achieved with with the following parameters:
+
+::
+
+   +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000 +x_0=330000 +y_0=-350000
+
+For more information see [Komsta2016]_ and [Verey2017]_.
+
+
+Forward projection
+==================
+
+.. math::
+
+   r = \cot \phi_0 - \tan (\phi - \phi_0)
+
+.. math::
+
+   x = r \sin (\lambda\sin\phi_0)
+
+.. math::
+
+   y = \cot \phi_0 - r \cos (\lambda\sin\phi_0) 
+
+
+Inverse projection
+==================
+
+.. math::
+
+   y = \cot \phi_0 - y
+
+.. math::
+
+   \phi = \phi_0 - \tan^{-1} ( \sqrt{x^2+y^2} - \cot \phi_0 )
+
+.. math::
+
+   \lambda = \frac{\tan^{-1} \sqrt{x^2+y^2}}{\sin \phi_0}
+
+Reference values
+==================
+
+For ATPOL to WGS84 test, run the following script:
+
+::
+
+   #!/bin/bash
+   cat << EOF | src/cs2cs -v -f "%E" +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000 +x_0=330000 +y_0=-350000 +to +proj=longlat +datum=WGS84 +no_defs
+   0 0
+   0 700000
+   700000 0
+   700000 700000
+   330000 350000
+   EOF
+
+It should result with
+
+::
+
+   1.384023E+01	5.503040E+01 0.000000E+00
+   1.451445E+01	4.877385E+01 0.000000E+00
+   2.478271E+01	5.500352E+01 0.000000E+00
+   2.402761E+01	4.875048E+01 0.000000E+00
+   1.900000E+01	5.200000E+01 0.000000E+00
+
+Analogous script can be run for reverse test:
+
+::
+
+   cat << EOF  | src/cs2cs -v -f "%E" +proj=longlat +datum=WGS84 +no_defs +to +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000 +x_0=330000 +y_0=-350000
+   24 55
+   15 49
+   24 49
+   19 52
+   EOF
+
+and it should give the following results:
+
+::
+
+   6.500315E+05	4.106162E+03 0.000000E+00
+   3.707419E+04	6.768262E+05 0.000000E+00
+   6.960534E+05	6.722946E+05 0.000000E+00
+   3.300000E+05	3.500000E+05 0.000000E+00
+
+
+
+
+
+
diff --git a/docs/source/projections/images/ccon.png b/docs/source/projections/images/ccon.png
new file mode 100644
index 00000000..ccc0ec9e
--- /dev/null
+++ b/docs/source/projections/images/ccon.png
diff --git a/docs/source/projections/index.rst b/docs/source/projections/index.rst
index ae8ef1ee..5ccf0045 100644
--- a/docs/source/projections/index.rst
+++ b/docs/source/projections/index.rst
@@ -20,6 +20,7 @@ Projections
    calcofi
    cass
    cc
+   ccon
    cea
    chamb
    collg
diff --git a/docs/source/references.rst b/docs/source/references.rst
index 26c4da08..31100ded 100644
--- a/docs/source/references.rst
+++ b/docs/source/references.rst
@@ -26,3 +26,9 @@ References
 .. [ONeilLaubscher1976] E.M. O'Neill and R.E. Laubscher, 1976, "Extended Studies of a Quadrilateralized Spherical Cube Earth Data Base". Tech. Rep. NEPRF 3-76 (CSC), Naval Environmental Prediction Research Facility.
 
 .. [LambersKolb2012] M. Lambers and A. Kolb, 2012, "Ellipsoidal Cube Maps for Accurate Rendering of Planetary-Scale Terrain Data", Proc. Pacfic Graphics (Short Papers).
+
+.. [Zajac1978] A. Zajac, 1978, "Atlas of distribution of vascular plants in Poland (ATPOL)". Taxon 27(5/6), 481–484. 
+
+.. [Komsta2016] L. Komsta, 2016, `ATPOL geobotanical grid revisited – a proposal of coordinate conversion algorithms <http://wydawnictwo.up.lublin.pl/annales/Agricultura/2016/1/03.pdf>`__. Annales UMCS Sectio E Agricultura 71(1), 31-37.
+
+.. [Verey2017] M. Verey, 2017, `Theoretical analysis and practical consequences of adopting an ATPOL grid model as a conical projection, defining the conversion of plane coordinates to the WGS - 84 ellipsoid <http://www.botany.pl/atpol/Siatka%20ATPOL%20w%20analitycznym%20ujeciu.pdf>`__. Fragmenta Floristica et Geobotanica Polonica (preprint submitted).
diff --git a/src/Makefile.am b/src/Makefile.am
index 3fdc26d7..096cc672 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,7 +46,7 @@ libproj_la_SOURCES = \
 	pj_list.h proj_internal.h\
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
-	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \
+	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c PJ_ccon.c\
 	PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
 	PJ_rpoly.c PJ_sconics.c proj_rouss.c \
 	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c PJ_gall.c PJ_geoc.c \
diff --git a/src/PJ_ccon.c b/src/PJ_ccon.c
new file mode 100644
index 00000000..cf8a9f68
--- /dev/null
+++ b/src/PJ_ccon.c
@@ -0,0 +1,102 @@
+/******************************************************************************
+ * Copyright (c) 2017, Lukasz Komsta
+ *
+ * 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 <errno.h>
+#include <proj.h>
+#include "projects.h"
+
+struct pj_opaque {
+    double phi1;
+    double ctgphi1;
+    double sinphi1;
+    double cosphi1;
+    double *en;	
+};
+
+PROJ_HEAD(ccon, "Central Conic")
+    "\n\tCentral Conic, Sph.\n\tlat_1=";
+
+
+
+static XY forward (LP lp, PJ *P) {          
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double r;
+
+    r = Q->ctgphi1 - tan(lp.phi - Q->phi1);
+    xy.x = r * sin(lp.lam * Q->sinphi1);
+    xy.y = Q->ctgphi1 - r * cos(lp.lam * Q->sinphi1);	
+
+    return xy;
+}
+
+
+static LP inverse (XY xy, PJ *P) {         
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    xy.y = Q->ctgphi1 - xy.y;
+    lp.phi = Q->phi1 - atan(hypot(xy.x,xy.y) - Q->ctgphi1);
+    lp.lam = atan2(xy.x,xy.y)/Q->sinphi1;
+
+    return lp;
+}
+
+
+static void *destructor (PJ *P, int errlev) {                       
+    if (0==P)
+        return 0;
+
+    if (0==P->opaque)
+        return pj_default_destructor (P, errlev);
+
+    pj_dealloc (P->opaque->en);
+    return pj_default_destructor (P, errlev);
+}
+
+
+PJ *PROJECTION(ccon) {
+
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+    P->destructor = destructor;
+
+    Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+
+    if (!(Q->en = pj_enfn(P->es)))
+        return pj_default_destructor(P, ENOMEM);
+
+    Q->sinphi1 = sin(Q->phi1);
+    Q->cosphi1 = cos(Q->phi1);
+    Q->ctgphi1 = Q->cosphi1/Q->sinphi1;
+
+
+    P->inv = inverse;
+    P->fwd = forward;
+    
+    return P;
+}
+
+
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 52abc665..88d88a97 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -54,6 +54,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_cart.c
         PJ_cass.c
         PJ_cc.c
+        PJ_ccon.c
         PJ_cea.c
         PJ_chamb.c
         PJ_collg.c
diff --git a/src/makefile.vc b/src/makefile.vc
index f2acf982..ef460719 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -11,7 +11,7 @@ azimuthal = \
 conic = \
 	PJ_aea.obj PJ_bipc.obj PJ_bonne.obj PJ_eqdc.obj \
 	PJ_imw_p.obj PJ_lcc.obj PJ_poly.obj \
-	PJ_rpoly.obj PJ_sconics.obj PJ_lcca.obj
+	PJ_rpoly.obj PJ_sconics.obj PJ_lcca.obj PJ_ccon.obj
 
 cylinder = \
 	PJ_cass.obj PJ_cc.obj PJ_cea.obj PJ_eqc.obj \
diff --git a/src/pj_list.h b/src/pj_list.h
index 093355bd..09a66034 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -21,6 +21,7 @@ PROJ_HEAD(calcofi, "Cal Coop Ocean Fish Invest Lines/Stations")
 PROJ_HEAD(cart,    "Geodetic/cartesian conversions")
 PROJ_HEAD(cass, "Cassini")
 PROJ_HEAD(cc, "Central Cylindrical")
+PROJ_HEAD(ccon, "Central Conic")
 PROJ_HEAD(cea, "Equal Area Cylindrical")
 PROJ_HEAD(chamb, "Chamberlin Trimetric")
 PROJ_HEAD(collg, "Collignon")
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index c6a801f8..6b820c01 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -502,6 +502,37 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+===============================================================================
+Central Conic
+	Sph
+	lat_1 
+===============================================================================
+
+-------------------------------------------------------------------------------
+operation +proj=pipeline +step +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +a=6390000 +x_0=330000 +y_0=-350000 +step +proj=axisswap +order=1,-2,3,4
+-------------------------------------------------------------------------------
+tolerance 0.00010 mm
+accept 24 55
+expect 650031.54109413219363 4106.1617770643609028
+accept 15 49
+expect 37074.189007307473069 676826.23559270039774
+accept 24 49
+expect 696053.36061617843913 672294.56795827199940
+accept 19 52
+expect 330000.00000000000000 350000.00000000000000
+
+direction inverse
+accept 0 0
+expect 13.840227318521004431 55.030403993648806391
+accept 0 700000
+expect 14.514453594615022781 48.773847834747808675
+accept 700000 0
+expect 24.782707184271129766 55.003515505218481835
+accept 700000 700000
+expect 24.027610763560529927 48.750476070495021286
+accept 330000 350000
+expect 19.000000000000000000 52.000000000000000000
+
 
 ===============================================================================
 Central Cylindrical