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 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