NEWS

@@ -162,7 +162,7 @@ o The PROJ.4 project also distributes the datum-grid package,
 
  o Updated EPSG database to version 9.2.0.
 
- o Geodesic library updated to version 1.49.1-c.
+ o Geodesic library updated to version 1.49.2-c.
 
  o Support for analytical partial derivatives has been removed.
 
@@ -232,9 +232,9 @@ o The PROJ.4 project also distributes the datum-grid package,
  o Equations for meridian convergence and partial derivatives have
    been corrected for non-conformal projections. #526.
 
- o Fixed scaling of cartesian coordiantes in pj_transform(). #726.
+ o Fixed scaling of cartesian coordinates in pj_transform(). #726.
 
- o Additional bug fixes courtesy of Googles OSS-Fuzz program:
+ o Additional bug fixes courtesy of Google's OSS-Fuzz program:
    https://bugs.chromium.org/p/oss-fuzz/issues/list?can=1&q=proj4
 
 

debian/changelog

+proj (5.0.0~rc5-1) UNRELEASED; urgency=medium
+
+  * New upstream release candidate.
+  * Drop install-nad-gtx.patch, applied upstream.
+  * Update copyright years for Charles Karney.
+
+ -- Bas Couwenberg <sebastic@debian.org>  Fri, 23 Feb 2018 20:40:32 +0100
+
 proj (5.0.0~rc4-1~exp2) experimental; urgency=medium
 
   * Update datumgrids.shar (XZ compressed) with proj-datumgrid-1.7RC1.

debian/copyright

@@ -34,7 +34,7 @@ Comment: Version 4.3.3 was authored by:
    I did the work.  Essentially all work was done by Gerald Evenden.
 Copyright:             2016-2018, Thomas Knudsen
                        2016-2018, SDFE
-                       2012-2017, Charles Karney <charles@karney.com>
+                       2012-2018, Charles Karney <charles@karney.com>
                        2016-2017, Kristian Evers
                             2017, Lukasz Komsta
                             2016, Karsten Engsager

debian/patches/install-nad-gtx.patch

