From 51967aa03df0bbe134c38475a2633430dddefc82 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Tue, 3 Apr 2018 17:45:23 +0200 Subject: Add --skip-lines option to cct With this it is possible to skip the header when transforming coordinates from a file. --- src/cct.c | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/cct.c b/src/cct.c index 6b8f35f3..2952f2b7 100644 --- a/src/cct.c +++ b/src/cct.c @@ -105,6 +105,7 @@ static const char usage[] = { " -z value Provide a fixed z value for all input data (e.g. -z 0)\n" " -t value Provide a fixed t value for all input data (e.g. -t 0)\n" " -I Do the inverse transformation\n" + " -s n Skip n first lines of a infile\n" " -v Verbose: Provide non-essential informational output.\n" " Repeat -v for more verbosity (e.g. -vv)\n" "--------------------------------------------------------------------------------\n" @@ -116,6 +117,7 @@ static const char usage[] = { " --time Alias for -t\n" " --verbose Alias for -v\n" " --inverse Alias for -I\n" + " --skip-lines Alias for -s\n" " --help Alias for -h\n" " --version Print version number\n" "--------------------------------------------------------------------------------\n" @@ -156,13 +158,19 @@ int main(int argc, char **argv) { OPTARGS *o; FILE *fout = stdout; char *buf; - int nfields = 4, direction = 1, verbose; + int nfields = 4, direction = 1, skip_lines = 0, verbose; double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL; int columns_xyzt[] = {1, 2, 3, 4}; const char *longflags[] = {"v=verbose", "h=help", "I=inverse", "version", 0}; - const char *longkeys[] = {"o=output", "c=columns", "z=height", "t=time", 0}; - - o = opt_parse (argc, argv, "hvI", "cozt", longflags, longkeys); + const char *longkeys[] = { + "o=output", + "c=columns", + "z=height", + "t=time", + "s=skip-lines", + 0}; + + o = opt_parse (argc, argv, "hvI", "cozts", longflags, longkeys); if (0==o) return 0; @@ -202,6 +210,10 @@ int main(int argc, char **argv) { nfields--; } + if (opt_given (o, "s")) { + skip_lines = atoi (opt_arg(o, "s")); + } + if (opt_given (o, "c")) { /* cppcheck-suppress invalidscanf */ int ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+2, columns_xyzt+3); @@ -265,6 +277,10 @@ int main(int argc, char **argv) { continue; } point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time); + if (skip_lines > 0) { + skip_lines--; + continue; + } if (HUGE_VAL==point.xyzt.x) { char *c = column (buf, 1); -- cgit v1.2.3 From ac78787bf3bfa15dabbbad0042ca762020542365 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 15:35:08 +0200 Subject: Pipeline: make sure geocentric/cartesian space transform is done with original ellipsoid parameters (before any projection mess with them) --- src/PJ_pipeline.c | 2 ++ src/proj_4D_api.c | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index b60e05a8..aa7d76f8 100644 --- a/src/PJ_pipeline.c +++ b/src/PJ_pipeline.c @@ -306,6 +306,8 @@ static void set_ellipsoid(PJ *P) { /* the PJ you provided". */ proj_errno_reset (P); } + P->a_orig = P->a; + P->es_orig = P->es; pj_calc_ellipsoid_params (P, P->a, P->es); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index 11a56ac0..aed4685b 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -485,7 +485,7 @@ Returns 1 on success, 0 on failure if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6]) { /* If the current ellipsoid is not WGS84, then make sure the */ /* change in ellipsoid is still done. */ - if (!(fabs(P->a - 6378137.0) < 1e-8 && fabs(P->f - 1./ 298.257223563) < 1e-15)) { + if (!(fabs(P->a_orig - 6378137.0) < 1e-8 && fabs(P->es_orig - 0.0066943799901413) < 1e-15)) { do_cart = 1; } break; @@ -512,7 +512,7 @@ Returns 1 on success, 0 on failure /* geocentric/cartesian space or we need to do a Helmert transform. */ if (P->is_geocent || P->helmert || do_cart) { char def[150]; - sprintf (def, "break_cs2cs_recursion proj=cart a=%40.20g f=%40.20g", P->a, P->f); + sprintf (def, "break_cs2cs_recursion proj=cart a=%40.20g es=%40.20g", P->a_orig, P->es_orig); Q = proj_create (P->ctx, def); if (0==Q) return 0; -- cgit v1.2.3 From 3a2ddb6c6efbccf388ea89e177ca51fd25946ecf Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 13:41:12 +0200 Subject: Add webmerc projection --- src/PJ_merc.c | 17 +++++++++++++++++ src/pj_list.h | 1 + 2 files changed, 18 insertions(+) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 18303c45..a3e5e846 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -1,8 +1,10 @@ #define PJ_LIB__ +#include "proj_internal.h" #include "proj.h" #include "projects.h" PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; +PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 @@ -76,3 +78,18 @@ PJ *PROJECTION(merc) { return P; } +PJ *PROJECTION(webmerc) { + + /* Overriding k_0, lat_0 and lon_0 with fixed parameters */ + P->k0 = 1.0; + P->phi0 = 0.0; + P->lam0 = 0.0; + + P->b = P->a; + /* Clean up the ellipsoidal parameters to reflect the sphere */ + P->es = P->e = P->f = 0; + pj_calc_ellipsoid_params (P, P->a, 0); + P->inv = s_inverse; + P->fwd = s_forward; + return P; +} diff --git a/src/pj_list.h b/src/pj_list.h index 09a66034..440f7972 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -156,6 +156,7 @@ PROJ_HEAD(wag4, "Wagner IV") PROJ_HEAD(wag5, "Wagner V") PROJ_HEAD(wag6, "Wagner VI") PROJ_HEAD(wag7, "Wagner VII") +PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") PROJ_HEAD(weren, "Werenskiold I") PROJ_HEAD(wink1, "Winkel I") PROJ_HEAD(wink2, "Winkel II") -- cgit v1.2.3 From f1ea5041d5049498e3cd88f9a215ebb0f4b86a34 Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Wed, 11 Apr 2018 17:56:08 +0200 Subject: Enhance the precision of Spherical Mercator projection near the equator This uses a linear approximation of tan(x+pi/4) for better precision at small latitudes. As a result points of latitude 0 maintain a 0 Y coordinate and 0,0 is transformed to 0,0 --- src/PJ_merc.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index a3e5e846..3801b828 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -8,6 +8,13 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 +static double _tan_near_fort_pi(double x) { + if (fabs(x) <= __DBL_EPSILON__) { + return 2*x + 1.0; + } + return tan(M_FORTPI + x); +} + static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ XY xy = {0.0,0.0}; if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { @@ -27,7 +34,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } xy.x = P->k0 * lp.lam; - xy.y = P->k0 * log(tan(M_FORTPI + .5 * lp.phi)); + xy.y = P->k0 * log(_tan_near_fort_pi(.5 * lp.phi)); return xy; } -- cgit v1.2.3 From bd2e1363b021e1a1c6a4ebbc0aa55e8dc98e2587 Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 08:50:59 +0200 Subject: Fix: use proper double epsilon declaration --- src/PJ_merc.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 3801b828..2c89fbc6 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -2,6 +2,7 @@ #include "proj_internal.h" #include "proj.h" #include "projects.h" +#include PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; @@ -9,7 +10,7 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 static double _tan_near_fort_pi(double x) { - if (fabs(x) <= __DBL_EPSILON__) { + if (fabs(x) <= DBL_EPSILON) { return 2*x + 1.0; } return tan(M_FORTPI + x); -- cgit v1.2.3 From f2b3604a94349697828baa7de63a30a97ac15b2b Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 09:47:43 +0200 Subject: Use log1p in forward spherical mercator --- src/PJ_merc.c | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 2c89fbc6..0bf98625 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -9,11 +9,31 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 -static double _tan_near_fort_pi(double x) { +#if !defined(HAVE_C99_MATH) +#define HAVE_C99_MATH 0 +#endif + +#if HAVE_C99_MATH +#define log1px log1p +#else +static double log1px(double x) { + volatile double + y = 1 + x, + z = y - 1; + /* Here's the explanation for this magic: y = 1 + z, exactly, and z + * approx x, thus log(y)/z (which is nearly constant near z = 0) returns + * a good approximation to the true log(1 + x)/x. The multiplication x * + * (log(y)/z) introduces little additional error. */ + return z == 0 ? x : x * log(y) / z; +} +#endif + +static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ if (fabs(x) <= DBL_EPSILON) { - return 2*x + 1.0; + /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */ + return log1px(x); } - return tan(M_FORTPI + x); + return log(tan(M_FORTPI + .5 * x)); } static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ @@ -35,7 +55,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } xy.x = P->k0 * lp.lam; - xy.y = P->k0 * log(_tan_near_fort_pi(.5 * lp.phi)); + xy.y = P->k0 * logtanpfpim1(lp.phi); return xy; } -- cgit v1.2.3 From 587d454391b867f5338c91de356ad40404efcc7f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 15:35:08 +0200 Subject: Pipeline: make sure geocentric/cartesian space transform is done with original ellipsoid parameters (before any projection mess with them) --- src/PJ_pipeline.c | 2 ++ src/proj_4D_api.c | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c index b60e05a8..aa7d76f8 100644 --- a/src/PJ_pipeline.c +++ b/src/PJ_pipeline.c @@ -306,6 +306,8 @@ static void set_ellipsoid(PJ *P) { /* the PJ you provided". */ proj_errno_reset (P); } + P->a_orig = P->a; + P->es_orig = P->es; pj_calc_ellipsoid_params (P, P->a, P->es); diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index 11a56ac0..aed4685b 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -485,7 +485,7 @@ Returns 1 on success, 0 on failure if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6]) { /* If the current ellipsoid is not WGS84, then make sure the */ /* change in ellipsoid is still done. */ - if (!(fabs(P->a - 6378137.0) < 1e-8 && fabs(P->f - 1./ 298.257223563) < 1e-15)) { + if (!(fabs(P->a_orig - 6378137.0) < 1e-8 && fabs(P->es_orig - 0.0066943799901413) < 1e-15)) { do_cart = 1; } break; @@ -512,7 +512,7 @@ Returns 1 on success, 0 on failure /* geocentric/cartesian space or we need to do a Helmert transform. */ if (P->is_geocent || P->helmert || do_cart) { char def[150]; - sprintf (def, "break_cs2cs_recursion proj=cart a=%40.20g f=%40.20g", P->a, P->f); + sprintf (def, "break_cs2cs_recursion proj=cart a=%40.20g es=%40.20g", P->a_orig, P->es_orig); Q = proj_create (P->ctx, def); if (0==Q) return 0; -- cgit v1.2.3 From a4f4228d055aac7bb66f77b10fe5128e69e9e09a Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 7 Apr 2018 13:41:12 +0200 Subject: Add webmerc projection --- src/PJ_merc.c | 17 +++++++++++++++++ src/pj_list.h | 1 + 2 files changed, 18 insertions(+) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 18303c45..a3e5e846 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -1,8 +1,10 @@ #define PJ_LIB__ +#include "proj_internal.h" #include "proj.h" #include "projects.h" PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; +PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 @@ -76,3 +78,18 @@ PJ *PROJECTION(merc) { return P; } +PJ *PROJECTION(webmerc) { + + /* Overriding k_0, lat_0 and lon_0 with fixed parameters */ + P->k0 = 1.0; + P->phi0 = 0.0; + P->lam0 = 0.0; + + P->b = P->a; + /* Clean up the ellipsoidal parameters to reflect the sphere */ + P->es = P->e = P->f = 0; + pj_calc_ellipsoid_params (P, P->a, 0); + P->inv = s_inverse; + P->fwd = s_forward; + return P; +} diff --git a/src/pj_list.h b/src/pj_list.h index 09a66034..440f7972 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -156,6 +156,7 @@ PROJ_HEAD(wag4, "Wagner IV") PROJ_HEAD(wag5, "Wagner V") PROJ_HEAD(wag6, "Wagner VI") PROJ_HEAD(wag7, "Wagner VII") +PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") PROJ_HEAD(weren, "Werenskiold I") PROJ_HEAD(wink1, "Winkel I") PROJ_HEAD(wink2, "Winkel II") -- cgit v1.2.3 From 647637f81247d341b8865d5d0eda1cad3c8cec8a Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Wed, 11 Apr 2018 17:56:08 +0200 Subject: Enhance the precision of Spherical Mercator projection near the equator This uses a linear approximation of tan(x+pi/4) for better precision at small latitudes. As a result points of latitude 0 maintain a 0 Y coordinate and 0,0 is transformed to 0,0 --- src/PJ_merc.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index a3e5e846..3801b828 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -8,6 +8,13 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 +static double _tan_near_fort_pi(double x) { + if (fabs(x) <= __DBL_EPSILON__) { + return 2*x + 1.0; + } + return tan(M_FORTPI + x); +} + static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ XY xy = {0.0,0.0}; if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) { @@ -27,7 +34,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } xy.x = P->k0 * lp.lam; - xy.y = P->k0 * log(tan(M_FORTPI + .5 * lp.phi)); + xy.y = P->k0 * log(_tan_near_fort_pi(.5 * lp.phi)); return xy; } -- cgit v1.2.3 From 74030eb65a21e43b08b32239e5ec7bb92b30eb77 Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 08:50:59 +0200 Subject: Fix: use proper double epsilon declaration --- src/PJ_merc.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 3801b828..2c89fbc6 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -2,6 +2,7 @@ #include "proj_internal.h" #include "proj.h" #include "projects.h" +#include PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; @@ -9,7 +10,7 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 static double _tan_near_fort_pi(double x) { - if (fabs(x) <= __DBL_EPSILON__) { + if (fabs(x) <= DBL_EPSILON) { return 2*x + 1.0; } return tan(M_FORTPI + x); -- cgit v1.2.3 From 60c709f19a5a582c5f3f6ee53af9fee2ddb7fe50 Mon Sep 17 00:00:00 2001 From: Javier Goizueta Date: Thu, 12 Apr 2018 09:47:43 +0200 Subject: Use log1p in forward spherical mercator --- src/PJ_merc.c | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 2c89fbc6..0bf98625 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -9,11 +9,31 @@ PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 -static double _tan_near_fort_pi(double x) { +#if !defined(HAVE_C99_MATH) +#define HAVE_C99_MATH 0 +#endif + +#if HAVE_C99_MATH +#define log1px log1p +#else +static double log1px(double x) { + volatile double + y = 1 + x, + z = y - 1; + /* Here's the explanation for this magic: y = 1 + z, exactly, and z + * approx x, thus log(y)/z (which is nearly constant near z = 0) returns + * a good approximation to the true log(1 + x)/x. The multiplication x * + * (log(y)/z) introduces little additional error. */ + return z == 0 ? x : x * log(y) / z; +} +#endif + +static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ if (fabs(x) <= DBL_EPSILON) { - return 2*x + 1.0; + /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */ + return log1px(x); } - return tan(M_FORTPI + x); + return log(tan(M_FORTPI + .5 * x)); } static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ @@ -35,7 +55,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } xy.x = P->k0 * lp.lam; - xy.y = P->k0 * log(_tan_near_fort_pi(.5 * lp.phi)); + xy.y = P->k0 * logtanpfpim1(lp.phi); return xy; } -- cgit v1.2.3 From f9d81882c473310e5a1785bdbee06d8658998051 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Tue, 3 Apr 2018 17:45:23 +0200 Subject: Add --skip-lines option to cct With this it is possible to skip the header when transforming coordinates from a file. --- src/cct.c | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/cct.c b/src/cct.c index 6b8f35f3..2952f2b7 100644 --- a/src/cct.c +++ b/src/cct.c @@ -105,6 +105,7 @@ static const char usage[] = { " -z value Provide a fixed z value for all input data (e.g. -z 0)\n" " -t value Provide a fixed t value for all input data (e.g. -t 0)\n" " -I Do the inverse transformation\n" + " -s n Skip n first lines of a infile\n" " -v Verbose: Provide non-essential informational output.\n" " Repeat -v for more verbosity (e.g. -vv)\n" "--------------------------------------------------------------------------------\n" @@ -116,6 +117,7 @@ static const char usage[] = { " --time Alias for -t\n" " --verbose Alias for -v\n" " --inverse Alias for -I\n" + " --skip-lines Alias for -s\n" " --help Alias for -h\n" " --version Print version number\n" "--------------------------------------------------------------------------------\n" @@ -156,13 +158,19 @@ int main(int argc, char **argv) { OPTARGS *o; FILE *fout = stdout; char *buf; - int nfields = 4, direction = 1, verbose; + int nfields = 4, direction = 1, skip_lines = 0, verbose; double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL; int columns_xyzt[] = {1, 2, 3, 4}; const char *longflags[] = {"v=verbose", "h=help", "I=inverse", "version", 0}; - const char *longkeys[] = {"o=output", "c=columns", "z=height", "t=time", 0}; - - o = opt_parse (argc, argv, "hvI", "cozt", longflags, longkeys); + const char *longkeys[] = { + "o=output", + "c=columns", + "z=height", + "t=time", + "s=skip-lines", + 0}; + + o = opt_parse (argc, argv, "hvI", "cozts", longflags, longkeys); if (0==o) return 0; @@ -202,6 +210,10 @@ int main(int argc, char **argv) { nfields--; } + if (opt_given (o, "s")) { + skip_lines = atoi (opt_arg(o, "s")); + } + if (opt_given (o, "c")) { /* cppcheck-suppress invalidscanf */ int ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+2, columns_xyzt+3); @@ -265,6 +277,10 @@ int main(int argc, char **argv) { continue; } point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time); + if (skip_lines > 0) { + skip_lines--; + continue; + } if (HUGE_VAL==point.xyzt.x) { char *c = column (buf, 1); -- cgit v1.2.3 From 2f082b70cbdafdea49bb123e027406089e7ad65b Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Thu, 5 Apr 2018 19:25:31 +0000 Subject: Move logging functions to proj.h API --- src/pj_internal.c | 9 ++------- src/proj.h | 15 ++++++++++++++- src/proj_internal.h | 16 ---------------- 3 files changed, 16 insertions(+), 24 deletions(-) (limited to 'src') diff --git a/src/pj_internal.c b/src/pj_internal.c index 4da47051..891e0b9f 100644 --- a/src/pj_internal.c +++ b/src/pj_internal.c @@ -364,11 +364,6 @@ to that context. return; } - - - - - /* logging */ /* pj_vlog resides in pj_log.c and relates to pj_log as vsprintf relates to sprintf */ @@ -376,11 +371,11 @@ void pj_vlog( projCtx ctx, int level, const char *fmt, va_list args ); /***************************************************************************************/ -enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level) { +PJ_LOG_LEVEL proj_log_level (PJ_CONTEXT *ctx, PJ_LOG_LEVEL log_level) { /**************************************************************************************** Set logging level 0-3. Higher number means more debug info. 0 turns it off ****************************************************************************************/ - enum proj_log_level previous; + PJ_LOG_LEVEL previous; if (0==ctx) ctx = pj_get_default_ctx(); if (0==ctx) diff --git a/src/proj.h b/src/proj.h index 7bc9b10e..3d64ed1a 100644 --- a/src/proj.h +++ b/src/proj.h @@ -268,6 +268,17 @@ struct PJ_INIT_INFO { char lastupdate[16]; /* Date of last update in YYYY-MM-DD format */ }; +typedef enum PJ_LOG_LEVEL { + PJ_LOG_NONE = 0, + PJ_LOG_ERROR = 1, + PJ_LOG_DEBUG = 2, + PJ_LOG_TRACE = 3, + PJ_LOG_TELL = 4, + PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */ + PJ_LOG_DEBUG_MINOR = 3 /* for proj_api.h compatibility */ +} PJ_LOG_LEVEL; + +typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *); /* The context type - properly namespaced synonym for projCtx */ @@ -342,7 +353,6 @@ double proj_xyz_dist (PJ_COORD a, PJ_COORD b); PJ_COORD proj_geod (const PJ *P, PJ_COORD a, PJ_COORD b); - /* Set or read error level */ int proj_context_errno (PJ_CONTEXT *ctx); int proj_errno (const PJ *P); @@ -351,6 +361,9 @@ int proj_errno_reset (const PJ *P); int proj_errno_restore (const PJ *P, int err); const char* proj_errno_string (int err); +PJ_LOG_LEVEL proj_log_level (PJ_CONTEXT *ctx, PJ_LOG_LEVEL level); +void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf); + /* Scaling and angular distortion factors */ PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp); diff --git a/src/proj_internal.h b/src/proj_internal.h index 3f6ccde0..1aca9201 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -97,25 +97,9 @@ double proj_vgrid_value(PJ *P, PJ_LP lp); PJ_LP proj_hgrid_value(PJ *P, PJ_LP lp); PJ_LP proj_hgrid_apply(PJ *P, PJ_LP lp, PJ_DIRECTION direction); -/* High level functionality for handling thread contexts */ -enum proj_log_level { - PJ_LOG_NONE = 0, - PJ_LOG_ERROR = 1, - PJ_LOG_DEBUG = 2, - PJ_LOG_TRACE = 3, - PJ_LOG_TELL = 4, - PJ_LOG_DEBUG_MAJOR = 2, /* for proj_api.h compatibility */ - PJ_LOG_DEBUG_MINOR = 3 /* for proj_api.h compatibility */ -}; - -/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */ -enum proj_log_level proj_log_level (PJ_CONTEXT *ctx, enum proj_log_level log_level); -typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *); - void proj_log_error (PJ *P, const char *fmt, ...); void proj_log_debug (PJ *P, const char *fmt, ...); void proj_log_trace (PJ *P, const char *fmt, ...); -void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf); int pj_ellipsoid (PJ *); void pj_inherit_ellipsoid_def (const PJ *src, PJ *dst); -- cgit v1.2.3 From c3e346b7b61bd18b30de9f09ad0b1b507210ef5f Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Thu, 5 Apr 2018 19:26:07 +0000 Subject: Use PROJ logging facility in cct --- src/cct.c | 98 +++++++++++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 71 insertions(+), 27 deletions(-) (limited to 'src') diff --git a/src/cct.c b/src/cct.c index 2952f2b7..c3d94ed2 100644 --- a/src/cct.c +++ b/src/cct.c @@ -76,6 +76,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26 #include #include #include +#include #include "proj.h" #include "proj_internal.h" @@ -86,13 +87,14 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-25/2017-10-26 /* Prototypes for functions in proj_strtod.c */ double proj_strtod(const char *str, char **endptr); double proj_atof(const char *str); +static void logger(void *data, int level, const char *msg); +static void print(PJ_LOG_LEVEL verbosity, const char *fmt, ...); /* Prototypes from functions in this file */ char *column (char *buf, int n); PJ_COORD parse_input_line (char *buf, int *columns, double fixed_height, double fixed_time); - static const char usage[] = { "--------------------------------------------------------------------------------\n" "Usage: %s [-options]... [+operator_specs]... infile...\n" @@ -151,12 +153,57 @@ static const char usage[] = { "--------------------------------------------------------------------------------\n" }; + +static void logger(void *data, int level, const char *msg) { + FILE *stream; + int log_tell = proj_log_level(PJ_DEFAULT_CTX, PJ_LOG_TELL); + + stream = (FILE *) data; + + /* if we use PJ_LOG_NONE we always want to print stuff to stream */ + if (level == PJ_LOG_NONE) { + fprintf(stream, "%s", msg); + return; + } + + /* should always print to stderr if level == PJ_LOG_ERROR */ + if (level == PJ_LOG_ERROR) { + fprintf(stderr, "%s", msg); + return; + } + + /* otherwise only print if log level set by user is high enough */ + if (level <= log_tell) + fprintf(stream, "%s", msg); +} + +FILE *fout; + +static void print(PJ_LOG_LEVEL log_level, const char *fmt, ...) { + + va_list args; + char *msg_buf; + + va_start( args, fmt ); + + msg_buf = (char *) malloc(100000); + if( msg_buf == NULL ) + return; + + vsprintf( msg_buf, fmt, args ); + + logger((void *) fout, log_level, msg_buf); + + va_end( args ); + free( msg_buf ); +} + + int main(int argc, char **argv) { PJ *P; PJ_COORD point; PJ_PROJ_INFO info; OPTARGS *o; - FILE *fout = stdout; char *buf; int nfields = 4, direction = 1, skip_lines = 0, verbose; double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL; @@ -179,26 +226,26 @@ int main(int argc, char **argv) { return 0; } - direction = opt_given (o, "I")? -1: 1; - verbose = opt_given (o, "v"); + + verbose = MIN(opt_given (o, "v"), 3); /* log level can't be larger than 3 */ + proj_log_level (PJ_DEFAULT_CTX, verbose); + proj_log_func (PJ_DEFAULT_CTX, (void *) fout, logger); if (opt_given (o, "version")) { - fprintf (stdout, "%s: %s\n", o->progname, pj_get_release ()); + print (PJ_LOG_NONE, "%s: %s\n", o->progname, pj_get_release ()); return 0; } if (opt_given (o, "o")) fout = fopen (opt_arg (o, "output"), "wt"); if (0==fout) { - fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output")); + print (PJ_LOG_ERROR, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output")); free (o); return 1; } - if (verbose > 3) - fprintf (fout, "%s: Running in very verbose mode\n", o->progname); - + print (PJ_LOG_TRACE, "%s: Running in very verbose mode\n", o->progname); if (opt_given (o, "z")) { fixed_z = proj_atof (opt_arg (o, "z")); @@ -218,7 +265,7 @@ int main(int argc, char **argv) { /* cppcheck-suppress invalidscanf */ int ncols = sscanf (opt_arg (o, "c"), "%d,%d,%d,%d", columns_xyzt, columns_xyzt+1, columns_xyzt+2, columns_xyzt+3); if (ncols != nfields) { - fprintf (stderr, "%s: Too few input columns given: '%s'\n", o->progname, opt_arg (o, "c")); + print (PJ_LOG_ERROR, "%s: Too few input columns given: '%s'\n", o->progname, opt_arg (o, "c")); free (o); if (stdout != fout) fclose (fout); @@ -229,7 +276,7 @@ int main(int argc, char **argv) { /* Setup transformation */ P = proj_create_argv (0, o->pargc, o->pargv); if ((0==P) || (0==o->pargc)) { - fprintf (stderr, "%s: Bad transformation arguments - (%s)\n '%s -h' for help\n", + print (PJ_LOG_ERROR, "%s: Bad transformation arguments - (%s)\n '%s -h' for help\n", o->progname, pj_strerrno (proj_errno(P)), o->progname); free (o); if (stdout != fout) @@ -237,15 +284,13 @@ int main(int argc, char **argv) { return 1; } - if (verbose > 4) { - info = proj_pj_info (P); - fprintf (stdout, "Final: %s argc=%d pargc=%d\n", info.definition, argc, o->pargc); - } + info = proj_pj_info (P); + print (PJ_LOG_TRACE, "Final: %s argc=%d pargc=%d\n", info.definition, argc, o->pargc); if (direction==-1) { /* fail if an inverse operation is not available */ - if (!proj_pj_info(P).has_inverse) { - fprintf (stderr, "Inverse operation not available\n"); + if (!info.has_inverse) { + print (PJ_LOG_ERROR, "Inverse operation not available\n"); if (stdout != fout) fclose (fout); return 1; @@ -258,7 +303,7 @@ int main(int argc, char **argv) { /* Allocate input buffer */ buf = calloc (1, 10000); if (0==buf) { - fprintf (stderr, "%s: Out of memory\n", o->progname); + print (PJ_LOG_ERROR, "%s: Out of memory\n", o->progname); pj_free (P); free (o); if (stdout != fout) @@ -273,7 +318,7 @@ int main(int argc, char **argv) { void *ret = fgets (buf, 10000, o->input); opt_eof_handler (o); if (0==ret) { - fprintf (stderr, "Read error in record %d\n", (int) o->record_index); + print (PJ_LOG_ERROR, "Read error in record %d\n", (int) o->record_index); continue; } point = parse_input_line (buf, columns_xyzt, fixed_z, fixed_time); @@ -286,14 +331,13 @@ int main(int argc, char **argv) { /* if it's a comment or blank line, we reflect it */ if (c && ((*c=='\0') || (*c=='#'))) { - fprintf (fout, "%s\n", buf); + print (PJ_LOG_NONE, "%s\n", buf); continue; } /* otherwise, it must be a syntax error */ - fprintf (fout, "# Record %d UNREADABLE: %s", (int) o->record_index, buf); - if (verbose) - fprintf (stderr, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o)); + print (PJ_LOG_NONE, "# Record %d UNREADABLE: %s", (int) o->record_index, buf); + print (PJ_LOG_ERROR, "%s: Could not parse file '%s' line %d\n", o->progname, opt_filename (o), opt_record (o)); continue; } @@ -306,8 +350,8 @@ int main(int argc, char **argv) { if (HUGE_VAL==point.xyzt.x) { /* transformation error */ - fprintf (fout, "# Record %d TRANSFORMATION ERROR: %s (%s)", - (int) o->record_index, buf, pj_strerrno (proj_errno(P))); + print (PJ_LOG_NONE, "# Record %d TRANSFORMATION ERROR: %s (%s)", + (int) o->record_index, buf, pj_strerrno (proj_errno(P))); proj_errno_restore (P, err); continue; } @@ -317,10 +361,10 @@ int main(int argc, char **argv) { if (proj_angular_output (P, direction)) { point.lpzt.lam = proj_todeg (point.lpzt.lam); point.lpzt.phi = proj_todeg (point.lpzt.phi); - fprintf (fout, "%14.10f %14.10f %12.4f %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); + print (PJ_LOG_NONE, "%14.10f %14.10f %12.4f %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); } else - fprintf (fout, "%13.4f %13.4f %12.4f %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); + print (PJ_LOG_NONE, "%13.4f %13.4f %12.4f %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t); } if (stdout != fout) -- cgit v1.2.3 From be27d12a028d07ed656080004c949c0da7a6012e Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Fri, 13 Apr 2018 18:11:29 +0200 Subject: Always ignore out-commented lines in cct Previous to this commit cct would return the following $ cct -c 2,3,4 -t 0 -I +proj=cart +ellps=GRS80 BLAH 3579685.56545 508396.50343 5236837.50646 8.0832413787 55.5578176654 99.9833 0.0000 8.0832413787 55.5578176654 99.9833 0.0000 where the second input should not be parsed as a valid coordinate. With this commit that no longer happens and the input is returned verbatim back to the user. Closes #932 --- src/cct.c | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) (limited to 'src') diff --git a/src/cct.c b/src/cct.c index 2952f2b7..00fc2127 100644 --- a/src/cct.c +++ b/src/cct.c @@ -271,6 +271,7 @@ int main(int argc, char **argv) { while (opt_input_loop (o, optargs_file_format_text)) { int err; void *ret = fgets (buf, 10000, o->input); + char *c = column (buf, 1); opt_eof_handler (o); if (0==ret) { fprintf (stderr, "Read error in record %d\n", (int) o->record_index); @@ -281,15 +282,14 @@ int main(int argc, char **argv) { skip_lines--; continue; } - if (HUGE_VAL==point.xyzt.x) { - char *c = column (buf, 1); - /* if it's a comment or blank line, we reflect it */ - if (c && ((*c=='\0') || (*c=='#'))) { - fprintf (fout, "%s\n", buf); - continue; - } + /* if it's a comment or blank line, we reflect it */ + if (c && ((*c=='\0') || (*c=='#'))) { + fprintf (fout, "%s", buf); + continue; + } + if (HUGE_VAL==point.xyzt.x) { /* otherwise, it must be a syntax error */ fprintf (fout, "# Record %d UNREADABLE: %s", (int) o->record_index, buf); if (verbose) -- cgit v1.2.3 From 299135153e33869157ecd432dd0194743e959170 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Fri, 13 Apr 2018 22:42:54 +0200 Subject: Fix segfault in deformation When transforming coordinates outside the grid model the deformation operation failed spectatularly. This is now fixed by checking that the coordinate is inside the grid. If it isn't an error is returned and a debugging log message is issued. Closes #934 --- src/PJ_deformation.c | 7 +++++++ src/pj_apply_gridshift.c | 7 ++++++- 2 files changed, 13 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c index d6376e40..58c0e4df 100644 --- a/src/PJ_deformation.c +++ b/src/PJ_deformation.c @@ -82,6 +82,7 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) { ********************************************************************************/ PJ_COORD geodetic, shift, temp; double sp, cp, sl, cl; + int previous_errno = proj_errno_reset(P); /* cartesian to geodetic */ geodetic.lpz = pj_inv3d(cartesian, P->opaque->cart); @@ -90,6 +91,10 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) { shift.lp = proj_hgrid_value(P, geodetic.lp); shift.enu.u = proj_vgrid_value(P, geodetic.lp); + if (proj_errno(P) == PJD_ERR_GRID_AREA) + proj_log_debug(P, "deformation: coordinate (%.3f, %.3f) outside deformation model", + proj_todeg(geodetic.lp.lam), proj_todeg(geodetic.lp.phi)); + /* grid values are stored as mm/yr, we need m/yr */ shift.xyz.x /= 1000; shift.xyz.y /= 1000; @@ -108,6 +113,8 @@ static XYZ get_grid_shift(PJ* P, XYZ cartesian) { shift.xyz = temp.xyz; + proj_errno_restore(P, previous_errno); + return shift.xyz; } diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c index c503ec0b..1ef39b20 100644 --- a/src/pj_apply_gridshift.c +++ b/src/pj_apply_gridshift.c @@ -306,13 +306,18 @@ int proj_hgrid_init(PJ* P, const char *grids) { /********************************************/ LP proj_hgrid_value(PJ *P, LP lp) { struct CTABLE *ct; - LP out; + LP out = proj_coord_error().lp; ct = find_ctable(P->ctx, lp, P->gridlist_count, P->gridlist); + if (ct == 0) { + pj_ctx_set_errno( P->ctx, PJD_ERR_GRID_AREA); + return out; + } /* normalize input to ll origin */ lp.lam -= ct->ll.lam; lp.phi -= ct->ll.phi; + lp.lam = adjlon(lp.lam - M_PI) + M_PI; out = nad_intr(lp, ct); -- cgit v1.2.3 From 837e1ecf6189e64326c5dc336bd9898f1d195468 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sat, 14 Apr 2018 12:41:30 +0200 Subject: Dont do false easting/northing on cartesian coords False easting and northing should only be applied to projected coordinates (PJ_IO_UNITS_PROJECTED). This commit removes the option of false easting/northing on operations suchs as helmert and deformation that both work on cartesian coordinates. --- src/pj_fwd.c | 3 --- src/pj_inv.c | 4 ---- 2 files changed, 7 deletions(-) (limited to 'src') diff --git a/src/pj_fwd.c b/src/pj_fwd.c index 2a064e58..347f8334 100644 --- a/src/pj_fwd.c +++ b/src/pj_fwd.c @@ -112,9 +112,6 @@ static PJ_COORD fwd_finalize (PJ *P, PJ_COORD coo) { coo = proj_trans (P->cart, PJ_FWD, coo); } - coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0); - coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0); - coo.xyz.z = P->fr_meter * (coo.xyz.z + P->z0); break; /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */ diff --git a/src/pj_inv.c b/src/pj_inv.c index f66fd7d1..327fef9f 100644 --- a/src/pj_inv.c +++ b/src/pj_inv.c @@ -98,10 +98,6 @@ static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) { /* de-scale and de-offset */ case PJ_IO_UNITS_CARTESIAN: - coo.xyz.x = P->to_meter * coo.xyz.x - P->x0; - coo.xyz.y = P->to_meter * coo.xyz.y - P->y0; - coo.xyz.z = P->to_meter * coo.xyz.z - P->z0; - if (P->is_geocent) coo = proj_trans (P->cart, PJ_INV, coo); -- cgit v1.2.3 From adce503f6597d7cf6921f31bfa2c66db0c07f319 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 16 Apr 2018 13:22:19 +0200 Subject: Silence complaints about missing tag --- src/gie.c | 1 + 1 file changed, 1 insertion(+) (limited to 'src') diff --git a/src/gie.c b/src/gie.c index 8c735b88..3831a8f1 100644 --- a/src/gie.c +++ b/src/gie.c @@ -964,6 +964,7 @@ Indicate that the remaining material should be skipped. Mostly for debugging. ******************************************************************************/ T.skip = 1; (void) args; + F->level = 2; /* Silence complaints about missing element */ return 0; } -- cgit v1.2.3 From 327a8f8b5a850167901a92747767e14ffe77c03b Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sat, 14 Apr 2018 16:47:53 +0200 Subject: Collect custom C99 math functions in proj_math.h We are relying more and more on C99 math functions. On C89 systems where those functions are not available our own custom versions of those functions are used instead. So far these has been spread across the code base. This commit gathers them in the same file and introduces the proj_math.h header. The build system checks for C99 math functions. If not found the proj_math.h header make sure that C99 functions are defined as their pj_ equivalent. Ideally proj_math.h is included instead of math.h. This removes the need for any checks against HAVE_C99_MATH in the code making it easier to read. For this commit the functions hypot, log1p and asinh has been taken care of. --- src/Makefile.am | 4 ++-- src/PJ_aea.c | 2 ++ src/PJ_aeqd.c | 1 + src/PJ_bipc.c | 1 + src/PJ_bonne.c | 2 ++ src/PJ_cart.c | 2 +- src/PJ_ccon.c | 1 + src/PJ_deformation.c | 1 + src/PJ_eqdc.c | 1 + src/PJ_geos.c | 1 + src/PJ_gnom.c | 1 + src/PJ_laea.c | 1 + src/PJ_lcc.c | 1 + src/PJ_merc.c | 23 ++---------------- src/PJ_mod_ster.c | 1 + src/PJ_nsper.c | 1 + src/PJ_oea.c | 1 + src/PJ_ortho.c | 1 + src/PJ_sconics.c | 1 + src/PJ_stere.c | 1 + src/PJ_sterea.c | 1 + src/PJ_tpeqd.c | 1 + src/gie.c | 1 + src/hypot.c | 36 ---------------------------- src/lib_proj.cmake | 2 ++ src/makefile.vc | 2 +- src/nad_cvt.c | 1 + src/pj_factors.c | 1 + src/pj_init.c | 1 + src/pj_math.c | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/proj_4D_api.c | 1 + src/proj_etmerc.c | 29 ++--------------------- src/proj_math.h | 53 +++++++++++++++++++++++++++++++++++++++++ src/projects.h | 5 ---- 34 files changed, 155 insertions(+), 93 deletions(-) delete mode 100644 src/hypot.c create mode 100644 src/pj_math.c create mode 100644 src/proj_math.h (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index f94e5f26..78ecfdf8 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -43,7 +43,7 @@ lib_LTLIBRARIES = libproj.la libproj_la_LDFLAGS = -no-undefined -version-info 13:0:0 libproj_la_SOURCES = \ - pj_list.h proj_internal.h\ + pj_list.h proj_internal.h proj_math.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_ccon.c\ @@ -83,7 +83,7 @@ libproj_la_SOURCES = \ pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \ geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \ jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \ - pj_strtod.c \ + pj_strtod.c pj_math.c\ \ proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \ PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \ diff --git a/src/PJ_aea.c b/src/PJ_aea.c index 3fe90524..640f1013 100644 --- a/src/PJ_aea.c +++ b/src/PJ_aea.c @@ -31,6 +31,8 @@ #include "proj.h" #include #include "projects.h" +#include "proj_math.h" + # define EPS10 1.e-10 # define TOL7 1.e-7 diff --git a/src/PJ_aeqd.c b/src/PJ_aeqd.c index 84554b91..5d8c3d38 100644 --- a/src/PJ_aeqd.c +++ b/src/PJ_aeqd.c @@ -30,6 +30,7 @@ #include "proj.h" #include #include "projects.h" +#include "proj_math.h" enum Mode { N_POLE = 0, diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c index 48fc93d0..2f60808d 100644 --- a/src/PJ_bipc.c +++ b/src/PJ_bipc.c @@ -2,6 +2,7 @@ #include "proj.h" #include #include "projects.h" +#include "proj_math.h" PROJ_HEAD(bipc, "Bipolar conic of western hemisphere") "\n\tConic Sph."; diff --git a/src/PJ_bonne.c b/src/PJ_bonne.c index 40bb28ee..d3d5b757 100644 --- a/src/PJ_bonne.c +++ b/src/PJ_bonne.c @@ -2,6 +2,8 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" + PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)") "\n\tConic Sph&Ell\n\tlat_1="; diff --git a/src/PJ_cart.c b/src/PJ_cart.c index ec273f1b..a4fd3254 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -43,8 +43,8 @@ #define PJ_LIB__ #include "proj_internal.h" #include "projects.h" +#include "proj_math.h" #include -#include #include PROJ_HEAD(cart, "Geodetic/cartesian conversions"); diff --git a/src/PJ_ccon.c b/src/PJ_ccon.c index 02b14388..aed59fb6 100644 --- a/src/PJ_ccon.c +++ b/src/PJ_ccon.c @@ -24,6 +24,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" #define EPS10 1e-10 diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c index 58c0e4df..5511eed4 100644 --- a/src/PJ_deformation.c +++ b/src/PJ_deformation.c @@ -55,6 +55,7 @@ grid-values in units of mm/year in ENU-space. #include #include "proj.h" #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" PROJ_HEAD(deformation, "Kinematic grid shift"); diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c index f36e0518..5f65e7d0 100644 --- a/src/PJ_eqdc.c +++ b/src/PJ_eqdc.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" struct pj_opaque { double phi1; diff --git a/src/PJ_geos.c b/src/PJ_geos.c index f0a9bc24..2988e62e 100644 --- a/src/PJ_geos.c +++ b/src/PJ_geos.c @@ -31,6 +31,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" struct pj_opaque { double h; diff --git a/src/PJ_gnom.c b/src/PJ_gnom.c index 8bb24b6a..b4bab0e2 100644 --- a/src/PJ_gnom.c +++ b/src/PJ_gnom.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph."; diff --git a/src/PJ_laea.c b/src/PJ_laea.c index 0422784b..bcf9c44d 100644 --- a/src/PJ_laea.c +++ b/src/PJ_laea.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell"; diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c index 863f6be9..34ed99cb 100644 --- a/src/PJ_lcc.c +++ b/src/PJ_lcc.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" PROJ_HEAD(lcc, "Lambert Conformal Conic") "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0"; diff --git a/src/PJ_merc.c b/src/PJ_merc.c index 0bf98625..9123eae1 100644 --- a/src/PJ_merc.c +++ b/src/PJ_merc.c @@ -1,6 +1,7 @@ #define PJ_LIB__ #include "proj_internal.h" #include "proj.h" +#include "proj_math.h" #include "projects.h" #include @@ -8,30 +9,10 @@ PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts="; PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Sph\n\t"; #define EPS10 1.e-10 - -#if !defined(HAVE_C99_MATH) -#define HAVE_C99_MATH 0 -#endif - -#if HAVE_C99_MATH -#define log1px log1p -#else -static double log1px(double x) { - volatile double - y = 1 + x, - z = y - 1; - /* Here's the explanation for this magic: y = 1 + z, exactly, and z - * approx x, thus log(y)/z (which is nearly constant near z = 0) returns - * a good approximation to the true log(1 + x)/x. The multiplication x * - * (log(y)/z) introduces little additional error. */ - return z == 0 ? x : x * log(y) / z; -} -#endif - static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */ if (fabs(x) <= DBL_EPSILON) { /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */ - return log1px(x); + return log1p(x); } return log(tan(M_FORTPI + .5 * x)); } diff --git a/src/PJ_mod_ster.c b/src/PJ_mod_ster.c index d807660c..5e6ce136 100644 --- a/src/PJ_mod_ster.c +++ b/src/PJ_mod_ster.c @@ -2,6 +2,7 @@ #define PJ_LIB__ #include #include "projects.h" +#include "proj_math.h" PROJ_HEAD(mil_os, "Miller Oblated Stereographic") "\n\tAzi(mod)"; PROJ_HEAD(lee_os, "Lee Oblated Stereographic") "\n\tAzi(mod)"; diff --git a/src/PJ_nsper.c b/src/PJ_nsper.c index 3c0b01f0..223cd75b 100644 --- a/src/PJ_nsper.c +++ b/src/PJ_nsper.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" enum Mode { N_POLE = 0, diff --git a/src/PJ_oea.c b/src/PJ_oea.c index 5112896f..2bfeffa8 100644 --- a/src/PJ_oea.c +++ b/src/PJ_oea.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" PROJ_HEAD(oea, "Oblated Equal Area") "\n\tMisc Sph\n\tn= m= theta="; diff --git a/src/PJ_ortho.c b/src/PJ_ortho.c index 7c27cd56..2b037819 100644 --- a/src/PJ_ortho.c +++ b/src/PJ_ortho.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph."; diff --git a/src/PJ_sconics.c b/src/PJ_sconics.c index 34a0f1fd..ce044c24 100644 --- a/src/PJ_sconics.c +++ b/src/PJ_sconics.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" enum Type { diff --git a/src/PJ_stere.c b/src/PJ_stere.c index 62e42460..82fd5c07 100644 --- a/src/PJ_stere.c +++ b/src/PJ_stere.c @@ -2,6 +2,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts="; PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth"; diff --git a/src/PJ_sterea.c b/src/PJ_sterea.c index 96b35609..eb4c9f2c 100644 --- a/src/PJ_sterea.c +++ b/src/PJ_sterea.c @@ -26,6 +26,7 @@ #define PJ_LIB__ #include #include "projects.h" +#include "proj_math.h" struct pj_opaque { diff --git a/src/PJ_tpeqd.c b/src/PJ_tpeqd.c index 21e3c5ed..87877ec1 100644 --- a/src/PJ_tpeqd.c +++ b/src/PJ_tpeqd.c @@ -1,6 +1,7 @@ #define PJ_LIB__ #include #include "proj.h" +#include "proj_math.h" #include "projects.h" diff --git a/src/gie.c b/src/gie.c index 8c735b88..10b5804e 100644 --- a/src/gie.c +++ b/src/gie.c @@ -114,6 +114,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2017-10-01/2017-10-08 #include "proj.h" #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" #include "optargpm.h" diff --git a/src/hypot.c b/src/hypot.c deleted file mode 100644 index 822c4595..00000000 --- a/src/hypot.c +++ /dev/null @@ -1,36 +0,0 @@ -/* hypot - sqrt(x * x + y * y) -** -** Because this was omitted from the ANSI standards, this version -** is included for those systems that do not include hypot as an -** extension to libm.a. Note: GNU version was not used because it -** was not properly coded to minimize potential overflow. -** -** The proper technique for determining hypot is to factor out the -** larger of the two terms, thus leaving a possible case of float -** overflow when max(x,y)*sqrt(2) > max machine value. This allows -** a wider range of numbers than the alternative of the sum of the -** squares < max machine value. For an Intel x87 IEEE double of -** approximately 1.8e308, only argument values > 1.27e308 are at -** risk of causing overflow. Whereas, not using this method limits -** the range to values less that 9.5e153 --- a considerable reduction -** in range! -*/ -extern double sqrt(double); - double -hypot(double x, double y) { - if ( x < 0.) - x = -x; - else if (x == 0.) - return (y < 0. ? -y : y); - if (y < 0.) - y = -y; - else if (y == 0.) - return (x); - if ( x < y ) { - x /= y; - return ( y * sqrt( 1. + x * x ) ); - } else { - y /= x; - return ( x * sqrt( 1. + y * y ) ); - } -} diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index c9e4d9e6..c0d0a6a9 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -199,6 +199,7 @@ SET(SRC_LIBPROJ_CORE pj_list.h pj_log.c pj_malloc.c + pj_math.c pj_mlfn.c pj_msfn.c pj_mutex.c @@ -218,6 +219,7 @@ SET(SRC_LIBPROJ_CORE pj_utils.c pj_zpoly1.c proj_mdist.c + proj_math.h proj_rouss.c rtodms.c vector1.c diff --git a/src/makefile.vc b/src/makefile.vc index ee89beb3..ef39b084 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -55,7 +55,7 @@ support = \ pj_utils.obj pj_gridlist.obj pj_gridinfo.obj \ proj_mdist.obj pj_mutex.obj pj_initcache.obj \ pj_ctx.obj pj_fileapi.obj pj_log.obj pj_apply_vgridshift.obj \ - pj_strtod.obj pj_internal.obj + pj_strtod.obj pj_internal.obj pj_math.obj pipeline = \ proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ diff --git a/src/nad_cvt.c b/src/nad_cvt.c index c913511f..d4330d54 100644 --- a/src/nad_cvt.c +++ b/src/nad_cvt.c @@ -1,5 +1,6 @@ #define PJ_LIB__ #include "projects.h" +#include "proj_math.h" #define MAX_ITERATIONS 10 #define TOL 1e-12 diff --git a/src/pj_factors.c b/src/pj_factors.c index e682a12d..373bd743 100644 --- a/src/pj_factors.c +++ b/src/pj_factors.c @@ -1,6 +1,7 @@ /* projection scale factors */ #define PJ_LIB__ #include "proj.h" +#include "proj_math.h" #include "projects.h" #include diff --git a/src/pj_init.c b/src/pj_init.c index 9b06ff28..bdaf64f7 100644 --- a/src/pj_init.c +++ b/src/pj_init.c @@ -37,6 +37,7 @@ #include #include "proj.h" #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" diff --git a/src/pj_math.c b/src/pj_math.c new file mode 100644 index 00000000..d5242636 --- /dev/null +++ b/src/pj_math.c @@ -0,0 +1,66 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Make C99 math functions available on C89 systems + * Author: Kristian Evers + * + ****************************************************************************** + * Copyright (c) 2018, Kristian Evers + * + * 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. + *****************************************************************************/ + +#include "proj_math.h" + +#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) + +/* Compute hypotenuse */ +double pj_hypot(double x, double y) { + x = fabs(x); + y = fabs(y); + if ( x < y ) { + x /= y; + return ( y * sqrt( 1. + x * x ) ); + } else { + y /= (x != 0.0 ? x : 1.0); + return ( x * sqrt( 1. + y * y ) ); + } +} + +/* Compute log(1+x) accurately */ +double pj_log1p(double x) { + volatile double + y = 1 + x, + z = y - 1; + /* Here's the explanation for this magic: y = 1 + z, exactly, and z + * approx x, thus log(y)/z (which is nearly constant near z = 0) returns + * a good approximation to the true log(1 + x)/x. The multiplication x * + * (log(y)/z) introduces little additional error. */ + return z == 0 ? x : x * log(y) / z; +} + +/* Compute asinh(x) accurately */ +double pj_asinh(double x) { + double y = fabs(x); /* Enforce odd parity */ + y = log1p(y * (1 + y/(hypot(1.0, y) + 1))); + return x > 0 ? y : (x < 0 ? -y : x); +} + +#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */ + + diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c index aed4685b..70110746 100644 --- a/src/proj_4D_api.c +++ b/src/proj_4D_api.c @@ -30,6 +30,7 @@ #include #include "proj.h" #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" #include "geodesic.h" diff --git a/src/proj_etmerc.c b/src/proj_etmerc.c index 99646d14..9a89e2b0 100644 --- a/src/proj_etmerc.c +++ b/src/proj_etmerc.c @@ -44,6 +44,7 @@ #include #include "proj.h" #include "projects.h" +#include "proj_math.h" struct pj_opaque { @@ -62,32 +63,6 @@ PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") #define PROJ_ETMERC_ORDER 6 - -#ifdef _GNU_SOURCE - inline -#endif -static double log1py(double x) { /* Compute log(1+x) accurately */ - volatile double - y = 1 + x, - z = y - 1; - /* Here's the explanation for this magic: y = 1 + z, exactly, and z - * approx x, thus log(y)/z (which is nearly constant near z = 0) returns - * a good approximation to the true log(1 + x)/x. The multiplication x * - * (log(y)/z) introduces little additional error. */ - return z == 0 ? x : x * log(y) / z; -} - - -#ifdef _GNU_SOURCE - inline -#endif -static double asinhy(double x) { /* Compute asinh(x) accurately */ - double y = fabs(x); /* Enforce odd parity */ - y = log1py(y * (1 + y/(hypot(1.0, y) + 1))); - return x < 0 ? -y : y; -} - - #ifdef _GNU_SOURCE inline #endif @@ -182,7 +157,7 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */ Ce = atan2 (sin_Ce*cos_Cn, hypot (sin_Cn, cos_Cn*cos_Ce)); /* compl. sph. N, E -> ell. norm. N, E */ - Ce = asinhy ( tan (Ce) ); /* Replaces: Ce = log(tan(FORTPI + Ce*0.5)); */ + Ce = asinh ( tan (Ce) ); /* Replaces: Ce = log(tan(FORTPI + Ce*0.5)); */ Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); Ce += dCe; if (fabs (Ce) <= 2.623395162778) { diff --git a/src/proj_math.h b/src/proj_math.h new file mode 100644 index 00000000..0d28cb30 --- /dev/null +++ b/src/proj_math.h @@ -0,0 +1,53 @@ +/****************************************************************************** + * Project: PROJ + * Purpose: Make C99 math functions available on C89 systems + * Author: Kristian Evers + * + ****************************************************************************** + * Copyright (c) 2018, Kristian Evers + * + * 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. + *****************************************************************************/ + +#include + +#ifndef PROJ_MATH_H +#define PROJ_MATH_H +#ifdef __cplusplus +extern "C" { +#endif + +#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) + + +double pj_hypot(double x, double y); +double pj_log1p(double x); +double pj_asinh(double x); + +#define hypot pj_hypot +#define log1p pj_log1p +#define asinh pj_asinh + + +#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */ + +#ifdef __cplusplus +} +#endif +#endif /*PROJ_MATH_H */ diff --git a/src/projects.h b/src/projects.h index 5c8d1301..d6d6e23e 100644 --- a/src/projects.h +++ b/src/projects.h @@ -92,11 +92,6 @@ typedef long pj_int32; #define MAX_PATH_FILENAME 1024 #endif -/* prototype hypot for systems where absent */ -#if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) -extern double hypot(double, double); -#endif - /* If we still haven't got M_PI*, we rely on our own defines. * For example, this is necessary when compiling with gcc and * the -ansi flag. -- cgit v1.2.3 From e197e40e9b9ae99b8dac569428af35036822dca2 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 16 Apr 2018 23:19:42 +0200 Subject: Prefer proj_math.h over math.h when geodesic.c is compiled as part of PROJ --- src/geodesic.c | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'src') diff --git a/src/geodesic.c b/src/geodesic.c index eb956c65..5cb273ac 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -24,7 +24,11 @@ */ #include "geodesic.h" +#ifdef PJ_LIB__ +#include "proj_math.h" +#else #include +#endif #if !defined(HAVE_C99_MATH) #define HAVE_C99_MATH 0 -- cgit v1.2.3 From 65d335589b381f9e7b79a1869f74389c3ff59ae5 Mon Sep 17 00:00:00 2001 From: Kurt Schwehr Date: Thu, 19 Apr 2018 09:53:35 -0700 Subject: Rename level -> log_level for proj_log_level to match function definition Found on https://github.com/OSGeo/proj.4/commit/2f082b70cbdafdea49bb123e027406089e7ad65b http://clang.llvm.org/extra/clang-tidy/checks/readability-inconsistent-declaration-parameter-name.html function 'proj_log_level' has a definition with different parameter names src/pj_internal.c:374:14: the definition seen here src/proj.h:364:14: differing parameters are named here: ('level'), in definition: ('log_level') --- src/proj.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/proj.h b/src/proj.h index 3d64ed1a..d254eade 100644 --- a/src/proj.h +++ b/src/proj.h @@ -361,7 +361,7 @@ int proj_errno_reset (const PJ *P); int proj_errno_restore (const PJ *P, int err); const char* proj_errno_string (int err); -PJ_LOG_LEVEL proj_log_level (PJ_CONTEXT *ctx, PJ_LOG_LEVEL level); +PJ_LOG_LEVEL proj_log_level (PJ_CONTEXT *ctx, PJ_LOG_LEVEL log_level); void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf); /* Scaling and angular distortion factors */ -- cgit v1.2.3 From 613d95ef99f29fcf42fab27d1ef26adb87745e35 Mon Sep 17 00:00:00 2001 From: Kurt Schwehr Date: Sun, 22 Apr 2018 07:52:59 -0700 Subject: PJ_isea.c: change local helper functions to void returns The return values are not used and do not mean anything. Similar to #423 --- src/PJ_isea.c | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) (limited to 'src') diff --git a/src/PJ_isea.c b/src/PJ_isea.c index 57412a38..6c870a2b 100644 --- a/src/PJ_isea.c +++ b/src/PJ_isea.c @@ -29,7 +29,7 @@ struct hex { /* y *must* be positive down as the xy /iso conversion assumes this */ ISEA_STATIC -int hex_xy(struct hex *h) { +void hex_xy(struct hex *h) { if (!h->iso) return 1; if (h->x >= 0) { h->y = -h->y - (h->x+1)/2; @@ -38,12 +38,10 @@ int hex_xy(struct hex *h) { h->y = -h->y - h->x/2; } h->iso = 0; - - return 1; } ISEA_STATIC -int hex_iso(struct hex *h) { +void hex_iso(struct hex *h) { if (h->iso) return 1; if (h->x >= 0) { @@ -55,12 +53,10 @@ int hex_iso(struct hex *h) { h->z = -h->x - h->y; h->iso = 1; - return 1; } ISEA_STATIC -int hexbin2(double width, double x, double y, - int *i, int *j) { +void hexbin2(double width, double x, double y, int *i, int *j) { double z, rx, ry, rz; double abs_dx, abs_dy, abs_dz; int ix, iy, iz, s; @@ -105,7 +101,6 @@ int hexbin2(double width, double x, double y, hex_xy(&h); *i = h.x; *j = h.y; - return ix * 100 + iy; } #ifndef ISEA_STATIC #define ISEA_STATIC @@ -1134,4 +1129,3 @@ PJ *PROJECTION(isea) { return P; } - -- cgit v1.2.3 From 454c1326e5f1cae87712b011434776375ac84fbb Mon Sep 17 00:00:00 2001 From: Kurt Schwehr Date: Sun, 22 Apr 2018 08:17:56 -0700 Subject: More cleanup and fix last commit Tried building before committing time --- src/PJ_isea.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/PJ_isea.c b/src/PJ_isea.c index 6c870a2b..3e941c77 100644 --- a/src/PJ_isea.c +++ b/src/PJ_isea.c @@ -30,7 +30,7 @@ struct hex { /* y *must* be positive down as the xy /iso conversion assumes this */ ISEA_STATIC void hex_xy(struct hex *h) { - if (!h->iso) return 1; + if (!h->iso) return; if (h->x >= 0) { h->y = -h->y - (h->x+1)/2; } else { @@ -42,7 +42,7 @@ void hex_xy(struct hex *h) { ISEA_STATIC void hex_iso(struct hex *h) { - if (h->iso) return 1; + if (h->iso) return; if (h->x >= 0) { h->y = (-h->y - (h->x+1)/2); @@ -600,27 +600,27 @@ isea_grid_init(struct isea_dgg * g) } ISEA_STATIC -int +void isea_orient_isea(struct isea_dgg * g) { if (!g) - return 0; + return; g->o_lat = ISEA_STD_LAT; g->o_lon = ISEA_STD_LON; g->o_az = 0.0; - return 1; + return; } ISEA_STATIC -int +void isea_orient_pole(struct isea_dgg * g) { if (!g) - return 0; + return; g->o_lat = M_PI / 2.0; g->o_lon = 0.0; g->o_az = 0; - return 1; + return; } ISEA_STATIC -- cgit v1.2.3 From 6f640d76e13cb643574b44b7498d51ecff6fb83e Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sun, 22 Apr 2018 22:39:57 +0200 Subject: Add isnan() to proj_math.h Code updated to use isnan() instead of pj_is_nan(). --- src/PJ_robin.c | 5 +++-- src/geodesic.c | 3 +-- src/nad_intr.c | 5 +++-- src/pj_apply_vgridshift.c | 4 ++-- src/pj_internal.c | 16 ---------------- src/pj_math.c | 7 ++++++- src/proj_math.h | 2 ++ 7 files changed, 17 insertions(+), 25 deletions(-) (limited to 'src') diff --git a/src/PJ_robin.c b/src/PJ_robin.c index 92aebfd3..a1d5ad4d 100644 --- a/src/PJ_robin.c +++ b/src/PJ_robin.c @@ -1,4 +1,5 @@ #define PJ_LIB__ +#include "proj_math.h" #include "proj_internal.h" #include "proj.h" #include "projects.h" @@ -82,7 +83,7 @@ static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */ (void) P; dphi = fabs(lp.phi); - i = pj_is_nan(lp.phi) ? -1 : (int)floor(dphi * C1); + i = isnan(lp.phi) ? -1 : (int)floor(dphi * C1); if( i < 0 ){ proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); return xy; @@ -117,7 +118,7 @@ static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */ } } else { /* general problem */ /* in Y space, reduce to table interval */ - i = pj_is_nan(lp.phi) ? -1 : (int)floor(lp.phi * NODES); + i = isnan(lp.phi) ? -1 : (int)floor(lp.phi * NODES); if( i < 0 || i >= NODES ) { proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); return lp; diff --git a/src/geodesic.c b/src/geodesic.c index 5cb273ac..9904c7fa 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -243,8 +243,7 @@ static void sincosdx(real x, real* sinx, real* cosx) { r = remquo(x, (real)(90), &q); #else r = fmod(x, (real)(360)); - /* check for NaN -- do not use pj_is_nan, since we want geodesic.c not to - * depend on the rest of proj.4 */ + /* check for NaN */ q = r == r ? (int)(floor(r / 90 + (real)(0.5))) : 0; r -= 90 * q; #endif diff --git a/src/nad_intr.c b/src/nad_intr.c index dc245831..80e428ed 100644 --- a/src/nad_intr.c +++ b/src/nad_intr.c @@ -1,6 +1,7 @@ /* Determine nad table correction value */ #define PJ_LIB__ #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" LP @@ -13,9 +14,9 @@ nad_intr(LP t, struct CTABLE *ct) { int in; t.lam /= ct->del.lam; - indx.lam = pj_is_nan(t.lam) ? 0 : (int)floor(t.lam); + indx.lam = isnan(t.lam) ? 0 : (int)floor(t.lam); t.phi /= ct->del.phi; - indx.phi = pj_is_nan(t.phi) ? 0 : (int)floor(t.phi); + indx.phi = isnan(t.phi) ? 0 : (int)floor(t.phi); frct.lam = t.lam - indx.lam; frct.phi = t.phi - indx.phi; diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c index 05803bc8..e7106b25 100644 --- a/src/pj_apply_vgridshift.c +++ b/src/pj_apply_vgridshift.c @@ -29,7 +29,7 @@ #define PJ_LIB__ #include -#include +#include "proj_math.h" #include "proj_internal.h" #include "projects.h" @@ -42,7 +42,7 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR float *cvs; /* do not deal with NaN coordinates */ /* cppcheck-suppress duplicateExpression */ - if( pj_is_nan(input.phi) || pj_is_nan(input.lam) ) + if( isnan(input.phi) || isnan(input.lam) ) itable = *gridlist_count_p; /* keep trying till we find a table that works */ diff --git a/src/pj_internal.c b/src/pj_internal.c index 891e0b9f..dc528649 100644 --- a/src/pj_internal.c +++ b/src/pj_internal.c @@ -438,19 +438,3 @@ void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf) { if (0!=logf) ctx->logger = logf; } - - -#if HAVE_C99_MATH -/* proj_internal.h defines pj_is_nan as isnan */ -#else -/*****************************************************************************/ -int pj_is_nan (double val) { -/****************************************************************************** - Returns 0 if not a NaN and non-zero if val is a NaN. - - Provides an equivalent to isnan(). -******************************************************************************/ - /* cppcheck-suppress duplicateExpression */ - return val != val; -} -#endif diff --git a/src/pj_math.c b/src/pj_math.c index d5242636..751c2750 100644 --- a/src/pj_math.c +++ b/src/pj_math.c @@ -61,6 +61,11 @@ double pj_asinh(double x) { return x > 0 ? y : (x < 0 ? -y : x); } -#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */ +/* Returns 0 if not a NaN and non-zero if val is a NaN */ +int pj_isnan (double x) { + /* cppcheck-suppress duplicateExpression */ + return x != x; +} +#endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */ diff --git a/src/proj_math.h b/src/proj_math.h index 0d28cb30..0108a439 100644 --- a/src/proj_math.h +++ b/src/proj_math.h @@ -39,10 +39,12 @@ extern "C" { double pj_hypot(double x, double y); double pj_log1p(double x); double pj_asinh(double x); +int pj_isnan(double x); #define hypot pj_hypot #define log1p pj_log1p #define asinh pj_asinh +#define isnan pj_isnan #endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */ -- cgit v1.2.3 From 0dc275a312c1b2ddefc066c66787eb8274d60f08 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 25 Apr 2018 22:03:41 +0200 Subject: pj_transform: reset error state before each call to pj_inv/pj_fwd Fixes issues raised in https://lists.osgeo.org/pipermail/gdal-dev/2018-April/048446.html The use case is that pj_transform() is called from geos projection to something else, with multiple coordinates. If one of the coordinate tuple fails the inverse transform of geos, it fails with PJD_ERR_TOLERANCE_CONDITION. Causing all following coordinates to fail since the error state is not reset. --- src/pj_transform.c | 6 ++++++ 1 file changed, 6 insertions(+) (limited to 'src') diff --git a/src/pj_transform.c b/src/pj_transform.c index 22b685b2..aaa8c1e4 100644 --- a/src/pj_transform.c +++ b/src/pj_transform.c @@ -42,6 +42,8 @@ enum PJ_DIRECTION { }; typedef enum PJ_DIRECTION PJ_DIRECTION; +/* Copied from proj.h FIXME */ +int proj_errno_reset (const PJ *P); static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag, @@ -199,6 +201,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double * if (geodetic_loc.u == HUGE_VAL) continue; + proj_errno_reset( P ); projected_loc = pj_fwd3d( geodetic_loc, P); if( P->ctx->last_errno != 0 ) { @@ -234,6 +237,7 @@ static int geographic_to_projected (PJ *P, long n, int dist, double *x, double * if( geodetic_loc.u == HUGE_VAL ) continue; + proj_errno_reset( P ); projected_loc = pj_fwd( geodetic_loc, P ); if( P->ctx->last_errno != 0 ) { @@ -301,6 +305,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double * if (projected_loc.u == HUGE_VAL) continue; + proj_errno_reset( P ); geodetic_loc = pj_inv3d(projected_loc, P); if( P->ctx->last_errno != 0 ) { @@ -337,6 +342,7 @@ static int projected_to_geographic (PJ *P, long n, int dist, double *x, double * if( projected_loc.u == HUGE_VAL ) continue; + proj_errno_reset( P ); geodetic_loc = pj_inv( projected_loc, P ); if( P->ctx->last_errno != 0 ) { -- cgit v1.2.3 From 56f089510a170164ffc4792a7f7c57ba8f64f1c9 Mon Sep 17 00:00:00 2001 From: Martin Desruisseaux Date: Fri, 27 Apr 2018 10:47:44 +0200 Subject: Remove src/org_proj4_Projections.h file. This file was not used anywhere, since it was replaced by org_proj4_PJ file since 2012. --- src/Makefile.am | 2 +- src/lib_proj.cmake | 3 +-- src/org_proj4_Projections.h | 37 ------------------------------------- 3 files changed, 2 insertions(+), 40 deletions(-) delete mode 100644 src/org_proj4_Projections.h (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index 78ecfdf8..971ac78c 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -10,7 +10,7 @@ AM_CPPFLAGS = -DPROJ_LIB=\"$(pkgdatadir)\" \ -DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ include_HEADERS = proj.h proj_api.h projects.h geodesic.h \ - org_proj4_Projections.h org_proj4_PJ.h + org_proj4_PJ.h EXTRA_DIST = makefile.vc proj.def bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \ bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake \ diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index c0d0a6a9..82a8fc81 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -258,8 +258,7 @@ if(JNI_SUPPORT) set(SRC_LIBPROJ_CORE ${SRC_LIBPROJ_CORE} jniproj.c ) set(HEADERS_LIBPROJ ${HEADERS_LIBPROJ} - org_proj4_PJ.h - org_proj4_Projections.h) + org_proj4_PJ.h) source_group("Source Files\\JNI" FILES ${SRC_LIBPROJ_JNI}) add_definitions(-DJNI_ENABLED) include_directories( ${JNI_INCLUDE_DIRS}) diff --git a/src/org_proj4_Projections.h b/src/org_proj4_Projections.h deleted file mode 100644 index 3841e057..00000000 --- a/src/org_proj4_Projections.h +++ /dev/null @@ -1,37 +0,0 @@ -/* DO NOT EDIT THIS FILE - it is machine generated */ -#include -/* Header for class org_proj4_Projections */ - -#ifndef _Included_org_proj4_Projections -#define _Included_org_proj4_Projections -#ifdef __cplusplus -extern "C" { -#endif -/* - * Class: org_proj4_Projections - * Method: getProjInfo - * Signature: (Ljava/lang/String;)Ljava/lang/String; - */ -JNIEXPORT jstring JNICALL Java_org_proj4_Projections_getProjInfo - (JNIEnv *, jobject, jstring); - -/* - * Class: org_proj4_Projections - * Method: getEllipsInfo - * Signature: (Ljava/lang/String;)Ljava/lang/String; - */ -JNIEXPORT jstring JNICALL Java_org_proj4_Projections_getEllipsInfo - (JNIEnv *, jobject, jstring); - -/* - * Class: org_proj4_Projections - * Method: transform - * Signature: ([D[D[DLjava/lang/String;Ljava/lang/String;JI)V - */ -JNIEXPORT void JNICALL Java_org_proj4_Projections_transform - (JNIEnv *, jobject, jdoubleArray, jdoubleArray, jdoubleArray, jstring, jstring, jlong, jint); - -#ifdef __cplusplus -} -#endif -#endif -- cgit v1.2.3 From 44b83cbb394639f6c393c349b56360af27bc3020 Mon Sep 17 00:00:00 2001 From: Martin Desruisseaux Date: Fri, 27 Apr 2018 15:42:36 +0200 Subject: Apply the https://github.com/opengeospatial/geoapi/issues/26 patch. Update contributor names in copyright header. --- src/jniproj.c | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/jniproj.c b/src/jniproj.c index 85d984da..072fa04c 100644 --- a/src/jniproj.c +++ b/src/jniproj.c @@ -5,7 +5,9 @@ * Martin Desruisseaux * ****************************************************************************** - * Copyright (c) 2005, Antonello Andrea + * Copyright (c) 2005, Andrea Antonello + * Copyright (c) 2011, Martin Desruisseaux + * Copyright (c) 2018, Even Rouault * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), @@ -351,11 +353,8 @@ void convertAngularOrdinates(PJ *pj, double* data, jint numPts, int dimension, d if (pj_is_latlong(pj)) { /* Convert only the 2 first ordinates and skip all the other dimensions. */ dimToSkip = dimension - 2; - } else if (pj_is_geocent(pj)) { - /* Convert only the 3 first ordinates and skip all the other dimensions. */ - dimToSkip = dimension - 3; } else { - /* Not a geographic or geocentric CRS: nothing to convert. */ + /* Not a geographic CRS: nothing to convert. */ return; } double *stop = data + dimension*numPts; -- cgit v1.2.3 From aef6d05bfe86b810dbd38dbdb9c757dae7b8368e Mon Sep 17 00:00:00 2001 From: Martin Desruisseaux Date: Fri, 27 Apr 2018 18:15:08 +0200 Subject: Update the Ant build.xml script for: - Generating the C header file during javac task. - Provide more accurate information in META-INF/MANIFEST.MF. - Merge tasks intended to be executed together. The src/org_proj4_PJ.h file has been re-generated with the most recent Java compiler. --- src/org_proj4_PJ.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/org_proj4_PJ.h b/src/org_proj4_PJ.h index 7e690e7a..be5d3f58 100644 --- a/src/org_proj4_PJ.h +++ b/src/org_proj4_PJ.h @@ -44,7 +44,7 @@ JNIEXPORT jstring JNICALL Java_org_proj4_PJ_getDefinition /* * Class: org_proj4_PJ * Method: getType - * Signature: ()Lorg/proj4/PJ$Type; + * Signature: ()Lorg/proj4/PJ/Type; */ JNIEXPORT jobject JNICALL Java_org_proj4_PJ_getType (JNIEnv *, jobject); -- cgit v1.2.3 From 4bb5771fbc813ce3a71413418f4e926978f21be8 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 30 Apr 2018 09:45:32 +0200 Subject: Remove unneeded pj_is_nan definition --- src/proj_internal.h | 6 ------ 1 file changed, 6 deletions(-) (limited to 'src') diff --git a/src/proj_internal.h b/src/proj_internal.h index 1aca9201..30a4a89e 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -118,12 +118,6 @@ void proj_fileapi_set (PJ *P, void *fileapi); const char * const *proj_get_searchpath(void); int proj_get_path_count(void); -#if HAVE_C99_MATH -#define pj_is_nan isnan -#else -int pj_is_nan (double val); -#endif - #ifdef __cplusplus } #endif -- cgit v1.2.3