From 44c611c574f001267cd39c904c081e8f1f0d8c96 Mon Sep 17 00:00:00 2001 From: Thomas Knudsen Date: Thu, 15 Dec 2016 16:46:54 +0100 Subject: Horner and helmert (#456) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Introducing the Horner polynomial evaluator also introduces the need for very long +init:tag arguments (a n'th order 2D polynomium has (n+1)(n+2)/2 coefficients, and n is typically in the range 5-10, i.e. up to around 60 coefficients for each polynomium, and there are 4 polynomia in a complete back/forward transformation set). Hence, in this commit, along with the first part of the Horner code, the code for reading +init files has been modified in a (for all practical purposes) backwards compatible way, by making it possible to introduce line continuations by escaping line breaks, i.e. preceding them with a backslash. An escaped line break works (as it would in TeX), by skipping all following whitespace, including interspersed #-comments. This simple extension makes it possible to create very long initialization elements without losing track of the structure (cf. s45b.pol and pj_init_test.c in the examples-directory for a demo). The s45b.pol file was created by hand-editing the output of the software doing the original constrained adjustment for the polynomial coefficients. The simple adding of the “skip following whitespace and comments” feature has made it possible to retain almost all metadata from the source material. This is considered very important, since 1) For the lack of a prior common file format for geodetic polynomial coefficients, there is a good chance that this will become THE standard, at least for the time being, and 2) Without the metadata represented, it will be very hard for a human to debug code involving a slightly misrepresented polynomium. Due to the current architecture of the pj_init.c code (mostly around the fill_buffer() function), it is next to impossible to implement the line continuation functionality in full generality. Hence, it has been necessary to limit this format extension to files smaller than 64 kB. * Correction of spherical HEALpix test case The first HEALpix test case in nad/testvarious is clearly intended to invoke the spherical form of HEALpix. It does, however, specify the spheroid using the +a=1 size parameter, without specifying any shape parameter. But since +no_defs is not specified either, a shape parameter is picked up from the nad/proj_def.dat file (where ellps=WGS84 is given in the section). It appears that this has not happened before I updated the pj_init code to support projection pipelines (see below). I do, however, believe that the present behaviour is the correct one, and rather than retrohacking the pj_init code, to (incorrectly, I believe) reproduce the prior behaviour, I have corrected the test case invocation in nad/testvarious to specify the spheroid using the +R=1 size parameter (which was already used in the following test case). * Repair scaling of projections stomping on value of semimajor axis * Workaround MSVC HUGE_VAL misimplementation. The "return const err object" idiom (i.e. const err = {HUGE_VAL,...}; ... if (bad) return err) is problematic to implement due to MSVC's misimplementation of HUGE_VAL as a non-const. Hence, we need to run-time initialize these. In the pj_inv functions, this was mistakenly done to the wrong object. For pj_fwdobs/invobs and the remaining part of the obs-based API, this is now worked around by providing functions returning a run time HUGE_VAL initialized PJ_OBS or PJ_COO resp. Obnoxious, but given MSVC's market penetration there is really not much else we can do. --- examples/pj_obs_api_test.c | 33 ++-- examples/s45b.pol | 402 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 418 insertions(+), 17 deletions(-) create mode 100644 examples/s45b.pol (limited to 'examples') diff --git a/examples/pj_obs_api_test.c b/examples/pj_obs_api_test.c index c75b4a21..bb792c62 100644 --- a/examples/pj_obs_api_test.c +++ b/examples/pj_obs_api_test.c @@ -44,8 +44,6 @@ int pj_pipeline_selftest (void); - - int main (void) { PJ *P; PJ_OBS a, b; @@ -71,7 +69,7 @@ int main (void) { return puts ("Oops"), 0; /* zero initialize everything, then set (longitude, latitude) to (12, 55) */ - a = pj_obs_null; + a = pj_obs_null (); /* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */ a.coo.lp.lam = TORAD(12); a.coo.lp.phi = TORAD(55); @@ -100,13 +98,12 @@ int main (void) { dist = pj_xyz_dist (a.coo.xyz, b.coo.xyz); printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9); - /* Invalid projection */ a = pj_trans (P, 42, a); if (a.coo.lpz.lam!=DBL_MAX) printf ("%15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi); - err = pj_error (P); - printf ("pj_error: %d\n", err); + err = pj_err_level (P, PJ_ERR_TELL); + printf ("pj_err_level: %d\n", err); /* Clean up */ pj_free (P); @@ -117,7 +114,7 @@ int main (void) { return puts ("Oops"), 0; /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */ - a = b = pj_obs_null; + a = b = pj_obs_null (); a.coo.lpz.lam = TORAD(12); a.coo.lpz.phi = TORAD(55); a.coo.lpz.z = 100; @@ -152,17 +149,18 @@ int main (void) { P I P E L I N E T E S T S - **************************************************************************** - - ***************************************************************************/ /* forward-reverse geo->utm->geo */ - P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +step +proj=utm +inv"); + P = pj_create ( + "+proj=pipeline +ellps=GRS80 +zone=32 +step " + "+proj=utm +step " + "+proj=utm +inv" + ); if (0==P) return puts ("Oops"), 0; /* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */ - a = b = pj_obs_null; + a = b = pj_obs_null (); a.coo.lpz.lam = TORAD(12); a.coo.lpz.phi = TORAD(55); a.coo.lpz.z = 100; @@ -180,12 +178,15 @@ int main (void) { /* And now the back-to-back situation utm->geo->utm */ - P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +inv +step +proj=utm"); + P = pj_create ( + "+proj=pipeline +ellps=GRS80 +zone=32 +step " + "+proj=utm +inv +step " + "+proj=utm"); if (0==P) return puts ("Oops"), 0; /* zero initialize everything, then set (easting, northing) to utm(12, 55) */ - a = b = pj_obs_null; + a = b = pj_obs_null (); a.coo.xy = cph_utm32; printf ("PRE: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y); @@ -207,7 +208,7 @@ int main (void) { return puts ("Oops"), 0; /* zero initialize everything, then set (easting, northing) to utm(12, 55) */ - a = b = pj_obs_null; + a = b = pj_obs_null (); a.coo.lpz.lam = TORAD(12); a.coo.lpz.phi = TORAD(55); printf ("PRE: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi)); @@ -221,7 +222,6 @@ int main (void) { a = pj_trans (P, PJ_INV, b); printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi)); - /* Geodesic distance between two points with angular 2D coordinates */ a.coo.lp.lam = TORAD(12); a.coo.lp.phi = TORAD(60); @@ -237,7 +237,6 @@ int main (void) { dist = pj_lp_dist (P, a.coo.lp, b.coo.lp); printf ("1 deg at Equator: %15.9f\n", dist); - pj_free (P); pj_pipeline_selftest (); return 0; diff --git a/examples/s45b.pol b/examples/s45b.pol new file mode 100644 index 00000000..9b70fbf7 --- /dev/null +++ b/examples/s45b.pol @@ -0,0 +1,402 @@ +################################################################################ +# +# Horner polynomial evaluation demo data +# +################################################################################ +# +# These are the polynomial coefficients used for transforming to and from +# the Danish legacy system s45b "System 45 Bornholm" +# +################################################################################ + + +proj=pipeline + +step init=./s45b.pol:s45b_tc32 +# step init=./s45b.pol:tc32_utm32 + +step proj=utm inv ellps=intl zone=32 +step proj=cart ellps=intl + +step proj=helmert ellps=GRS80 + x=-81.0703 y=-89.3603 z=-115.7526 + rx=-484.88 ry=-24.36 rz=-413.21 s=-0.540645 + +step proj=cart inv ellps=GRS80 +step proj=utm ellps=GRS80 zone=32 + + + + +################################################################################ +# +# S 4 5 B -> T C 3 2 +# +################################################################################ + + +proj=horner +ellps=intl +range=500000 + +fwd_origin=47022.563745,51779.260103 +inv_origin=878354.943082,6125305.175366 +deg=6 + + +# static double C_ttb[] +# tc32_ed50 -> s45b +# m_lim_gen: 0.153 red = 0 OBS = 1074 +# m = 1.51 cm my_loss = +3 y_enp = +8.4 +# m = 1.53 cm mx_loss = +4 x_enp = +8.4 + +# mht C_ttb er +# fwd-inv ombyttet ifht original Poder/Engsager-kode +# For at opnå at to fwd transform fører fra s45b->tc32->utm32 (->ETRS89) + +inv_v=\ +# Poly NORTH :: e-degree = 0 : n-degree = 6 + 5.1779004699e+04,\ + 9.9508320295e-01,\ + -2.9453823207e-10,\ + 1.9995084102e-14,\ + -1.4895751366e-18,\ + -9.9734812211e-23,\ + 1.1194218845e-26,\ + +# Poly NORTH :: e-degree = 1 : n-degree = 5 + -8.4285679515e-02,\ + -7.9623049286e-09,\ + -3.7190046062e-14,\ + -2.3324127411e-18,\ + -1.1150449763e-22,\ + 2.8703154270e-27,\ + +# Poly NORTH :: e-degree = 2 : n-degree = 4 + 8.7160434140e-10,\ + -3.3634602927e-14,\ + -5.5718245313e-18,\ + 6.2611750909e-23,\ + -2.1011243838e-26,\ + +# Poly NORTH :: e-degree = 3 : n-degree = 3 + 1.0905463989e-14,\ + -4.3960034360e-18,\ + 3.6121595001e-22,\ + -1.3493066011e-27,\ + +# Poly NORTH :: e-degree = 4 : n-degree = 2 + -1.3360171462e-18,\ + 1.0780850646e-22,\ + 4.5118286607e-26,\ + +# Poly NORTH :: e-degree = 5 : n-degree = 1 + -1.3718883973e-22,\ + 1.6263920750e-26,\ + +# Poly NORTH :: e-degree = 6 : n-degree = 0 + -5.1004217526e-27 + +# tcy 6125305.175366 + +inv_u=\ +# Poly EAST :: n-degree = 0 : e-degree = 6 + 4.7022495967e+04,\ + -9.9508282498e-01,\ + 3.2436283039e-09,\ + -2.6276394334e-15,\ + 8.6318533291e-18,\ + -3.8327518550e-23,\ + -2.5704924282e-26,\ + +# Poly EAST :: n-degree = 1 : e-degree = 5 + -8.4285975934e-02,\ + 5.7098765263e-10,\ + -6.0863955939e-14,\ + 2.3608788740e-18,\ + 6.8899581969e-24,\ + -1.1429511179e-26,\ + +# Poly EAST :: n-degree = 2 : e-degree = 4 + -4.6079778412e-09,\ + 1.5072604543e-14,\ + 5.4063862378e-18,\ + 1.2591327827e-22,\ + 7.9336388691e-27,\ + +# Poly EAST :: n-degree = 3 : e-degree = 3 + -2.9479268638e-14,\ + 1.7090049434e-18,\ + 2.8413337985e-22,\ + -3.3577391552e-27,\ + +# Poly EAST :: n-degree = 4 : e-degree = 2 + 3.0434879273e-18,\ + -1.8081673510e-22,\ + -2.3651419850e-26,\ + +# Poly EAST :: n-degree = 5 : e-degree = 1 + 9.2060044804e-23,\ + 3.7807953325e-27,\ + +# Poly EAST :: n-degree = 6 : e-degree = 0 + -4.9415665221e-27 + +# tcx 878354.943082 + + +# static double C_btt[] +# s45b -> tc32_ed50 +# m_lim_gen: 0.154 red = 0 OBS = 1074 +# m = 1.50 cm my_loss = +3 y_enp = +8.5 +# m = 1.54 cm mx_loss = +4 x_enp = +8.3 + +fwd_v=\ +# Poly NORTH :: e-degree = 0 : n-degree = 6 + 6.1253054245e+06,\ + 9.9778251908e-01,\ + -7.7346152025e-10,\ + -2.5359789369e-14,\ + 1.5614918228e-18,\ + 9.8091134295e-23,\ + -1.1092581145e-26,\ + +# Poly NORTH :: e-degree = 1 : n-degree = 5 + -8.4514352088e-02,\ + -7.9847579284e-09,\ + -2.6865560962e-14,\ + -2.0731372756e-18,\ + -1.3660341123e-22,\ + 1.1244836340e-26,\ + +# Poly NORTH :: e-degree = 2 : n-degree = 4 + 8.0551988135e-11,\ + 3.6661500679e-14,\ + 5.4247705403e-18,\ + 8.4494604807e-23,\ + 1.3334858516e-26,\ + +# Poly NORTH :: e-degree = 3 : n-degree = 3 + 8.3889821184e-15,\ + -4.8124202237e-18,\ + 2.9088188830e-22,\ + -2.0129874264e-26,\ + +# Poly NORTH :: e-degree = 4 : n-degree = 2 + 2.4716463766e-18,\ + -2.1717177513e-22,\ + -3.2828537638e-26,\ + +# Poly NORTH :: e-degree = 5 : n-degree = 1 + -1.2080655753e-22,\ + 2.5050435391e-26,\ + +# Poly NORTH :: e-degree = 6 : n-degree = 0 + 1.1383483826e-27 + +# tcy 51779.260103, + +fwd_u=\ +# Poly EAST :: n-degree = 0 : e-degree = 6 + 8.7835485387e+05,\ + -9.9778289691e-01,\ + 3.2537215213e-09,\ + 6.9217640616e-15,\ + 8.6268883840e-18,\ + 4.6748156909e-23,\ + -2.6492402009e-26,\ + +# Poly EAST :: n-degree = 1 : e-degree = 5 + -8.4514648771e-02,\ + 1.4399520180e-09,\ + -6.0423329711e-14,\ + 6.9816167332e-20,\ + 6.7729233542e-23,\ + -5.3308251880e-27,\ + +# Poly EAST :: n-degree = 2 : e-degree = 4 + -4.5697800099e-09,\ + -1.5194038814e-14,\ + 5.1112653016e-18,\ + -2.0307532869e-22,\ + 1.0374125432e-26,\ + +# Poly EAST :: n-degree = 3 : e-degree = 3 + -2.8983003841e-14,\ + -1.6414425785e-18,\ + 1.7874983379e-22,\ + 1.5492164174e-26,\ + +# Poly EAST :: n-degree = 4 : e-degree = 2 + 2.7919197366e-18,\ + 1.9218613279e-22,\ + -2.1007264634e-26,\ + +# Poly EAST :: n-degree = 5 : e-degree = 1 + 1.0032412389e-22,\ + -5.9007997846e-27,\ + +# Poly EAST :: n-degree = 6 : e-degree = 0 + -4.4410970979e-27 + +# tcx 47022.563745 + + + + + +################################################################################ +# +# T C 3 2 -> U T M 3 2 +# +################################################################################ + + + +proj=horner +ellps=intl +range=500000 + + +fwd_origin=877605.269066,6125810.306769 +inv_origin=877605.760036,6125811.281773 + + +# tc32_ed50 -> utm32_ed50 : Bornholm + +deg=4 + +# ttu_n and ttu_e are based on static double C_ttu_b[] +# m_lim_gen: 0.086 red = 0 OBS = 852 +# m = 1.38 cm my_loss = +2 y_enp = +10.5 +# m = 1.44 cm mx_loss = +2 x_enp = +10.4 +# static double ttu_n[] + +fwd_v=\ +# Poly NORTH :: e-degree = 0 : n-degree = 0..4 + 6.1258112678e+06,\ + 9.9999971567e-01,\ + 1.5372750011e-10,\ + 5.9300860915e-15,\ + 2.2609497633e-19,\ + +# Poly NORTH :: e-degree = 1 : n-degree = 0..3 + 4.3188227445e-05,\ + 2.8225130416e-10,\ + 7.8740007114e-16,\ + -1.7453997279e-19,\ + +# Poly NORTH :: e-degree = 2 : n-degree = 0..2 + 1.6877465415e-10,\ + -1.1234649773e-14,\ + -1.7042333358e-18,\ + +# Poly NORTH :: e-degree = 3 : n-degree = 0..1 + -7.9303467953e-15,\ + -5.2906832535e-19,\ + +# Poly NORTH :: e-degree = 4 : n-degree = 0 + 3.9984284847e-19 + +# tcy 6125810.306769 + + + +# static double ttu_e[] +fwd_u=\ +# Poly EAST :: n-degree = 0 : e-degree = 0..4 + 8.7760574982e+05,\ + 9.9999752475e-01,\ + 2.8817299305e-10,\ + 5.5641310680e-15,\ + -1.5544700949e-18,\ + +# Poly EAST :: n-degree = 1 : e-degree = 3 + -4.1357045890e-05,\ + 4.2106213519e-11,\ + 2.8525551629e-14,\ + -1.9107771273e-18,\ + +# Poly EAST :: n-degree = 2 : e-degree = 2 + 3.3615590093e-10,\ + 2.4380247154e-14,\ + -2.0241230315e-18,\ + +# Poly EAST :: n-degree = 3 : e-degree = 1 + 1.2429019719e-15,\ + 5.3886155968e-19,\ + +# Poly EAST :: n-degree = 4 : e-degree = 0 + -1.0167505000e-18 + +# tcx 877605.760036 + + +# utt_n and utt_e are based on static double C_utt_b[] +# utm32_ed50 -> tc32_ed50 : Bornholm +# m_lim_gen: 0.086 red = 0 OBS = 852 +# m = 1.38 cm my_loss = +2 y_enp = +10.8 +# m = 1.44 cm mx_loss = +2 x_enp = +10.7 +# static double utt_n[] + +inv_v=\ +# Poly NORTH :: e-degree = 0 : n-degree = 4 + 6.1258103208e+06,\ + 1.0000002826e+00,\ + -1.5372762184e-10,\ + -5.9304261011e-15,\ + -2.2612705361e-19,\ + +# Poly NORTH :: e-degree = 1 : n-degree = 3 + -4.3188331419e-05,\ + -2.8225549995e-10,\ + -7.8529116371e-16,\ + 1.7476576773e-19,\ + +# Poly NORTH :: e-degree = 2 : n-degree = 2 + -1.6875687989e-10,\ + 1.1236475299e-14,\ + 1.7042518057e-18,\ + +# Poly NORTH :: e-degree = 3 : n-degree = 1 + 7.9300735257e-15,\ + 5.2881862699e-19,\ + +# Poly NORTH :: e-degree = 4 : n-degree = 0 + -3.9990736798e-19 + +# tcy 6125811.281773 + + + + +# static double utt_e[] +inv_u=\ +# Poly EAST :: n-degree = 0 : e-degree = 0..4 + 8.7760527928e+05,\ + 1.0000024735e+00,\ + -2.8817540032e-10,\ + -5.5627059451e-15,\ + 1.5543637570e-18,\ + +# Poly EAST :: n-degree = 1 : e-degree = 0..3 + 4.1357152105e-05,\ + -4.2114813612e-11,\ + -2.8523713454e-14,\ + 1.9109017837e-18,\ + +# Poly EAST :: n-degree = 2 : e-degree = 0..2 + -3.3616407783e-10,\ + -2.4382678126e-14,\ + 2.0245020199e-18,\ + +# Poly EAST :: n-degree = 3 : e-degree = 0..1 + -1.2441377565e-15,\ + -5.3885232238e-19,\ + +# Poly EAST :: n-degree = 4 : e-degree = 0 + 1.0167203661e-18 + +# tcx 877605.760036 + +<> -- cgit v1.2.3