-Description: Also install .gtx files from proj-datumgrid.
-Author: Bas Couwenberg <sebastic@debian.org>
-Forwarded: https://github.com/OSGeo/proj.4/pull/813
-Applied-Upstream: https://github.com/OSGeo/proj.4/commit/b65e8ef7b29f19c6b4b4a49cf95766ec1f752491
-
---- a/nad/Makefile.am
-+++ b/nad/Makefile.am
-@@ -50,7 +50,7 @@ install-data-local: process-nad2bin
- 	else \
- 	  echo "nad2nad NADCON source files not present"; \
- 	fi
--	@for gridfile in $(NADPATH)/*.gsb $(NADPATH)/ntv1_can.dat dummy \
-+	@for gridfile in $(NADPATH)/*.gsb $(NADPATH)/*.gtx $(NADPATH)/ntv1_can.dat dummy \
- 	                 $(NADPATH)/alaska $(NADPATH)/conus $(NADPATH)/hawaii $(NADPATH)/null \
- 	                 $(NADPATH)/prvi $(NADPATH)/stgeorge $(NADPATH)/stlrnc $(NADPATH)/stpaul \
- 	                 $(NADPATH)/FL $(NADPATH)/MD $(NADPATH)/TN $(NADPATH)/WI $(NADPATH)/WO; do \
-@@ -63,7 +63,7 @@ install-data-local: process-nad2bin
- check-local: process-nad2bin
- # Check if we are doing an out-of-tree build
- 	@if test ! -f epsg; then \
--	  for gridfile in $(NADPATH)/*.gsb $(NADPATH)/ntv1_can.dat dummy ; do \
-+	  for gridfile in $(NADPATH)/*.gsb $(NADPATH)/*.gtx $(NADPATH)/ntv1_can.dat dummy ; do \
- 	    if test "$$gridfile" != "dummy" -a -f "$$gridfile" ; then \
- 	      cp $$gridfile .; \
- 	    fi; \
---- a/nad/CMakeLists.txt
-+++ b/nad/CMakeLists.txt
-@@ -23,7 +23,8 @@ set(PROJ_DICTIONARY  epsg
- #
- 
- file(GLOB GSB_FILES  *.gsb)
--set(GRIDSHIFT_FILES ${GSB_FILES})
-+file(GLOB GTX_FILES  *.gtx)
-+set(GRIDSHIFT_FILES ${GSB_FILES} ${GTX_FILES})
- set(GRIDSHIFT_FILES ${GRIDSHIFT_FILES}
-                      )
- option(CONVERT_DATA "convert some ascii file to binary file for use in proj4" OFF)

debian/patches/series

-install-nad-gtx.patch

man/man1/cct.1

@@ -30,7 +30,7 @@ professor of Geodesy at the University of Copenhagen, mentor and advisor
 for a generation of Danish geodesists, colleague and collaborator for
 two generations of global geodesists, Secretary General for the
 International Association of Geodesy, IAG (1995--2007), fellow of the
-Amercan Geophysical Union (1991), recipient of the IAG Levallois Medal
+American Geophysical Union (1991), recipient of the IAG Levallois Medal
 (2007), the European Geosciences Union Vening Meinesz Medal (2008), and
 of numerous other honours.
 .PP
@@ -40,7 +40,7 @@ the development of geodesy - both through his scientific contributions,
 comprising more than 250 publications, and by his mentoring and teaching
 of the next generations of geodesists.
 .PP
-As Christian was an avid Fortran programmer, and a keen Unix connoiseur,
+As Christian was an avid Fortran programmer, and a keen Unix connoisseur,
 he would have enjoyed to know that his initials would be used to name a
 modest Unix style transformation filter, hinting at the tireless aspect
 of his personality, which was certainly one of the reasons he accomplished

man/man1/proj.1

@@ -119,7 +119,7 @@ output from x-y or longitude-latitude to y-x or latitude-longitude.
 .TP
 .BI \-S
 Causes estimation of
-.I meridinal
+.I meridional
 and
 .I parallel
 scale factors,
@@ -131,7 +131,7 @@ and
 and
 .I minimum
 scale factors to be listed between <> for each input point.
-For conformal projections meridinal and parallel scales factors
+For conformal projections meridional and parallel scales factors
 will be equal and angular distortion zero.
 Equal area projections will have an area factor of 1.
 .TP

nad/CMakeLists.txt

@@ -23,7 +23,8 @@ set(PROJ_DICTIONARY  epsg
 #
 
 file(GLOB GSB_FILES  *.gsb)
-set(GRIDSHIFT_FILES ${GSB_FILES})
+file(GLOB GTX_FILES  *.gtx)
+set(GRIDSHIFT_FILES ${GSB_FILES} ${GTX_FILES})
 set(GRIDSHIFT_FILES ${GRIDSHIFT_FILES}
                      )
 option(CONVERT_DATA "convert some ascii file to binary file for use in proj4" OFF)

nad/Makefile.am

@@ -50,7 +50,7 @@ install-data-local: process-nad2bin
 	else \
 	  echo "nad2nad NADCON source files not present"; \
 	fi
-	@for gridfile in $(NADPATH)/*.gsb $(NADPATH)/ntv1_can.dat dummy \
+	@for gridfile in $(NADPATH)/*.gsb $(NADPATH)/*.gtx $(NADPATH)/ntv1_can.dat dummy \
 	                 $(NADPATH)/alaska $(NADPATH)/conus $(NADPATH)/hawaii $(NADPATH)/null \
 	                 $(NADPATH)/prvi $(NADPATH)/stgeorge $(NADPATH)/stlrnc $(NADPATH)/stpaul \
 	                 $(NADPATH)/FL $(NADPATH)/MD $(NADPATH)/TN $(NADPATH)/WI $(NADPATH)/WO; do \
@@ -63,7 +63,7 @@ install-data-local: process-nad2bin
 check-local: process-nad2bin
 # Check if we are doing an out-of-tree build
 	@if test ! -f epsg; then \
-	  for gridfile in $(NADPATH)/*.gsb $(NADPATH)/ntv1_can.dat dummy ; do \
+	  for gridfile in $(NADPATH)/*.gsb $(NADPATH)/*.gtx $(NADPATH)/ntv1_can.dat dummy ; do \
 	    if test "$$gridfile" != "dummy" -a -f "$$gridfile" ; then \
 	      cp $$gridfile .; \
 	    fi; \

nad/Makefile.in

@@ -543,7 +543,7 @@ install-data-local: process-nad2bin
 	else \
 	  echo "nad2nad NADCON source files not present"; \
 	fi
-	@for gridfile in $(NADPATH)/*.gsb $(NADPATH)/ntv1_can.dat dummy \
+	@for gridfile in $(NADPATH)/*.gsb $(NADPATH)/*.gtx $(NADPATH)/ntv1_can.dat dummy \
 	                 $(NADPATH)/alaska $(NADPATH)/conus $(NADPATH)/hawaii $(NADPATH)/null \
 	                 $(NADPATH)/prvi $(NADPATH)/stgeorge $(NADPATH)/stlrnc $(NADPATH)/stpaul \
 	                 $(NADPATH)/FL $(NADPATH)/MD $(NADPATH)/TN $(NADPATH)/WI $(NADPATH)/WO; do \
@@ -556,7 +556,7 @@ install-data-local: process-nad2bin
 check-local: process-nad2bin
 # Check if we are doing an out-of-tree build
 	@if test ! -f epsg; then \
-	  for gridfile in $(NADPATH)/*.gsb $(NADPATH)/ntv1_can.dat dummy ; do \
+	  for gridfile in $(NADPATH)/*.gsb $(NADPATH)/*.gtx $(NADPATH)/ntv1_can.dat dummy ; do \
 	    if test "$$gridfile" != "dummy" -a -f "$$gridfile" ; then \
 	      cp $$gridfile .; \
 	    fi; \

nad/README

@@ -67,10 +67,3 @@ TN.lla.Z    --- Tennessee
 WI.lla.Z    --- Wisconsin
 WO.lla.Z    --- Washington, Oregon, N. California
 
-	Other grid shift files
-ntv1_can.dat	--- Canadian NTv1 grid shift file (NAD27-->NAD83)
-ntf_r93.gsb	--- French NTv2 grid shift file (NTF-->RGF93)
-BETA2007.gsb	--- German NTv2 grid shift file (DE_DHDN-->ETRS89) :
-					http://crs.bkg.bund.de/crseu/crs/descrtrans/BeTA/de_dhdn2etrs_beta.php
-					Confirmed with Uwe Schmitz <uwe.schmitz@bezreg-koeln.nrw.de> that free
-					redistribution is allowed and welcome.

nad/test27

@@ -4,7 +4,7 @@
 #   the second pair of numbers are respective easting and northing output.
 #
 # Proj will vary in the .001ft range with projections using Transverse
-#   Mercator due to greater precision of meridinal distance function.
+#   Mercator due to greater precision of meridional distance function.
 #
 NAD_DIR=`dirname $0`
 EXE=$1

nad/test83

@@ -5,7 +5,7 @@
 #   the second pair of numbers are respective easting and northing output.
 #
 # Proj will vary in the .001ft range with projections using Transverse
-#   Mercator due to greater precision of meridinal distance function.
+#   Mercator due to greater precision of meridional distance function.
 #
 NAD_DIR=`dirname $0`
 EXE=$1

src/PJ_calcofi.c

@@ -13,7 +13,7 @@ PROJ_HEAD(calcofi,
 
 /* Conversions for the California Cooperative Oceanic Fisheries Investigations
 Line/Station coordinate system following the algorithm of:
-Eber, L.E., and  R.P. Hewitt. 1979. Conversion algorithms for the CALCOFI
+Eber, L.E., and  R.P. Hewitt. 1979. Conversion algorithms for the CalCOFI
 station grid. California Cooperative Oceanic Fisheries Investigations Reports
 20:135-137. (corrected for typographical errors).
 http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf

src/PJ_horner.c

@@ -96,6 +96,8 @@ static HORNER *horner_alloc (size_t order, int complex_polynomia);
 static void    horner_free (HORNER *h);
 
 struct horner {
+    int    uneg;     /* u axis negated? */
+    int    vneg;     /* v axis negated? */
     int    order;    /* maximum degree of polynomium */
     int    coefs;    /* number of coefficients for each polynomium  */
     double range;    /* radius of the region of validity */
@@ -340,11 +342,19 @@ polynomial evaluation engine.
         c  =  cb + sz;
         e  =  position.u - transformation->fwd_origin->u;
         n  =  position.v - transformation->fwd_origin->v;
+        if (transformation->uneg)
+            e  =  -e;
+        if (transformation->vneg)
+            n  =  -n;
     } else {                                              /* inverse */
         cb =  transformation->inv_c;
         c  =  cb + sz;
         e  =  position.u - transformation->inv_origin->u;
         n  =  position.v - transformation->inv_origin->v;
+        if (transformation->uneg)
+            e  =  -e;
+        if (transformation->vneg)
+            n  =  -n;
     }
 
     if ((fabs(n) > range) || (fabs(e) > range)) {
@@ -454,6 +464,10 @@ PJ *PROJECTION(horner) {
     P->opaque = (void *) Q;
 
     if (complex_horner) {
+        /* Westings and/or southings? */
+        Q->uneg = pj_param_exists (P->params, "uneg") ? 1 : 0;
+        Q->vneg = pj_param_exists (P->params, "vneg") ? 1 : 0;
+
         n = 2*degree + 2;
         if (0==parse_coefs (P, Q->fwd_c, "fwd_c", n))
             return horner_freeup (P, PJD_ERR_MISSING_ARGS);
@@ -461,8 +475,8 @@ PJ *PROJECTION(horner) {
             return horner_freeup (P, PJD_ERR_MISSING_ARGS);
         P->fwd4d = complex_horner_forward_4d;
         P->inv4d = complex_horner_reverse_4d;
-
     }
+
     else {
         n = horner_number_of_coefficients (degree);
         if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n))

src/cct.c

@@ -24,7 +24,7 @@ professor of Geodesy at the University of Copenhagen, mentor and advisor
 for a generation of Danish geodesists, colleague and collaborator for
 two generations of global geodesists, Secretary General for the
 International Association of Geodesy, IAG (1995--2007), fellow of the
-Amercan Geophysical Union (1991), recipient of the IAG Levallois Medal
+American Geophysical Union (1991), recipient of the IAG Levallois Medal
 (2007), the European Geosciences Union Vening Meinesz Medal (2008), and
 of numerous other honours.
 
@@ -34,7 +34,7 @@ the development of geodesy - both through his scientific contributions,
 comprising more than 250 publications, and by his mentoring and teaching
 of the next generations of geodesists.
 
-As Christian was an avid Fortran programmer, and a keen Unix connoiseur,
+As Christian was an avid Fortran programmer, and a keen Unix connoisseur,
 he would have enjoyed to know that his initials would be used to name a
 modest Unix style transformation filter, hinting at the tireless aspect
 of his personality, which was certainly one of the reasons he accomplished

src/geodesic.c

@@ -18,7 +18,7 @@
  *
  * See the comments in geodesic.h for documentation.
  *
- * Copyright (c) Charles Karney (2012-2017) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2018) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  */
@@ -89,10 +89,14 @@ static void Init() {
     tolb = tol0 * tol2;
     xthresh = 1000 * tol2;
     degree = pi/180;
+    #if defined(NAN)
+    NaN = NAN;
+    #else
     {
       real minus1 = -1;
       NaN = sqrt(minus1);
     }
+    #endif
     init = 1;
   }
 }
@@ -1809,13 +1813,13 @@ int transitdirect(real lon1, real lon2) {
 #if HAVE_C99_MATH
   lon1 = remainder(lon1, (real)(720));
   lon2 = remainder(lon2, (real)(720));
-  return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) -
-           (lon1 >= 0 && lon1 < 360 ? 0 : 1) );
+  return ( (lon2 <= 0 && lon2 > -360 ? 1 : 0) -
+           (lon1 <= 0 && lon1 > -360 ? 1 : 0) );
 #else
   lon1 = fmod(lon1, (real)(720));
   lon2 = fmod(lon2, (real)(720));
-  return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) -
-           ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) );
+  return ( ((lon2 <= 0 && lon2 > -360) || lon2 > 360 ? 1 : 0) -
+           ((lon1 <= 0 && lon1 > -360) || lon1 > 360 ? 1 : 0) );
 #endif
 }
 

src/geodesic.h

@@ -107,7 +107,7 @@
  * twice about restructuring the internals of the C code since this may make
  * porting fixes from the C++ code more difficult.
  *
- * Copyright (c) Charles Karney (2012-2017) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2012-2018) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  *
@@ -132,7 +132,7 @@
  * The patch level of the geodesic library.  (This tracks the version of
  * GeographicLib.)
  **********************************************************************/
-#define GEODESIC_VERSION_PATCH 1
+#define GEODESIC_VERSION_PATCH 2
 
 /**
  * Pack the version components into a single integer.  Users should not rely on

src/geodtest.c

@@ -4,7 +4,7 @@
  *
  * Run these tests by configuring with cmake and running "make test".
  *
- * Copyright (c) Charles Karney (2015-2017) <charles@karney.com> and licensed
+ * Copyright (c) Charles Karney (2015-2018) <charles@karney.com> and licensed
  * under the MIT/X11 License.  For more information, see
  * https://geographiclib.sourceforge.io/
  **********************************************************************/
@@ -688,6 +688,48 @@ static int GeodSolve78() {
   return result;
 }
 
+static int GeodSolve80() {
+  /* Some tests to add code coverage: computing scale in special cases + zero
+     length geodesic (includes GeodSolve80 - GeodSolve83) + using an incapable
+     line. */
+  double a12, s12, azi1, azi2, m12, M12, M21, S12;
+  struct geod_geodesic g;
+  struct geod_geodesicline l;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_geninverse(&g, 0, 0, 0, 90, 0, 0, 0, 0, &M12, &M21, 0);
+  result += assertEquals(M12, -0.0052842753, 0.5e-10);
+  result += assertEquals(M21, -0.0052842753, 0.5e-10);
+  geod_geninverse(&g, 0, 0, 1e-6, 1e-6, 0, 0, 0, 0, &M12, &M21, 0);
+  result += assertEquals(M12, 1, 0.5e-10);
+  result += assertEquals(M21, 1, 0.5e-10);
+  a12 = geod_geninverse(&g, 20.001, 0, 20.001, 0,
+                        &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
+  result += a12 == 0 ? 0 : 1;
+  result += s12 == 0 ? 0 : 1;
+  result += azi1 == 180 ? 0 : 1;
+  result += azi2 == 180 ? 0 : 1;
+  result += m12 == 0 ? 0 : 1;
+  result += assertEquals(M12, 1, 1e-15);
+  result += assertEquals(M21, 1, 1e-15);
+  result += S12 == 0 ? 0 : 1;
+  a12 = geod_geninverse(&g, 90, 0, 90, 180,
+                        &s12, &azi1, &azi2, &m12, &M12, &M21, &S12);
+  result += a12 == 0 ? 0 : 1;
+  result += s12 == 0 ? 0 : 1;
+  result += azi1 == 0 ? 0 : 1;
+  result += azi2 == 180 ? 0 : 1;
+  result += m12 == 0 ? 0 : 1;
+  result += assertEquals(M12, 1, 1e-15);
+  result += assertEquals(M21, 1, 1e-15);
+  result += assertEquals(S12, 127516405431022, 0.5);
+  /* An incapable line which can't take distance as input */
+  geod_lineinit(&l, &g, 1, 2, 90, GEOD_LATITUDE);
+  a12 = geod_genposition(&l, 0, 1000, 0, 0, 0, 0, 0, 0, 0, 0);
+  result += a12 != a12 ? 0 : 1;
+  return result;
+}
+
 static int Planimeter0() {
   /* Check fix for pole-encircling bug found 2011-03-16 */
   double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}};
@@ -787,6 +829,188 @@ static int Planimeter13() {
   return result;
 }
 
+static int Planimeter15() {
+  /* Coverage tests, includes Planimeter15 - Planimeter18 (combinations of
+     reverse and sign) + calls to testpoint, testedge, geod_polygonarea. */
+  struct geod_geodesic g;
+  struct geod_polygon p;
+  double lat[] = {2, 1, 3}, lon[] = {1, 2, 3};
+  double area, s12, azi1;
+  double r = 18454562325.45119,
+    a0 = 510065621724088.5093;  /* ellipsoid area */
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_polygon_init(&p, 0);
+  geod_polygon_addpoint(&g, &p, lat[0], lon[0]);
+  geod_polygon_addpoint(&g, &p, lat[1], lon[1]);
+  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 0, 1, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 0, 0, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 1, 1, &area, 0);
+  result += assertEquals(area, -r, 0.5);
+  geod_polygon_testpoint(&g, &p, lat[2], lon[2], 1, 0, &area, 0);
+  result += assertEquals(area, a0-r, 0.5);
+  geod_inverse(&g, lat[1], lon[1], lat[2], lon[2], &s12, &azi1, 0);
+  geod_polygon_testedge(&g, &p, azi1, s12, 0, 1, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  geod_polygon_testedge(&g, &p, azi1, s12, 0, 0, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  geod_polygon_testedge(&g, &p, azi1, s12, 1, 1, &area, 0);
+  result += assertEquals(area, -r, 0.5);
+  geod_polygon_testedge(&g, &p, azi1, s12, 1, 0, &area, 0);
+  result += assertEquals(area, a0-r, 0.5);
+  geod_polygon_addpoint(&g, &p, lat[2], lon[2]);
+  geod_polygon_compute(&g, &p, 0, 1, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  geod_polygon_compute(&g, &p, 0, 0, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  geod_polygon_compute(&g, &p, 1, 1, &area, 0);
+  result += assertEquals(area, -r, 0.5);
+  geod_polygon_compute(&g, &p, 1, 0, &area, 0);
+  result += assertEquals(area, a0-r, 0.5);
+  geod_polygonarea(&g, lat, lon, 3, &area, 0);
+  result += assertEquals(area, r, 0.5);
+  return result;
+}
+
+static int Planimeter19() {
+  /* Coverage tests, includes Planimeter19 - Planimeter20 (degenerate
+     polygons) + extra cases.  */
+  struct geod_geodesic g;
+  struct geod_polygon p;
+  double area, perim;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_polygon_init(&p, 0);
+  geod_polygon_compute(&g, &p, 0, 1, &area, &perim);
+  result += area == 0 ? 0 : 1;
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, &area, &perim);
+  result += area == 0 ? 0 : 1;
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, &area, &perim);
+  result += area != area ? 0 : 1;
+  result += perim != perim ? 0 : 1;
+  geod_polygon_addpoint(&g, &p, 1, 1);
+  geod_polygon_compute(&g, &p, 0, 1, &area, &perim);
+  result += area == 0 ? 0 : 1;
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_init(&p, 1);
+  geod_polygon_compute(&g, &p, 0, 1, 0, &perim);
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, 0, &perim);
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, 0, &perim);
+  result += perim != perim ? 0 : 1;
+  geod_polygon_addpoint(&g, &p, 1, 1);
+  geod_polygon_compute(&g, &p, 0, 1, 0, &perim);
+  result += perim == 0 ? 0 : 1;
+  return result;
+}
+
+static int Planimeter21() {
+  /* Some test to add code coverage: multiple circlings of pole (includes
+     Planimeter21 - Planimeter28) + invocations via testpoint and testedge.
+     Some of the results for i = 4 in the loop are wrong because we don't
+     reduce the area to the allowed range correctly.  However these cases are
+     not "simple" polygons, so we'll defer fixing the problem for now.
+   */
+  struct geod_geodesic g;
+  struct geod_polygon p;
+  double area, lat = 45,
+    a = 39.2144607176828184218, s = 8420705.40957178156285,
+    r = 39433884866571.4277,    /* Area for one circuit */
+    a0 = 510065621724088.5093;  /* Ellipsoid area */
+  int result = 0, i;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_polygon_init(&p, 0);
+  geod_polygon_addpoint(&g, &p, lat,  60);
+  geod_polygon_addpoint(&g, &p, lat, 180);
+  geod_polygon_addpoint(&g, &p, lat, -60);
+  geod_polygon_addpoint(&g, &p, lat,  60);
+  geod_polygon_addpoint(&g, &p, lat, 180);
+  geod_polygon_addpoint(&g, &p, lat, -60);
+  for (i = 3; i <= 4; ++i) {
+    geod_polygon_addpoint(&g, &p, lat,  60);
+    geod_polygon_addpoint(&g, &p, lat, 180);
+    geod_polygon_testpoint(&g, &p, lat, -60, 0, 1, &area, 0);
+    if (i != 4) result += assertEquals(area,  i*r, 0.5);
+    geod_polygon_testpoint(&g, &p, lat, -60, 0, 0, &area, 0);
+    if (i != 4) result += assertEquals(area,  i*r, 0.5);
+    geod_polygon_testpoint(&g, &p, lat, -60, 1, 1, &area, 0);
+    if (i != 4) result += assertEquals(area, -i*r, 0.5);
+    geod_polygon_testpoint(&g, &p, lat, -60, 1, 0, &area, 0);
+    result += assertEquals(area, -i*r + a0, 0.5);
+    geod_polygon_testedge(&g, &p, a, s, 0, 1, &area, 0);
+    if (i != 4) result += assertEquals(area,  i*r, 0.5);
+    geod_polygon_testedge(&g, &p, a, s, 0, 0, &area, 0);
+    if (i != 4) result += assertEquals(area,  i*r, 0.5);
+    geod_polygon_testedge(&g, &p, a, s, 1, 1, &area, 0);
+    if (i != 4) result += assertEquals(area, -i*r, 0.5);
+    geod_polygon_testedge(&g, &p, a, s, 1, 0, &area, 0);
+    result += assertEquals(area, -i*r + a0, 0.5);
+    geod_polygon_addpoint(&g, &p, lat, -60);
+    geod_polygon_compute(&g, &p, 0, 1, &area, 0);
+    if (i != 4) result += assertEquals(area,  i*r, 0.5);
+    geod_polygon_compute(&g, &p, 0, 0, &area, 0);
+    if (i != 4) result += assertEquals(area,  i*r, 0.5);
+    geod_polygon_compute(&g, &p, 1, 1, &area, 0);
+    if (i != 4) result += assertEquals(area, -i*r, 0.5);
+    geod_polygon_compute(&g, &p, 1, 0, &area, 0);
+    result += assertEquals(area, -i*r + a0, 0.5);
+  }
+  return result;
+}
+
+static int AddEdge1() {
+  /* Check fix to transitdirect vs transit zero handling inconsistency */
+  struct geod_geodesic g;
+  struct geod_polygon p;
+  double area;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_polygon_init(&p, 0);
+  geod_polygon_addpoint(&g, &p, 0, 0);
+  geod_polygon_addedge(&g, &p,  90, 1000);
+  geod_polygon_addedge(&g, &p,   0, 1000);
+  geod_polygon_addedge(&g, &p, -90, 1000);
+  geod_polygon_compute(&g, &p, 0, 1, &area, 0);
+  result += assertEquals(area, 1000000.0, 0.01);
+  return result;
+}
+
+static int EmptyPoly() {
+  struct geod_geodesic g;
+  struct geod_polygon p;
+  double perim, area;
+  int result = 0;
+  geod_init(&g, wgs84_a, wgs84_f);
+  geod_polygon_init(&p, 0);
+  geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, &area, &perim);
+  result += area == 0 ? 0 : 1;
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, &area, &perim);
+  result += area != area ? 0 : 1;
+  result += perim != perim ? 0 : 1;
+  geod_polygon_compute(&g, &p, 0, 1, &area, &perim);
+  result += area == 0 ? 0 : 1;
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_init(&p, 1);
+  geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, 0, &perim);
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, 0, &perim);
+  result += perim != perim ? 0 : 1;
+  geod_polygon_compute(&g, &p, 0, 1, 0, &perim);
+  result += perim == 0 ? 0 : 1;
+  geod_polygon_addpoint(&g, &p, 1, 1);
+  geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, 0, &perim);
+  result += assertEquals(perim, 1000, 1e-10);
+  geod_polygon_testpoint(&g, &p, 2, 2, 0, 1, 0, &perim);
+  result += assertEquals(perim, 156876.149, 0.5e-3);
+  return result;
+}
+
 int main() {
   int n = 0, i;
   if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);}
@@ -818,11 +1042,17 @@ int main() {
   if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);}
   if ((i = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);}
   if ((i = GeodSolve78())) {++n; printf("GeodSolve78 fail: %d\n", i);}
+  if ((i = GeodSolve80())) {++n; printf("GeodSolve80 fail: %d\n", i);}
   if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);}
   if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);}
   if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);}
   if ((i = Planimeter12())) {++n; printf("Planimeter12 fail: %d\n", i);}
   if ((i = Planimeter13())) {++n; printf("Planimeter13 fail: %d\n", i);}
+  if ((i = Planimeter15())) {++n; printf("Planimeter15 fail: %d\n", i);}
+  if ((i = Planimeter19())) {++n; printf("Planimeter19 fail: %d\n", i);}
+  if ((i = Planimeter21())) {++n; printf("Planimeter21 fail: %d\n", i);}
+  if ((i = AddEdge1())) {++n; printf("AddEdge1 fail: %d\n", i);}
+  if ((i = EmptyPoly())) {++n; printf("EmptyPoly fail: %d\n", i);}
   return n;
 }
 

src/gie.c

@@ -146,7 +146,7 @@ static ffio *ffio_create (const char **tags, size_t n_tags, size_t max_record_si
 
 static const char *gie_tags[] = {
     "<gie>", "operation", "accept", "expect", "roundtrip", "banner", "verbose",
-    "direction", "tolerance", "ignore", "builtins", "echo", "</gie>"
+    "direction", "tolerance", "ignore", "builtins", "echo", "skip", "</gie>"
 };
 
 static const size_t n_gie_tags = sizeof gie_tags / sizeof gie_tags[0];
@@ -177,6 +177,7 @@ typedef struct {
     PJ_COORD a, b, c, e;
     PJ_DIRECTION dir;
     int verbosity;
+    int skip;
     int op_id;
     int op_ok,    op_ko,    op_skip;
     int total_ok, total_ko, total_skip;
@@ -192,6 +193,8 @@ typedef struct {
 ffio *F = 0;
 
 static gie_ctx T;
+int tests=0, succs=0, succ_fails=0, fail_fails=0, succ_rtps=0, fail_rtps=0;
+int succ_builtins=0, fail_builtins=0;
 
 
 static const char delim[] = {"-------------------------------------------------------------------------------\n"};
@@ -299,6 +302,13 @@ int main (int argc, char **argv) {
         fprintf (T.fout, "%sGrand total: %d. Success: %d, Skipped: %d, Failure: %d\n",
                  delim, T.grand_ok+T.grand_ko+T.grand_skip, T.grand_ok, T.grand_skip, T.grand_ko);
         fprintf (T.fout, "%s", delim);
+        if (T.verbosity > 1) {
+            fprintf (T.fout, "Failing roundtrips: %4d,    Succeeding roundtrips: %4d\n", fail_rtps, succ_rtps);
+            fprintf (T.fout, "Failing failures:   %4d,    Succeeding failures:   %4d\n", fail_fails, succ_fails);
+            fprintf (T.fout, "Failing builtins:   %4d,    Succeeding builtins:   %4d\n", fail_builtins, succ_builtins);
+            fprintf (T.fout, "Internal counters:                            %4.4d(%4.4d)\n", tests, succs);
+            fprintf (T.fout, "%s", delim);
+        }
     }
     else
         if (T.grand_ko)
@@ -315,21 +325,54 @@ int main (int argc, char **argv) {
 static int another_failure (void) {
     T.op_ko++;
     T.total_ko++;
+    proj_errno_reset (T.P);
     return 0;
 }
 
 static int another_skip (void) {
     T.op_skip++;
     T.total_skip++;
+    proj_errno_reset (T.P);
     return 0;
 }
 
 static int another_success (void) {
     T.op_ok++;
     T.total_ok++;
+    proj_errno_reset (T.P);
     return 0;
 }
 
+static int another_succeeding_failure (void) {
+    succ_fails++;
+    return another_success ();
+}
+
+static int another_failing_failure (void) {
+    fail_fails++;
+    return another_failure ();
+}
+
+static int another_succeeding_roundtrip (void) {
+    succ_rtps++;
+    return another_success ();
+}
+
+static int another_failing_roundtrip (void) {
+    fail_rtps++;
+    return another_failure ();
+}
+
+static int another_succeeding_builtin (void) {
+    succ_builtins++;
+    return another_success ();
+}
+
+static int another_failing_builtin (void) {
+    fail_builtins++;
+    return another_failure ();
+}
+
 
 static int process_file (const char *fname) {
     FILE *f;
@@ -339,6 +382,9 @@ static int process_file (const char *fname) {
     T.op_ko = T.total_ko = 0;
     T.op_skip = T.total_skip = 0;
 
+    if (T.skip)
+        return proj_destroy (T.P), T.P = 0, 0;
+
     f = fopen (fname, "rt");
     if (0==f) {
         if (T.verbosity > 0) {
@@ -556,27 +602,27 @@ using the "builtins" command verb.
     i = pj_unitconvert_selftest ();
     if (i!=0) {
         fprintf (T.fout, "pj_unitconvert_selftest fails with %d\n", i);
-        another_failure();
+        another_failing_builtin();
     }
     else
-        another_success ();
+        another_succeeding_builtin ();
 
 
     i = pj_cart_selftest ();
     if (i!=0) {
         fprintf (T.fout, "pj_cart_selftest fails with %d\n", i);
-        another_failure();
+        another_failing_builtin();
     }
     else
-        another_success ();
+        another_succeeding_builtin ();
 
     i = pj_horner_selftest ();
     if (i!=0) {
         fprintf (T.fout, "pj_horner_selftest fails with %d\n", i);
-        another_failure();
+        another_failing_builtin();
     }
     else
-        another_success ();
+        another_succeeding_builtin ();
 
     return 0;
 }
@@ -674,7 +720,7 @@ back/forward transformation pairs.
 
     r = proj_roundtrip (T.P, T.dir, ntrips, &coo);
     if (r <= d)
-        return another_success ();
+        return another_succeeding_roundtrip ();
 
     if (T.verbosity > -1) {
         if (0==T.op_ko && T.verbosity < 2)
@@ -683,7 +729,7 @@ back/forward transformation pairs.
         fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) F->lineno);
         fprintf (T.fout, "     roundtrip deviation: %.6f mm, expected: %.6f mm\n", 1000*r, 1000*d);
     }
-    return another_failure ();
+    return another_failing_roundtrip ();
 }
 
 
@@ -723,7 +769,7 @@ static int expect_message_cannot_parse (const char *args) {
 }
 
 static int expect_failure_with_errno_message (int expected, int got) {
-    another_failure ();
+    another_failing_failure ();
 
     if (T.verbosity < 0)
         return 1;
@@ -779,14 +825,16 @@ Tell GIE what to expect, when transforming the ACCEPTed input
             if (expect_failure_with_errno && proj_errno (T.P)!=expect_failure_with_errno)
                 return expect_failure_with_errno_message (expect_failure_with_errno, proj_errno(T.P));
 
-            return another_success ();
+            return another_succeeding_failure ();
         }
 
         /* Otherwise, it's a true failure */
         banner (T.operation);
-        errmsg(3, "%sInvalid operation definition in line no. %d: %s\n",
-            delim, (int) T.operation_lineno, pj_strerrno(proj_errno(T.P)));
-        return another_failure ();
+        errmsg (3, "%sInvalid operation definition in line no. %d:\n       %s (errno=%s/%d)\n",
+            delim, (int) T.operation_lineno, pj_strerrno(proj_errno(T.P)),
+            err_const_from_errno (proj_errno(T.P)), proj_errno(T.P)
+        );
+        return another_failing_failure ();
     }
 
     /* We may still successfully fail even if the proj_create succeeded */
@@ -797,20 +845,24 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
         co = expect_trans_n_dim (ci);
 
-        /* Failed to fail? - that's a failure */
-        if (co.xyz.x!=HUGE_VAL)
-            return another_failure ();
-
         if (expect_failure_with_errno) {
-            printf ("errno=%d, expected=%d\n", proj_errno (T.P), expect_failure_with_errno);
             if (proj_errno (T.P)==expect_failure_with_errno)
-                return another_success ();
-
-            return another_failure ();
+                return another_succeeding_failure ();
+            printf ("errno=%d, expected=%d\n", proj_errno (T.P), expect_failure_with_errno);
+            return another_failing_failure ();
         }
 
-        /* Yes, we failed successfully */
-        return another_success ();
+
+        /* Succeeded in failing? - that's a success */
+        if (co.xyz.x==HUGE_VAL)
+            return another_succeeding_failure ();
+
+        /* Failed to fail? - that's a failure */
+        banner (T.operation);
+        errmsg (3, "%sFailed to fail. Operation definition in line no. %d\n",
+            delim, (int) T.operation_lineno
+        );
+        return another_failing_failure ();
     }
 
 
@@ -822,10 +874,12 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         printf ("left: %d   right:  %d\n", T.P->left, T.P->right);
     }
 
+    tests++;
     T.e  =  parse_coord (args);
     if (HUGE_VAL==T.e.v[0])
         return expect_message_cannot_parse (args);
 
+
     /* expected angular values, probably in degrees */
     ce = proj_angular_output (T.P, T.dir)? torad_coord (T.P, T.dir, T.e): T.e;
     if (T.verbosity > 3)
@@ -836,8 +890,14 @@ Tell GIE what to expect, when transforming the ACCEPTed input
     if (T.verbosity > 3)
         printf ("ACCEPTS  %.12f  %.12f  %.12f  %.12f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
 
-    /* angular output from proj_trans comes in radians */
+    /* do the transformation, but mask off dimensions not given in expect-ation */
     co = expect_trans_n_dim (ci);
+    if (T.dimensions_given < 4)
+        co.v[3] = 0;
+    if (T.dimensions_given < 3)
+        co.v[2] = 0;
+
+    /* angular output from proj_trans comes in radians */
     T.b = proj_angular_output (T.P, T.dir)? todeg_coord (T.P, T.dir, co): co;
     if (T.verbosity > 3)
         printf ("GOT      %.12f  %.12f  %.12f  %.12f\n", co.v[0],co.v[1],co.v[2],co.v[3]);
@@ -856,6 +916,7 @@ Tell GIE what to expect, when transforming the ACCEPTed input
 
     if (d > T.tolerance)
         return expect_message (d, args);
+    succs++;
 
     another_success ();
     return 0;
@@ -893,7 +954,20 @@ fprintf (T.fout, "%s\n", args);
 
 
 
+/*****************************************************************************/
+static int skip (const char *args) {
+/*****************************************************************************
+Indicate that the remaining material should be skipped. Mostly for debugging.
+******************************************************************************/
+    T.skip = 1;
+    (void) args;
+    return 0;
+}
+
+
 static int dispatch (const char *cmnd, const char *args) {
+    if (T.skip)
+        return SKIP;
     if  (0==strcmp (cmnd, "operation")) return  operation ((char *) args);
     if  (0==strcmp (cmnd, "accept"))    return  accept    (args);
     if  (0==strcmp (cmnd, "expect"))    return  expect    (args);
@@ -905,6 +979,7 @@ static int dispatch (const char *cmnd, const char *args) {
     if  (0==strcmp (cmnd, "ignore"))    return  ignore    (args);
     if  (0==strcmp (cmnd, "builtins"))  return  builtins  (args);
     if  (0==strcmp (cmnd, "echo"))      return  echo      (args);
+    if  (0==strcmp (cmnd, "skip"))      return  skip      (args);
 
     return 0;
 }
@@ -1223,6 +1298,8 @@ static int nextline (ffio *G) {
 Read next line of input file. Returns 1 on success, 0 on failure.
 ****************************************************************************************/
     G->next_args[0] = 0;
+    if (T.skip)
+        return 0;
     if (0==fgets (G->next_args, (int) G->next_args_size - 1, G->f))
         return 0;
     if (feof (G->f))
@@ -1570,7 +1647,7 @@ static int pj_cart_selftest (void) {
 
     /* Forward projection: Ellipsoidal-to-3D-Cartesian */
     dist = proj_roundtrip (P, PJ_FWD, 1, &a);
-    if (dist > 1e-12)
+    if (dist > 1e-9)
         return 8;
 
     /* Test at the South Pole */
@@ -1582,7 +1659,7 @@ static int pj_cart_selftest (void) {
 
     /* Forward projection: Ellipsoidal-to-3D-Cartesian */
     dist = proj_roundtrip (P, PJ_FWD, 1, &a);
-    if (dist > 1e-12)
+    if (dist > 1e-9)
         return 9;