aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas Knudsen <busstoptaktik@users.noreply.github.com>2016-12-15 16:46:54 +0100
committerGitHub <noreply@github.com>2016-12-15 16:46:54 +0100
commit44c611c574f001267cd39c904c081e8f1f0d8c96 (patch)
tree708136e458bfa83e51f4ea4fc287bcbcb295ae74
parentc64461525d3f08ca93fcf950b76a665cfe142082 (diff)
downloadPROJ-44c611c574f001267cd39c904c081e8f1f0d8c96.tar.gz
PROJ-44c611c574f001267cd39c904c081e8f1f0d8c96.zip
Horner and helmert (#456)
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 <general> 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 <type> 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.
Diffstat
-rw-r--r--examples/pj_obs_api_test.c33
-rw-r--r--examples/s45b.pol402
-rwxr-xr-xnad/testvarious4
-rw-r--r--src/Makefile.am2
-rw-r--r--src/PJ_calcofi.c4
-rw-r--r--src/PJ_cart.c110
-rw-r--r--src/PJ_helmert.c492
-rw-r--r--src/PJ_horner.c953
-rw-r--r--src/PJ_krovak.c2
-rw-r--r--src/PJ_pipeline.c54
-rw-r--r--src/lib_proj.cmake2
-rw-r--r--src/makefile.vc3
-rw-r--r--src/pj_fwd.c85
-rw-r--r--src/pj_fwd3d.c96
-rw-r--r--src/pj_init.c76
-rw-r--r--src/pj_inv.c60
-rw-r--r--src/pj_inv3d.c92
-rw-r--r--src/pj_list.h2
-rw-r--r--src/pj_obs_api.c81
-rw-r--r--src/pj_param.c12
-rw-r--r--src/pj_strerrno.c1
-rw-r--r--src/proj.def41
-rw-r--r--src/proj.h19
23 files changed, 2307 insertions, 319 deletions
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"
+#
+################################################################################
+
+<s45b>
+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
+#
+################################################################################
+<s45b_tc32>
+
+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
+#
+################################################################################
+
+<tc32_utm32>
+
+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
+
+<>
diff --git a/nad/testvarious b/nad/testvarious
index 58770347..f11bc806 100755
--- a/nad/testvarious
+++ b/nad/testvarious
@@ -283,8 +283,8 @@ $EXE +proj=hammer +datum=WGS84 \
EOF
echo "##############################################################" >> ${OUT}
echo "Test healpix forward projection on sphere" >> ${OUT}
-$EXE +proj=latlong +a=1 +lon_0=0 \
- +to +proj=healpix +a=1 +lon_0=0 -f '%.'5'f' \
+$EXE +proj=latlong +R=1 +lon_0=0 \
+ +to +proj=healpix +R=1 +lon_0=0 -f '%.'5'f' \
-E >>${OUT} <<EOF
0 41.81031
-90 0
diff --git a/src/Makefile.am b/src/Makefile.am
index 69e0e44e..a2508fe8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -74,7 +74,7 @@ libproj_la_SOURCES = \
jniproj.c pj_mutex.c pj_initcache.c pj_apply_vgridshift.c geodesic.c \
pj_strtod.c \
\
- pj_obs_api.c PJ_cart.c PJ_pipeline.c
+ pj_obs_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c
install-exec-local:
rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT)
diff --git a/src/PJ_calcofi.c b/src/PJ_calcofi.c
index b9a8ca22..0f206155 100644
--- a/src/PJ_calcofi.c
+++ b/src/PJ_calcofi.c
@@ -174,8 +174,8 @@ int pj_calcofi_selftest (void) {
double tolerance_lp = 1e-10;
double tolerance_xy = 1e-7;
- char e_args[] = {"+proj=calcofi +ellps=GRS80 +lat_1=0.5 +lat_2=2"};
- char s_args[] = {"+proj=calcofi +a=6400000 +lat_1=0.5 +lat_2=2"};
+ char e_args[] = {"+proj=calcofi +ellps=GRS80 +lat_1=0.5 +lat_2=2 +no_defs"};
+ char s_args[] = {"+proj=calcofi +R=6400000 +lat_1=0.5 +lat_2=2 +no_defs"};
LP fwd_in[] = {
{ 2, 1},
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 4ec20ced..c585fbeb 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -88,6 +88,16 @@ PROJ_HEAD(cart, "Geodetic/cartesian conversions");
(TF, below).
+ Close to the poles, we avoid singularities by switching to an
+ approximation requiring knowledge of the geocentric radius
+ at the given latitude. For this, we use an adaptation of the
+ formula given in:
+
+ Wikipedia: Earth Radius
+ https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude
+ (Derivation and commentary at http://gis.stackexchange.com/questions/20200/how-do-you-compute-the-earths-radius-at-a-given-geodetic-latitude)
+
+ (WP2, below)
These routines are probably not as robust at those in
geocent.c, at least thay haven't been through as heavy
@@ -103,13 +113,6 @@ static void freeup (PJ *P) {
return;
}
-
-/*********************************************************************/
-static double semiminor_axis (double a, double es) {
-/*********************************************************************/
- return a * sqrt (1 - es);
-}
-
/*********************************************************************/
static double normal_radius_of_curvature (double a, double es, double phi) {
/*********************************************************************/
@@ -120,15 +123,11 @@ static double normal_radius_of_curvature (double a, double es, double phi) {
return a / sqrt (1 - es*s*s);
}
-
/*********************************************************************/
-static double second_eccentricity_squared (double a, double es) {
-/**********************************************************************
- Follows the definition, e'2 = (a2-b2)/b2, but uses the identity
- (a2-b2) = (a-b)(a+b), for improved numerical precision.
-***********************************************************************/
- double b = semiminor_axis (a, es);
- return (a - b) * (a + b) / (b*b);
+static double geocentric_radius (double a, double b, double phi) {
+/*********************************************************************/
+ /* This is from WP2, but uses hypot() for potentially better numerical robustness */
+ return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi));
}
@@ -148,21 +147,6 @@ static XYZ cartesian (LPZ geodetic, PJ *P) {
xyz.y = (N + h) * cosphi * sin(lam);
xyz.z = (N * (1 - P->es) + h) * sin(phi);
- /*********************************************************************/
- /* */
- /* For historical reasons, in proj, plane coordinates are measured */
- /* in units of the semimajor axis. Since 3D handling is grafted on */
- /* later, this is not the case for heights. And even though this */
- /* coordinate really is 3D cartesian, the z-part looks like a height */
- /* to proj. Hence, we have the somewhat unusual situation of having */
- /* a point coordinate with differing units between dimensions. */
- /* */
- /* The scaling and descaling is handled by the pj_fwd/inv functions. */
- /* */
- /*********************************************************************/
- xyz.x /= P->a;
- xyz.y /= P->a;
-
return xyz;
}
@@ -173,15 +157,12 @@ static LPZ geodetic (XYZ cartesian, PJ *P) {
double N, b, p, theta, c, s, e2s;
LPZ lpz;
- cartesian.x *= P->a;
- cartesian.y *= P->a;
-
/* Perpendicular distance from point to Z-axis (HM eq. 5-28) */
p = hypot (cartesian.x, cartesian.y);
/* Ancillary ellipsoidal parameters */
- b = semiminor_axis (P->a, P->es);
- e2s = second_eccentricity_squared (P->a, P->es);
+ b = P->b;
+ e2s = P->e2s;
/* HM eq. (5-37) */
theta = atan2 (cartesian.z * P->a, p * b);
@@ -196,11 +177,12 @@ static LPZ geodetic (XYZ cartesian, PJ *P) {
c = cos(lpz.phi);
if (fabs(c) < 1e-6) {
- /* poleward of 89.99994 deg, we avoid division by zero */
- /* by computing the height as the cartesian z value */
- /* minus the semiminor axis length */
- /* (potential improvement: approx. by 2nd order poly) */
- lpz.z = cartesian.z - (cartesian.z > 0? b: -b);
+ /* poleward of 89.99994 deg, we avoid division by zero */
+ /* by computing the height as the cartesian z value */
+ /* minus the geocentric radius of the Earth at the given */
+ /* latitude */
+ double r = geocentric_radius (P->a, b, lpz.phi);
+ lpz.z = fabs (cartesian.z) - r;
}
else
lpz.z = p / c - N;
@@ -239,6 +221,8 @@ PJ *PROJECTION(cart) {
P->inv3d = geodetic;
P->fwd = cart_forward;
P->inv = cart_reverse;
+ P->left = PJ_IO_UNITS_RADIANS;
+ P->right = PJ_IO_UNITS_METERS;
return P;
}
@@ -255,9 +239,6 @@ int pj_cart_selftest (void) {
char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
- /* Log everything libproj offers to log for you */
- pj_log_level (0, PJ_LOG_TRACE);
-
/* An utm projection on the GRS80 ellipsoid */
P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==P)
@@ -269,7 +250,7 @@ int pj_cart_selftest (void) {
/* Same projection, now using argc/argv style initialization */
P = pj_create_argv (3, args);
if (0==P)
- return puts ("Oops"), 0;
+ return 2;
/* Turn off logging */
pj_log_level (0, PJ_LOG_NONE);
@@ -294,18 +275,18 @@ int pj_cart_selftest (void) {
dist = pj_xy_dist (a.coo.xy, b.coo.xy);
if (dist > 2e-9)
- return 2;
+ return 3;
/* Invalid projection */
a = pj_trans (P, 42, a);
- if (a.coo.lpz.lam!=DBL_MAX)
- return 3;
- err = pj_error (P);
- if (0==err)
+ if (a.coo.lpz.lam!=HUGE_VAL)
return 4;
+ err = pj_err_level (P, PJ_ERR_TELL);
+ if (0==err)
+ return 5;
/* Clear error */
- pj_error_set (P, 0);
+ pj_err_level (P, 0);
/* Clean up */
pj_free (P);
@@ -313,7 +294,7 @@ int pj_cart_selftest (void) {
/* Now do some 3D transformations */
P = pj_create ("+proj=cart +ellps=GRS80");
if (0==P)
- return 5;
+ return 6;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
a = b = pj_obs_null;
@@ -328,9 +309,32 @@ int pj_cart_selftest (void) {
dist = pj_roundtrip (P, PJ_FWD, 10000, a);
dist = pj_roundtrip (P, PJ_INV, 10000, b);
if (dist > 2e-9)
- return 6;
+ return 7;
+
+
+ /* Test at the North Pole */
+ a = b = pj_obs_null;
+ a.coo.lpz.lam = TORAD(0);
+ a.coo.lpz.phi = TORAD(90);
+ a.coo.lpz.z = 100;
+
+ /* Forward projection: Ellipsoidal-to-3D-Cartesian */
+ dist = pj_roundtrip (P, PJ_FWD, 1, a);
+ if (dist > 1e-12)
+ return 8;
+
+ /* Test at the South Pole */
+ a = b = pj_obs_null;
+ a.coo.lpz.lam = TORAD(0);
+ a.coo.lpz.phi = TORAD(-90);
+ a.coo.lpz.z = 100;
+
+ /* Forward projection: Ellipsoidal-to-3D-Cartesian */
+ dist = pj_roundtrip (P, PJ_FWD, 1, a);
+ if (dist > 1e-12)
+ return 9;
- /* Inverse projection: Ellipsoidal-to-3D-Cartesian */
+ /* Inverse projection: 3D-Cartesian-to-Ellipsoidal */
b = pj_trans (P, PJ_INV, b);
/* Move p to another context */
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
new file mode 100644
index 00000000..72cefe90
--- /dev/null
+++ b/src/PJ_helmert.c
@@ -0,0 +1,492 @@
+/***********************************************************************
+
+ Helmert transformations, 3- and 7-parameter
+
+ Thomas Knudsen, 2016-05-24
+
+************************************************************************
+
+ Implements 3- and 7-parameter Helmert transformations for 3D data.
+
+ Primarily useful for implementation of datum shifts in transformation
+ pipelines.
+
+************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-24/06-05
+Last update: 2016-11-24
+
+************************************************************************
+* Copyright (c) 2016, Thomas Knudsen / SDFE
+*
+* Permission is hereby granted, free of charge, to any person obtaining a
+* copy of this software and associated documentation files (the "Software"),
+* to deal in the Software without restriction, including without limitation
+* the rights to use, copy, modify, merge, publish, distribute, sublicense,
+* and/or sell copies of the Software, and to permit persons to whom the
+* Software is furnished to do so, subject to the following conditions:
+*
+* The above copyright notice and this permission notice shall be included
+* in all copies or substantial portions of the Software.
+*
+* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+* DEALINGS IN THE SOFTWARE.
+*
+***********************************************************************/
+
+#define PJ_LIB__
+#include <proj.h>
+#include <projects.h>
+#include <geocent.h>
+#include <assert.h>
+#include <stddef.h>
+#include <errno.h>
+PROJ_HEAD(helmert, "3- and 7-parameter Helmert shift");
+
+static XYZ helmert_forward_3d (LPZ lpz, PJ *P);
+static LPZ helmert_reverse_3d (XYZ xyz, PJ *P);
+
+
+
+static void *freeup_msg (PJ *P, int errlev) { /* Destructor */
+ if (0==P)
+ return 0;
+
+ if (0!=P->ctx)
+ pj_ctx_set_errno (P->ctx, errlev);
+
+ if (0==P->opaque)
+ return pj_dealloc (P);
+
+ pj_dealloc (P->opaque);
+ return pj_dealloc(P);
+}
+
+
+/* Adapts pipeline_freeup to the format defined for the PJ object */
+static void freeup (PJ *P) {
+ freeup_msg (P, 0);
+ return;
+}
+
+
+/***********************************************************************/
+struct pj_opaque_helmert {
+/************************************************************************
+ Projection specific elements for the "helmert" PJ object
+************************************************************************/
+ XYZ xyz;
+ PJ_OPK opk;
+ double scale;
+ double R[3][3];
+ int no_rotation;
+};
+
+
+/* Make the maths of the rotation operations somewhat more readable and textbook like */
+#define R00 (Q->R[0][0])
+#define R01 (Q->R[0][1])
+#define R02 (Q->R[0][2])
+
+#define R10 (Q->R[1][0])
+#define R11 (Q->R[1][1])
+#define R12 (Q->R[1][2])
+
+#define R20 (Q->R[2][0])
+#define R21 (Q->R[2][1])
+#define R22 (Q->R[2][2])
+
+
+
+
+
+
+/***********************************************************************/
+static XY helmert_forward (LP lp, PJ *P) {
+/************************************************************************
+
+ The 2D implementation is a slightly silly stop-gap, returning a
+ 2D sub-space of the 3D transformation.
+
+ In order to follow the "principle of least astonishment", the 2D
+ interface should provide an implementation of the 4 parameter
+ 2D Helmert shift.
+
+************************************************************************/
+ PJ_TRIPLET point;
+ point.lp = lp;
+ point.lpz.z = 0;
+ point.xyz = helmert_forward_3d (point.lpz, P);
+ return point.xy;
+}
+
+
+/***********************************************************************/
+static LP helmert_reverse (XY xy, PJ *P) {
+/***********************************************************************/
+ PJ_TRIPLET point;
+ point.xy = xy;
+ point.xyz.z = 0;
+ point.lpz = helmert_reverse_3d (point.xyz, P);
+ return point.lp;
+}
+
+
+/***********************************************************************/
+static XYZ helmert_forward_3d (LPZ lpz, PJ *P) {
+/***********************************************************************/
+ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
+ PJ_TRIPLET point;
+ double X, Y, Z, scale;
+ if (Q->no_rotation) {
+ point.xyz.x = lpz.lam + Q->xyz.x;
+ point.xyz.y = lpz.phi + Q->xyz.y;
+ point.xyz.z = lpz.z + Q->xyz.z;
+ return point.xyz;
+ }
+
+ scale = 1 + Q->scale * 1e-6;
+
+ X = lpz.lam;
+ Y = lpz.phi;
+ Z = lpz.z;
+
+ point.lpz = lpz;
+
+ point.xyz.x = scale * ( R00 * X + R01 * Y + R02 * Z);
+ point.xyz.y = scale * ( R10 * X + R11 * Y + R12 * Z);
+ point.xyz.z = scale * ( R20 * X + R21 * Y + R22 * Z);
+
+ point.xyz.x += Q->xyz.x;
+ point.xyz.y += Q->xyz.y;
+ point.xyz.z += Q->xyz.z;
+
+ return point.xyz;
+}
+
+
+/***********************************************************************/
+static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) {
+/***********************************************************************/
+ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
+ PJ_TRIPLET point;
+ double X, Y, Z, scale;
+
+ point.xyz = xyz;
+
+ if (Q->no_rotation) {
+ point.xyz.x = xyz.x - Q->xyz.x;
+ point.xyz.y = xyz.y - Q->xyz.y;
+ point.xyz.z = xyz.z - Q->xyz.z;
+ return point.lpz;
+ }
+
+ scale = 1 + Q->scale * 1e-6;
+
+ /* Unscale and deoffset */
+ X = (xyz.x - Q->xyz.x) / scale;
+ Y = (xyz.y - Q->xyz.y) / scale;
+ Z = (xyz.z - Q->xyz.z) / scale;
+
+ /* Inverse rotation through transpose multiplication */
+ point.xyz.x = ( R00 * X + R10 * Y + R20 * Z);
+ point.xyz.y = ( R01 * X + R11 * Y + R21 * Z);
+ point.xyz.z = ( R02 * X + R12 * Y + R22 * Z);
+
+ return point.lpz;
+}
+
+
+static PJ_OBS helmert_forward_obs (PJ_OBS point, PJ *P) {
+ point.coo.xyz = helmert_forward_3d (point.coo.lpz, P);
+ return point;
+}
+
+
+static PJ_OBS helmert_reverse_obs (PJ_OBS point, PJ *P) {
+ point.coo.lpz = helmert_reverse_3d (point.coo.xyz, P);
+ return point;
+}
+
+
+/* Milliarcsecond to radians */
+#define MAS_TO_RAD (DEG_TO_RAD / 3600000.0)
+
+
+/***********************************************************************/
+PJ *PROJECTION(helmert) {
+/***********************************************************************/
+ XYZ xyz = {0, 0, 0};
+ PJ_OPK opk = {0, 0, 0};
+ double scale = 0;
+ int approximate = 0, transpose = 0;
+
+ struct pj_opaque_helmert *Q = pj_calloc (1, sizeof (struct pj_opaque_helmert));
+ if (0==Q)
+ return freeup_msg (P, ENOMEM);
+ P->opaque = (void *) Q;
+
+ P->fwdobs = helmert_forward_obs;
+ P->invobs = helmert_reverse_obs;
+ P->fwdobs = 0;
+ P->invobs = 0;
+ P->fwd3d = helmert_forward_3d;
+ P->inv3d = helmert_reverse_3d;
+ P->fwd = helmert_forward;
+ P->inv = helmert_reverse;
+
+ P->left = PJ_IO_UNITS_METERS;
+ P->right = PJ_IO_UNITS_METERS;
+
+ /* Translations */
+ if (pj_param (P->ctx, P->params, "tx").i)
+ xyz.x = pj_param (P->ctx, P->params, "dx").f;
+
+ if (pj_param (P->ctx, P->params, "ty").i)
+ xyz.y = pj_param (P->ctx, P->params, "dy").f;
+
+ if (pj_param (P->ctx, P->params, "tz").i)
+ xyz.z = pj_param (P->ctx, P->params, "dz").f;
+
+ /* Rotations */
+ if (pj_param (P->ctx, P->params, "trx").i)
+ opk.o = pj_param (P->ctx, P->params, "drx").f * MAS_TO_RAD;
+
+ if (pj_param (P->ctx, P->params, "try").i)
+ opk.p = pj_param (P->ctx, P->params, "dry").f * MAS_TO_RAD;
+
+ if (pj_param (P->ctx, P->params, "trz").i)
+ opk.k = pj_param (P->ctx, P->params, "drz").f * MAS_TO_RAD;
+
+
+ /* Scale */
+ if (pj_param (P->ctx, P->params, "ts").i)
+ scale = pj_param (P->ctx, P->params, "ds").f;
+
+ /* Use small angle approximations? */
+ if (pj_param (P->ctx, P->params, "bapprox").i)
+ approximate = 1;
+
+ /* Use "other" rotation sign convention? */
+ if (pj_param (P->ctx, P->params, "ttranspose").i)
+ transpose = 1;
+
+ Q->xyz = xyz;
+ Q->opk = opk;
+ Q->scale = scale;
+
+ if ((opk.o==0) && (opk.p==0) && (opk.k==0) && (scale==0)) {
+ Q->no_rotation = 1;
+ return P;
+ }
+
+
+ /**************************************************************************
+
+ Build rotation matrix.
+ ----------------------
+
+ Here we rename rotation indices from omega, phi, kappa (opk), to
+ fi (i.e. phi), theta, psi (ftp), in order to reduce the mental agility
+ needed to implement the expression for the rotation matrix derived over
+ at https://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions
+
+ If the option "approximate" is set, small angle approximations are used:
+ The matrix elements are approximated by expanding the trigonometric
+ functions to linear order (i.e. cos(x) = 1, sin(x) = x), and discarding
+ products of second order.
+
+ This was a useful hack when calculating by hand was the only option,
+ but in general, today, should be avoided because:
+
+ 1. It does not save much computation time, as the rotation matrix
+ is built only once, and probably used many times.
+
+ 2. The error induced may be too large for ultra high accuracy
+ applications: the Earth is huge and the linear error is
+ approximately the angular error multiplied by the Earth radius.
+
+ However, in many cases the approximation is necessary, since it has
+ been used historically: Rotation angles from older published datum
+ shifts may actually be a least squares fit to the linearized rotation
+ approximation, hence not being strictly valid for deriving the full
+ rotation matrix.
+
+ So in order to fit historically derived coordinates, the access to
+ the approximate rotation matrix is necessary - at least in principle.
+
+ Also, when using any published datum transformation information, one
+ should always check which convention (exact or approximate rotation
+ matrix) is expected, and whether the induced error for selecting
+ the opposite convention is acceptable (which it often is).
+
+
+ Sign conventions
+ ----------------
+
+ Take care: Two different sign conventions exist.
+
+ Conceptually they relate to whether we rotate the coordinate system
+ or the "position vector" (the vector going from the coordinate system
+ origin to the point being transformed, i.e. the point coordinates
+ interpreted as vector coordinates).
+
+ Switching between the "position vector" and "coordinate system"
+ conventions is simply a matter of switching the sign of the rotation
+ angles, which algebraically also translates into a transposition of
+ the rotation matrix.
+
+ Hence, as geodetic constants should preferably be referred to exactly
+ as published, the "transpose" option provides the ability to switch
+ between the conventions.
+
+ ***************************************************************************/
+ do {
+ double f, t, p; /* phi/fi , theta, psi */
+ double cf, ct, cp; /* cos (fi, theta, psi) */
+ double sf, st, sp; /* sin (fi, theta, psi) */
+
+ /* rename (omega, phi, kappa) to (fi, theta, psi) */
+ f = opk.o;
+ t = opk.p;
+ p = opk.k;
+
+ if (approximate) {
+ R00 = 1;
+ R01 = p;
+ R02 = -t;
+
+ R10 = -p;
+ R11 = 1;
+ R12 = f;
+
+ R20 = t;
+ R21 = -f;
+ R22 = 1;
+ break;
+ }
+
+ cf = cos(f);
+ sf = sin(f);
+ ct = cos(t);
+ st = sin(t);
+ cp = cos(p);
+ sp = sin(p);
+
+
+ R00 = ct*cp;
+ R01 = cf*sp + sf*st*cp;
+ R02 = sf*sp - cf*st*cp;
+
+ R10 = -ct*sp;
+ R11 = cf*cp - sf*st*sp;
+ R12 = sf*cp + cf*st*sp;
+
+ R20 = st;
+ R21 = -sf*ct;
+ R22 = cf*ct;
+
+ /*
+ For comparison: Description from Engsager/Poder implementation
+ in set_dtm_1.c (trlib)
+
+ DATUM SHIFT:
+ TO = scale * ROTZ * ROTY * ROTX * FROM + TRANSLA
+
+ ( cz sz 0) (cy 0 -sy) (1 0 0)
+ ROTZ=(-sz cz 0), ROTY=(0 1 0), ROTX=(0 cx sx)
+ ( 0 0 1) (sy 0 cy) (0 -sx cx)
+
+ trp->r11 = cos_ry*cos_rz;
+ trp->r12 = cos_rx*sin_rz + sin_rx*sin_ry*cos_rz;
+ trp->r13 = sin_rx*sin_rz - cos_rx*sin_ry*cos_rz;
+
+ trp->r21 = -cos_ry*sin_rz;
+ trp->r22 = cos_rx*cos_rz - sin_rx*sin_ry*sin_rz;
+ trp->r23 = sin_rx*cos_rz + cos_rx*sin_ry*sin_rz;
+
+ trp->r31 = sin_ry;
+ trp->r32 = -sin_rx*cos_ry;
+ trp->r33 = cos_rx*cos_ry;
+
+ trp->scale = 1.0 + scale;
+ */
+
+ } while (0);
+
+
+ if (transpose) {
+ double r;
+ r = R01; R01 = R10; R10 = r;
+ r = R02; R02 = R20; R20 = r;
+ r = R12; R12 = R21; R21 = r;
+ }
+
+ return P;
+}
+
+
+
+
+#ifndef PJ_SELFTEST
+
+int pj_helmert_selftest (void) {return 0;} /* self test stub */
+
+#else
+
+
+static int test (char *args, PJ_TRIPLET in, PJ_TRIPLET expect) {
+ PJ_TRIPLET out;
+ PJ *P = pj_init_plus (args);
+
+ if (0==P)
+ return 5;
+
+ out.xyz = pj_fwd3d (in.lpz, P);
+ if (pj_xyz_dist (out.xyz, expect.xyz) > 1e-4)
+ return 1;
+
+ out.lpz = pj_inv3d (out.xyz, P);
+ if (pj_xyz_dist (out.xyz, in.xyz) > 1e-4)
+ return 2;
+
+ return 0;
+}
+
+
+
+int pj_helmert_selftest (void) {
+
+ /* This example is from
+ Lotti Jivall:
+ Simplified transformations from ITRF2008/IGS08 to ETRS89 for maritime applications */
+ PJ_TRIPLET in1 = {{3565285.0000, 855949.0000, 5201383.0000}};
+ PJ_TRIPLET expect1 = {{3565285.41342351, 855948.67986759, 5201382.72939791}};
+ char args1[] = {
+ " +proj=helmert +ellps=GRS80"
+ " +x=0.67678 +y=0.65495 +z=-0.52827"
+ " +rx=-22.742 +ry=12.667 +rz=22.704 +s=-0.01070"
+ };
+
+
+ /* This example is a random point, transformed from ED50 to ETRS89 using KMStrans2 */
+ PJ_TRIPLET in2 = {{3494994.3012, 1056601.9725, 5212382.1666}};
+ PJ_TRIPLET expect2 = {{3494909.84026368, 1056506.78938633, 5212265.66699761}};
+ char args2[] = {
+ " +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"
+ };
+ int ret;
+
+ ret = test (args1, in1, expect1); if (ret) return ret;
+ ret = test (args2, in2, expect2); if (ret) return ret + 10;
+ return 0;
+}
+
+#endif
diff --git a/src/PJ_horner.c b/src/PJ_horner.c
new file mode 100644
index 00000000..f39d2c2c
--- /dev/null
+++ b/src/PJ_horner.c
@@ -0,0 +1,953 @@
+#define PJ_LIB__
+#include <proj.h>
+#include <projects.h>
+#include <assert.h>
+#include <stddef.h>
+#include <math.h>
+#include <errno.h>
+PROJ_HEAD(horner, "Horner polynomial evaluation");
+#define HORNER_SILENCE
+
+/* The next few hundred lines comprises a direct cut-and-paste from the horner.h header library */
+
+/***********************************************************************
+
+ Interfacing to a classic piece of geodetic software
+
+************************************************************************
+
+ gen_pol is a highly efficient, classic implementation of a generic
+ 2D Horner's Scheme polynomial evaluation routine by Knud Poder and
+ Karsten Engsager, originating in the vivid geodetic environment at
+ what was then (1960-ish) the Danish Geodetic Institute.
+
+ The original Poder/Engsager gen_pol implementation (where
+ the polynomial degree and two sets of polynomial coefficients
+ are packed together in one compound array, handled via a simple
+ double pointer) is compelling and "true to the code history":
+
+ It has a beautiful classical 1960s ring to it, not unlike the
+ original fft implementations, which revolutionized spectral
+ analysis in twenty lines of code.
+
+ The Poder coding sound, as classic 1960s as Phil Spector's Wall
+ of Sound, is beautiful and inimitable.
+
+ On the other hand: For the uninitiated, the gen_pol code is hard
+ to follow, despite being compact.
+
+ Also, since adding metadata and improving maintainability
+ of the code are among the implied goals of a current SDFE/DTU Space
+ project, the material in this file introduces a version with a more
+ more modern (or at least 1990s) look, introducing a "double 2D
+ polynomial" data type, HORNER.
+
+ Despite introducing a new data type for handling the polynomial
+ coefficients, great care has been taken to keep the coefficient
+ array organization identical to that of gen_pol.
+
+ Hence, on one hand, the HORNER data type helps improving the
+ long term maintainability of the code by making the data
+ organization more mentally accessible.
+
+ On the other hand, it allows us to preserve the business end of
+ the original gen_pol implementation - although not including the
+ famous "Poder dual autocheck" in all its enigmatic elegance.
+
+ The original code has, however, been included in the conditionally
+ compiled TEST-section.
+
+ This is partially for validation of the revised version, partially
+ to enable more generic enjoyment of an interesting piece of
+ ingenious geodetic code - simplistic and enigmatic at the same time.
+
+ **********************************************************************
+
+ The material included here was written by Knud Poder, starting
+ around 1960, and Karsten Engsager, starting around 1970. It was
+ originally written in Algol 60, later (1980s) reimplemented in C.
+
+ The HORNER data type interface, and the organization as a header
+ library was implemented by Thomas Knudsen, starting around 2015.
+
+ ***********************************************************************
+ *
+ * The gen_pol routine comes with this legal statement (ISC/OpenBSD):
+ *
+ * Copyright (c) 2011, National Survey and Cadastre, Denmark
+ * (Kort- og Matrikelstyrelsen), kms@kms.dk
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ *
+ *
+ *******************************************************************************
+ *
+ * The remaining parts are...
+ *
+ * Copyright (c) 2016, Thomas Knudsen / Karsten Engsager / SDFE http://www.sdfe.dk
+ *
+ * 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.
+ *
+ *****************************************************************************/
+
+
+ #ifndef HORNER_H
+ #define HORNER_H
+ #ifdef __cplusplus
+ extern "C" {
+ #endif
+
+
+
+#if defined(PROJ_H) || defined(PROJECTS_H)
+#define horner_dealloc(x) pj_dealloc(x)
+#define horner_calloc(n,x) pj_calloc(n,x)
+#else
+#define horner_dealloc(x) free(x)
+#define horner_calloc(n,x) calloc(n,x)
+typedef struct {double u,v;} UV;
+#endif
+
+struct horner;
+typedef struct horner HORNER;
+static UV horner (const HORNER *transformation, int direction, UV position);
+static HORNER *horner_alloc (size_t order);
+static void horner_free (HORNER *h);
+
+struct horner {
+ int order; /* maximum degree of polynomium */
+ int coefs; /* number of coefficients for each polynomium */
+ double range; /* radius of the region of validity */
+
+ double *fwd_u; /* coefficients for the forward transformations */
+ double *fwd_v; /* i.e. latitude/longitude to northing/easting */
+
+ double *inv_u; /* coefficients for the inverse transformations */
+ double *inv_v; /* i.e. northing/easting to latitude/longitude */
+
+ UV *fwd_origin; /* False longitude/latitude */
+ UV *inv_origin; /* False easting/northing */
+};
+
+/* e.g. degree = 2: a + bx + cy + dxx + eyy + fxy, i.e. 6 coefficients */
+#define horner_number_of_coefficients(order) \
+ (((order + 1)*(order + 2)/2))
+
+static int horner_degree_u (int order, int index);
+static int horner_degree_v (int order, int index);
+static int horner_index (int order, int degree_u, int degree_v);
+
+#ifndef HORNER_HEADER_ONLY
+
+/***************************************************************************/
+static int horner_index (int order, int degree_1, int degree_2) {
+/****************************************************************************
+
+ Returns the index of the polynomial coefficient, C, for the element
+
+ C * pow (c_1, degree_2) * pow (c_2, degree_2),
+
+ given that degree_1 > -1, degree_2 > -1, degree_1 + degree_2 <= order.
+
+ Otherwise returns -1 and sets errno to EDOM.
+
+ The range of the index is [0 : (order + 1) * (order + 2) / 2 - 1].
+
+ A very important thing to note is that the order of the coordinates
+ c_1 and c_2 depend on the polynomium:
+
+ For the fwd and inv polynomia for the "u" coordinate,
+ u is first (degree_1), v is second (degree_2).
+ For the fwd and inv polynomia for the "v" coordinate,
+ v is first (degree_1), u is second (degree_2).
+
+****************************************************************************/
+ if ( (degree_1 < 0) || (degree_2 < 0) || (degree_1 + degree_2 > order) ) {
+ errno = EDOM;
+ return -1;
+ }
+
+ return ( horner_number_of_coefficients(order) - 1
+ - (order - degree_1)*(order - degree_1 + 1)/2
+ - (order - degree_1 - degree_2));
+}
+
+#define index_u(h, u, v) horner_index (h->order, u, v)
+#define index_v(h, u, v) horner_index (h->order, v, u)
+
+
+static int horner_degree_u (int order, int index) {
+ int n = horner_number_of_coefficients(order);
+ int i, j;
+ if ((order < 0) || (index >= n)) {
+ errno = EDOM;
+ return -1;
+ }
+ for (i = 0; i <= order; i++)
+ for (j = 0; j <= order - i; j++)
+ if (index == horner_index (order, i, j))
+ return i;
+ return -1;
+}
+
+
+static int horner_degree_v (int order, int index) {
+ int n = horner_number_of_coefficients(order);
+ int i, j;
+ if ((order < 0) || (index >= n)) {
+ errno = EDOM;
+ return -1;
+ }
+ for (i = 0; i <= order; i++)
+ for (j = 0; j <= order - i; j++)
+ if (index == horner_index (order, i, j))
+ return j;
+ return -1;
+}
+
+
+
+
+static void horner_free (HORNER *h) {
+ horner_dealloc (h->inv_v);
+ horner_dealloc (h->inv_u);
+ horner_dealloc (h->fwd_v);
+ horner_dealloc (h->fwd_u);
+ horner_dealloc (h);
+}
+
+
+static HORNER *horner_alloc (size_t order) {
+ /* size_t is unsigned, so we need not check for order > 0 */
+ int n = horner_number_of_coefficients(order);
+
+ HORNER *h = horner_calloc (1, sizeof (HORNER));
+
+ if (0==h)
+ return 0;
+
+ h->order = order;
+ h->coefs = n;
+
+ h->fwd_u = horner_calloc (n, sizeof(double));
+ h->fwd_v = horner_calloc (n, sizeof(double));
+ h->inv_u = horner_calloc (n, sizeof(double));
+ h->inv_v = horner_calloc (n, sizeof(double));
+
+ h->fwd_origin = horner_calloc (1, sizeof(UV));
+ h->inv_origin = horner_calloc (1, sizeof(UV));
+
+ if (h->fwd_u && h->fwd_v && h->inv_u && h->inv_v && h->fwd_origin && h->inv_origin)
+ return h;
+
+ /* safe, since all pointers are null-initialized (by calloc) */
+ horner_free (h);
+ return 0;
+}
+
+
+
+
+/**********************************************************************/
+static UV horner (const HORNER *transformation, int direction, UV position) {
+/***********************************************************************
+
+A reimplementation of the classic Engsager/Poder 2D Horner polynomial
+evaluation engine "gen_pol".
+
+This version omits the inimitable Poder "dual autocheck"-machinery,
+which here is intended to be implemented at a higher level of the
+library: We separate the polynomial evaluation from the quality
+control (which, given the limited MTBF for "computing machinery",
+typical when Knud Poder invented the dual autocheck method,
+was not defensible at that time).
+
+Another difference from the original version is that we return the
+result on the stack, rather than accepting pointers to result variables
+as input. This results in code that is easy to read:
+
+ projected = horner (s34j, 1, geographic);
+ geographic = horner (s34j, -1, projected );
+
+and experiments have shown that on contemporary architectures, the time
+taken for returning even comparatively large objects on the stack (and
+the UV is not that large - typically only 16 bytes) is negligibly
+different from passing two pointers (i.e. typically also 16 bytes) the
+other way.
+
+The polynomium has the form:
+
+P = sum (i = [0 : order])
+ sum (j = [0 : order - i])
+ pow(par_1, i) * pow(par_2, j) * coef(index(order, i, j))
+
+For numerical stability, the summation is carried out backwards,
+summing the tiny high order elements first.
+
+***********************************************************************/
+
+ /* These variable names follow the Engsager/Poder implementation */
+ int sz; /* Number of coefficients per polynomial */
+ double *tcx, *tcy; /* Coefficient pointers */
+ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */
+ double n, e;
+ UV uv_error;
+ uv_error.u = uv_error.v = HUGE_VAL;
+
+ if (0==transformation)
+ return uv_error;
+
+ /* Check for valid value of direction (-1, 0, 1) */
+ switch (direction) {
+ case 0: /* no-op */
+ return position;
+ case 1: /* forward */
+ case -1: /* inverse */
+ break;
+ default: /* invalid */
+ errno = EINVAL;
+ return uv_error;
+ }
+
+ /* Prepare for double Horner */
+ sz = horner_number_of_coefficients(transformation->order);
+ range = transformation->range;
+
+
+ if (direction==1) { /* forward */
+ tcx = transformation->fwd_u + sz;
+ tcy = transformation->fwd_v + sz;
+ e = position.u - transformation->fwd_origin->u;
+ n = position.v - transformation->fwd_origin->v;
+ } else { /* inverse */
+ tcx = transformation->inv_u + sz;
+ tcy = transformation->inv_v + sz;
+ e = position.u - transformation->inv_origin->u;
+ n = position.v - transformation->inv_origin->v;
+ }
+
+ if ((fabs(n) > range) || (fabs(e) > range)) {
+ errno = EDOM;
+ return uv_error;
+ }
+
+ /* The melody of this block is straight out of the great Engsager/Poder songbook */
+ else {
+ int g = transformation->order;
+ int r = g, c;
+ double u, v, N, E;
+
+ /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */
+ for (N = *--tcy, E = *--tcx; r > 0; r--) {
+ for (c = g, u = *--tcy, v = *--tcx; c >= r; c--) {
+ u = n*u + *--tcy;
+ v = e*v + *--tcx;
+ }
+ N = e*N + u;
+ E = n*E + v;
+ }
+
+ position.u = E;
+ position.v = N;
+ }
+
+ return position;
+}
+
+
+#ifdef HORNER_SILENCE
+/**********************************************************************/
+static int horner_silence (int i) {
+/***********************************************************************
+ useless function that silences coompiler warnings about unused stuff
+***********************************************************************/
+ HORNER *Q;
+ UV uv_error;
+ if (i==0)
+ return i;
+ uv_error.u = uv_error.v = HUGE_VAL;
+ horner(0, 1, uv_error);
+ Q = horner_alloc (2);
+ if (Q)
+ horner_free (Q);
+ if (horner_degree_u (2,1))
+ return horner_degree_v (2,1);
+ return i;
+}
+#endif /* def HORNER_SILENCE */
+
+
+
+
+#ifdef HORNER_TEST_ORIGINAL_GEN_POL_CODE
+
+double fwd_u[] = {1,2,3};
+double fwd_v[] = {4,5,6};
+double inv_u[] = {4,6,5};
+double inv_v[] = {1,3,2};
+UV uv_origin_fwd = {0, 0};
+UV uv_origin_inv = {0, 0};
+HORNER uv = {1, 3, 500000.0, fwd_u, fwd_v, inv_u, inv_v, &uv_origin_fwd, &uv_origin_inv};
+
+void show_test (void);
+void tuut_b_test (void);
+
+
+int main (int argc, char **argv) {
+ int i, j, order = 10;
+ UV res = {1,1};
+
+ if (argc==2)
+ order = atoi (argv[1]);
+
+ printf ("Testing %d combinations\n", horner_number_of_coefficients(order));
+
+ for (i = 0; i <= order; i++)
+ for (j = 0; j <= order - i; j++) {
+ int hh = horner_index (order, i, j);
+ int ii = horner_degree_u (order, hh);
+ int jj = horner_degree_v (order, hh);
+ assert (i==ii);
+ assert (j==jj);
+ printf ("%2.2d %1d%1du %1d%1dv\n", hh, i, ii, j, jj);
+ }
+
+ tuut_b_test ();
+
+ puts ("Forward...");
+
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, 1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v);
+ assert ( 6==res.u); /* fwd_u: a + bu + cv = 1 + 2*1 + 3*1 = 6 */
+ assert (15==res.v); /* fwd_v: a + bu + cv = 4 + 5*1 + 6*1 = 15 */
+
+ res.u = 2;
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, 1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v); /* u = 2, v = 15 */
+ assert (50==res.u); /* fwd_u: a + bu + cv = 1 + 2*(u=2) + 3*(v=15) = 50 */
+ assert (91==res.v); /* fwd_v: a + bu + cv = 4 + 5*(v=15) + 6*(u=2) = 91 */
+
+ puts ("Inverse...");
+
+ res.u = 1;
+ res.v = 1;
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, -1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v); /* u = 1, v = 1 */
+ assert (15==res.u); /* inv_u: a + bu + cv = 4 + 5*1 + 6*1 = 15 */
+ assert ( 6==res.v); /* inv_v: a + bu + cv = 1 + 2*1 + 3*1 = 6 */
+
+ res.u = 2;
+ printf ("inp = {%.4f, %.4f}\n", res.u, res.v);
+ res = horner (&uv, -1, res);
+ printf ("res = {%.4f, %.4f}\n", res.u, res.v); /* u = 2, v = 6 */
+ assert (46==res.u); /* inv_u: a + bu + cv = 4 + 6*(u=2) + 5*(v=6) = 46 */
+ assert (23==res.v); /* inv_v: a + bu + cv = 1 + 2*(u=2) + 3*(v=6) = 23 */
+
+/* show_test ();*/
+ return 0;
+}
+
+
+
+
+
+
+/* ttu_n and ttu_e are based on static double C_ttu_b[] */
+static double ttu_n[] = {
+ /* tc32_ed50 -> utm32_ed50 : Bornholm */
+ /* 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 */
+
+ /*deg 4.0,*/
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258112678e+06, 9.9999971567e-01, 1.5372750011e-10,
+ /* 3*/ 5.9300860915e-15, 2.2609497633e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ 4.3188227445e-05, 2.8225130416e-10, 7.8740007114e-16,
+ /* 8*/ -1.7453997279e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ 1.6877465415e-10, -1.1234649773e-14, -1.7042333358e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ -7.9303467953e-15, -5.2906832535e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ 3.9984284847e-19,
+
+ /*tcy 6125810.306769, */
+};
+
+static double ttu_e[] = {
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760574982e+05, 9.9999752475e-01, 2.8817299305e-10,
+ /* 3*/ 5.5641310680e-15, -1.5544700949e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ -4.1357045890e-05, 4.2106213519e-11, 2.8525551629e-14,
+ /* 8*/ -1.9107771273e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ 3.3615590093e-10, 2.4380247154e-14, -2.0241230315e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ 1.2429019719e-15, 5.3886155968e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ -1.0167505000e-18,
+
+ /* tcx 877605.269066 */
+ };
+
+
+/* utt_n and utt_e are based on static double C_utt_b[] */
+static double utt_n[] = {
+ /* 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 */
+
+ /*deg 4.0,*/
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258103208e+06, 1.0000002826e+00, -1.5372762184e-10,
+ /* 3*/ -5.9304261011e-15, -2.2612705361e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ -4.3188331419e-05, -2.8225549995e-10, -7.8529116371e-16,
+ /* 8*/ 1.7476576773e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ -1.6875687989e-10, 1.1236475299e-14, 1.7042518057e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ 7.9300735257e-15, 5.2881862699e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ -3.9990736798e-19,
+
+ /*tcy 6125811.281773,*/
+};
+
+static double utt_e[] = {
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760527928e+05, 1.0000024735e+00, -2.8817540032e-10,
+ /* 3*/ -5.5627059451e-15, 1.5543637570e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ 4.1357152105e-05, -4.2114813612e-11, -2.8523713454e-14,
+ /* 8*/ 1.9109017837e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ -3.3616407783e-10, -2.4382678126e-14, 2.0245020199e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ -1.2441377565e-15, -5.3885232238e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ 1.0167203661e-18,
+
+ /*tcx 877605.760036 */
+};
+
+UV tuut_b_origin_fwd = {877605.269066, 6125810.306769};
+UV tuut_b_origin_inv = {877605.760036, 6125811.281773};
+HORNER tuut_b = {4, 15, 500000.0, ttu_e, ttu_n, utt_e, utt_n, &tuut_b_origin_fwd, &tuut_b_origin_inv};
+
+
+/* Prototype and forward declarations of the material needed for cross-check with original implementation */
+int gen_pol(double *C_f, double *C_r, double N_in, double E_in, double *Nout, double *Eout);
+static double C_ttu_b[];
+static double C_utt_b[];
+
+
+void gen_pol_roundtrip (double *C_ttu_b, double *C_utt_b, UV fwd) {
+ UV res, hrn;
+ int ret;
+
+ ret = gen_pol(C_ttu_b, C_utt_b, fwd.v, fwd.u, &res.v, &res.u);
+ printf ("\n------\n\n");
+ if (0!=ret) printf ("ret: %d\n", ret);
+ printf ("inp: %11.3f %11.3f\n", fwd.u, fwd.v);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ hrn = horner (&tuut_b, 1, fwd);
+ printf ("hrn: %11.3f %11.3f\n", hrn.u, hrn.v);
+ assert (hrn.u==res.u);
+ assert (hrn.v==res.v);
+
+ ret = gen_pol(C_utt_b, C_ttu_b, res.v, res.u, &res.v, &res.u);
+ hrn = horner (&tuut_b, -1, hrn);
+ printf ("hrn: %11.3f %11.3f\n", hrn.u, hrn.v);
+ if (0!=ret) printf ("ret: %d\n", ret);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ assert (hrn.u==res.u);
+ assert (hrn.v==res.v);
+ printf ("inp: %11.3f %11.3f (%.3g mm)\n", fwd.u, fwd.v, 1e3*hypot(fwd.u-res.u, fwd.v-res.v));
+}
+
+void tuut_b_test (void) {
+ UV fwd = tuut_b_origin_fwd;
+ UV res;
+
+ puts ("Bornholm");
+ printf ("fwd: %11.3f %11.3f\n", fwd.u, fwd.v);
+ res = horner (&tuut_b, 1, fwd);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ res = horner (&tuut_b, -1, res);
+ printf ("res: %11.3f %11.3f\n", res.u, res.v);
+ printf ("fwd: %11.3f %11.3f\n", fwd.u, fwd.v);
+
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+ fwd.u = 877000;
+ fwd.v = 6125000;
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+ fwd.u = 800000;
+ fwd.v = 6100000;
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+ fwd.u = 850000;
+ fwd.v = 6200000;
+ gen_pol_roundtrip (C_ttu_b, C_utt_b, fwd);
+
+}
+
+
+
+
+
+#define FLOATLIMIT 5.0e5
+#define TRF_AREA_ EDOM
+
+/********************************************************************************************/
+int gen_pol(double *C_f, double *C_r, double N_in, double E_in, double *Nout, double *Eout) {
+/*********************************************************************************************
+
+ This is the original Poder/Engsager implementation of gen_pol.
+ It is included here for test-comparison with the horner() routine.
+
+*********************************************************************************************/
+ double N, E, n, e;
+ double *Cp, *tcy, *tcx;
+ double tol = 1.0e-4;
+ int i;
+ int g;
+ int sz;
+ int r, c;
+ int res = 0;
+
+ /* Preserve input for reverse check */
+ N = N_in;
+ E = E_in;
+ Cp = C_f;
+
+ /* Transformation loop */
+ for (i = -1; i <= 1 && res == 0; ++i)
+ if (i) {
+
+ /* Prepare for double Horner */
+ g = (int) *Cp;
+ sz = (g + 1)*(g + 2)/2 + 1;
+ tcy = Cp + sz;
+ tcx = tcy + sz;
+ /* Double Horner's scheme */
+ /* N = n*Cy*e -> yout, E = e*Cx*n -> xout */
+ n = N - *tcy;
+ e = E - *tcx;
+ if ((fabs(n) < FLOATLIMIT) && (fabs(e) < FLOATLIMIT)) {
+
+ for ( r = g, N = *--tcy, E = *--tcx; r > 0; r--) {
+ double u, v;
+ for (c = g, u = *--tcy, v = *--tcx; c >= r; c--) {
+ u = n*u + *--tcy;
+ v = e*v + *--tcx;
+ }
+ N = e*N + u;
+ E = n*E + v;
+ }
+ } else res = TRF_AREA_;
+ }
+ else { /* collect output coord,switch to reverse checking */
+ *Nout = N;
+ *Eout = E;
+ Cp = C_r;
+ }
+
+ /* tol-check of results*/
+ if (res == 0 && (fabs(N - N_in) < tol && fabs(E - E_in) < tol))
+ return (0);
+ else if (res == 0) res = TRF_AREA_;
+ return(res);
+
+#undef FLOATLIMIT
+
+}
+
+
+/* s45b polynomia */
+ static double C_ttu_b[] = {
+ /* tc32_ed50 -> utm32_ed50 : Bornholm */
+ /* 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 */
+
+ /*deg*/ 4.0,
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258112678e+06, 9.9999971567e-01, 1.5372750011e-10,
+ /* 3*/ 5.9300860915e-15, 2.2609497633e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ 4.3188227445e-05, 2.8225130416e-10, 7.8740007114e-16,
+ /* 8*/ -1.7453997279e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ 1.6877465415e-10, -1.1234649773e-14, -1.7042333358e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ -7.9303467953e-15, -5.2906832535e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ 3.9984284847e-19,
+
+ /*tcy*/ 6125810.306769,
+
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760574982e+05, 9.9999752475e-01, 2.8817299305e-10,
+ /* 3*/ 5.5641310680e-15, -1.5544700949e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ -4.1357045890e-05, 4.2106213519e-11, 2.8525551629e-14,
+ /* 8*/ -1.9107771273e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ 3.3615590093e-10, 2.4380247154e-14, -2.0241230315e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ 1.2429019719e-15, 5.3886155968e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ -1.0167505000e-18,
+
+ /*tcx*/ 877605.269066
+ };
+
+ 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 */
+
+ /*deg*/ 4.0,
+ /*Poly NORTH :: e-degree = 0 : n-degree = 4 */
+ /* 0*/ 6.1258103208e+06, 1.0000002826e+00, -1.5372762184e-10,
+ /* 3*/ -5.9304261011e-15, -2.2612705361e-19,
+
+ /*Poly NORTH :: e-degree = 1 : n-degree = 3 */
+ /* 5*/ -4.3188331419e-05, -2.8225549995e-10, -7.8529116371e-16,
+ /* 8*/ 1.7476576773e-19,
+
+ /*Poly NORTH :: e-degree = 2 : n-degree = 2 */
+ /* 9*/ -1.6875687989e-10, 1.1236475299e-14, 1.7042518057e-18,
+ /*Poly NORTH :: e-degree = 3 : n-degree = 1 */
+ /* 12*/ 7.9300735257e-15, 5.2881862699e-19,
+ /*Poly NORTH :: e-degree = 4 : n-degree = 0 */
+ /* 14*/ -3.9990736798e-19,
+
+ /*tcy*/ 6125811.281773,
+
+ /*Poly EAST :: n-degree = 0 : e-degree = 4 */
+ /* 0*/ 8.7760527928e+05, 1.0000024735e+00, -2.8817540032e-10,
+ /* 3*/ -5.5627059451e-15, 1.5543637570e-18,
+
+ /*Poly EAST :: n-degree = 1 : e-degree = 3 */
+ /* 5*/ 4.1357152105e-05, -4.2114813612e-11, -2.8523713454e-14,
+ /* 8*/ 1.9109017837e-18,
+
+ /*Poly EAST :: n-degree = 2 : e-degree = 2 */
+ /* 9*/ -3.3616407783e-10, -2.4382678126e-14, 2.0245020199e-18,
+ /*Poly EAST :: n-degree = 3 : e-degree = 1 */
+ /* 12*/ -1.2441377565e-15, -5.3885232238e-19,
+ /*Poly EAST :: n-degree = 4 : e-degree = 0 */
+ /* 14*/ 1.0167203661e-18,
+
+ /*tcx*/ 877605.760036
+ };
+
+#endif /* HORNER_TEST_ORIGINAL_GEN_POL_CODE */
+
+
+#endif /* ndef HORNER_HEADER_ONLY */
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* ndef HORNER_H */
+
+
+
+
+
+
+
+static PJ_OBS horner_forward_obs (PJ_OBS point, PJ *P) {
+ point.coo.uv = horner ((HORNER *) P->opaque, 1, point.coo.uv);
+ return point;
+}
+
+static PJ_OBS horner_reverse_obs (PJ_OBS point, PJ *P) {
+ point.coo.uv = horner ((HORNER *) P->opaque, -1, point.coo.uv);
+ return point;
+}
+
+
+
+
+
+
+static void *horner_freeup (PJ *P) { /* Destructor */
+ if (0==P)
+ return 0;
+ if (0==P->opaque)
+ return pj_dealloc (P);
+ horner_free ((HORNER *) P->opaque);
+ return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+ horner_freeup (P);
+ return;
+}
+
+
+static int parse_coefs (PJ *P, double *coefs, char *param, int ncoefs) {
+ char buf[20], *init, *next;
+ int i;
+ sprintf (buf, "t%s", param);
+ if (0==pj_param (P->ctx, P->params, buf).i)
+ return 0;
+ sprintf (buf, "s%s", param);
+ init = pj_param(P->ctx, P->params, buf).s;
+
+ for (i = 0; i < ncoefs; i++) {
+ if (i > 0) {
+ if (','!=*next) {
+ pj_log_error (P, "Horner: Malformed polynomium set %s. need %d coefs", param, ncoefs);
+ return 0;
+ }
+ init = ++next;
+ }
+ coefs[i] = pj_strtod (init, &next);
+ }
+ return 1;
+}
+
+
+/*********************************************************************/
+PJ *PROJECTION(horner) {
+/*********************************************************************/
+ int degree = 0, n;
+ HORNER *Q;
+ P->fwdobs = horner_forward_obs;
+ P->invobs = horner_reverse_obs;
+ P->fwd3d = 0;
+ P->inv3d = 0;
+ P->fwd = 0;
+ P->inv = 0;
+ P->left = P->right = PJ_IO_UNITS_METERS;
+ /* silence a few compiler warnings */
+ horner_silence (0);
+
+ /* Polynomial degree specified? */
+ if (pj_param (P->ctx, P->params, "tdeg").i) /* degree specified? */
+ degree = pj_param(P->ctx, P->params, "ideg").i;
+ else {
+ pj_log_debug (P, "Horner: Need to specify polynomial degree, (+deg=n)");
+ return horner_freeup (P);
+ }
+
+ Q = horner_alloc (degree);
+ P->opaque = (void *) Q;
+
+ n = horner_number_of_coefficients (degree);
+
+ if (0==parse_coefs (P, Q->fwd_u, "fwd_u", n))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, Q->fwd_v, "fwd_v", n))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, Q->inv_u, "inv_u", n))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, Q->inv_v, "inv_v", n))
+ return horner_freeup (P);
+
+ if (0==parse_coefs (P, (double *)(Q->fwd_origin), "fwd_origin", 2))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2))
+ return horner_freeup (P);
+ if (0==parse_coefs (P, &Q->range, "range", 1))
+ Q->range = 500000;
+
+ return P;
+}
+
+#ifndef PJ_SELFTEST
+/* selftest stub */
+int pj_horner_selftest (void) {return 0;}
+#else
+
+char tc32_utm32[] = {
+ " +proj=horner"
+ " +ellps=intl"
+ " +range=500000"
+ " +fwd_origin=877605.269066,6125810.306769"
+ " +inv_origin=877605.760036,6125811.281773"
+ " +deg=4"
+ " +fwd_v=6.1258112678e+06,9.9999971567e-01,1.5372750011e-10,5.9300860915e-15,2.2609497633e-19,4.3188227445e-05,2.8225130416e-10,7.8740007114e-16,-1.7453997279e-19,1.6877465415e-10,-1.1234649773e-14,-1.7042333358e-18,-7.9303467953e-15,-5.2906832535e-19,3.9984284847e-19"
+ " +fwd_u=8.7760574982e+05,9.9999752475e-01,2.8817299305e-10,5.5641310680e-15,-1.5544700949e-18,-4.1357045890e-05,4.2106213519e-11,2.8525551629e-14,-1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1.2429019719e-15,5.3886155968e-19,-1.0167505000e-18"
+ " +inv_v=6.1258103208e+06,1.0000002826e+00,-1.5372762184e-10,-5.9304261011e-15,-2.2612705361e-19,-4.3188331419e-05,-2.8225549995e-10,-7.8529116371e-16,1.7476576773e-19,-1.6875687989e-10,1.1236475299e-14,1.7042518057e-18,7.9300735257e-15,5.2881862699e-19,-3.9990736798e-19"
+ " +inv_u=8.7760527928e+05,1.0000024735e+00,-2.8817540032e-10,-5.5627059451e-15,1.5543637570e-18,4.1357152105e-05,-4.2114813612e-11,-2.8523713454e-14,1.9109017837e-18,-3.3616407783e-10,-2.4382678126e-14,2.0245020199e-18,-1.2441377565e-15,-5.3885232238e-19,1.0167203661e-18"
+};
+
+
+int pj_horner_selftest (void) {
+ PJ *P;
+ PJ_OBS a, b;
+ double dist;
+
+ /* The TC32 for "System 45 Bornholm" */
+ /* TC32 -> UTM32" */
+ P = pj_create (tc32_utm32);
+ if (0==P)
+ return 10;
+
+ a = b = pj_obs_null;
+ a.coo.uv.v = 6125305.4245;
+ a.coo.uv.u = 878354.8539;
+
+ /* Forward projection */
+ b = pj_trans (P, PJ_FWD, a);
+ b = pj_trans (P, PJ_INV, b);
+
+ /* Check roundtrip precision for 1 iteration each way */
+ dist = pj_roundtrip (P, PJ_FWD, 1, a);
+ if (dist > 0.01)
+ return 1;
+ return 0;
+}
+#endif
diff --git a/src/PJ_krovak.c b/src/PJ_krovak.c
index 136978a0..c469346f 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -239,7 +239,7 @@ int pj_krovak_selftest (void) {
double tolerance_lp = 1e-10;
double tolerance_xy = 1e-7;
- char e_args[] = {"+proj=krovak +ellps=GRS80"};
+ char e_args[] = {"+proj=krovak +ellps=GRS80 +no_defs"};
LP fwd_in[] = {
{ 2, 1},
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index ef86ccd3..27fe6ec9 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -114,6 +114,8 @@ struct pj_opaque {
int *reverse_step;
int *omit_forward;
int *omit_inverse;
+ char **argv;
+ char **current_argv;
PJ_OBS stack[PIPELINE_STACK_SIZE];
PJ **pipeline;
};
@@ -175,13 +177,17 @@ static LP pipeline_reverse (XY xyz, PJ *P);
static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
- int i, first_step, last_step, incr;
+ int i, first_step, last_step;
first_step = 1;
last_step = P->opaque->steps + 1;
- incr = 1;
- for (i = first_step; i != last_step; i += incr) {
+ for (i = first_step; i != last_step; i++) {
+ pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s",
+ i-first_step, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
+ P->opaque->pipeline[i]->descr
+ );
+
if (P->opaque->omit_forward[i])
continue;
if (P->opaque->reverse_step[i])
@@ -191,19 +197,23 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
if (P->opaque->depth < PIPELINE_STACK_SIZE)
P->opaque->stack[P->opaque->depth++] = point;
}
+ pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
P->opaque->depth = 0; /* Clear the stack */
return point;
}
+
static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
- int i, first_step, last_step, incr;
+ int i, first_step, last_step;
first_step = P->opaque->steps;
last_step = 0;
- incr = -1;
-
- for (i = first_step; i != last_step; i += incr) {
+ for (i = first_step; i != last_step; i--) {
+ pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f",
+ i - 1, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
+ P->opaque->pipeline[i]->a, P->opaque->pipeline[i]->rf
+ );
if (P->opaque->omit_inverse[i])
continue;
if (P->opaque->reverse_step[i])
@@ -213,6 +223,7 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
if (P->opaque->depth < PIPELINE_STACK_SIZE)
P->opaque->stack[P->opaque->depth++] = point;
}
+ pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
P->opaque->depth = 0; /* Clear the stack */
return point;
@@ -264,7 +275,7 @@ static void *pipeline_freeup (PJ *P, int errlev) { /* Destructor */
if (0==P)
return 0;
- pj_error_set (P, errlev);
+ pj_err_level (P, errlev);
if (0==P->opaque)
return pj_dealloc (P);
@@ -275,6 +286,8 @@ static void *pipeline_freeup (PJ *P, int errlev) { /* Destructor */
pj_dealloc (P->opaque->reverse_step);
pj_dealloc (P->opaque->omit_forward);
pj_dealloc (P->opaque->omit_inverse);
+ pj_dealloc (P->opaque->argv);
+ pj_dealloc (P->opaque->current_argv);
pj_dealloc (P->opaque->pipeline);
pj_dealloc (P->opaque);
@@ -313,8 +326,9 @@ static PJ *pj_create_pipeline (PJ *P, size_t steps) {
return P;
}
+
/* count the number of args in pipeline definition */
-size_t argc_params (paralist *params) {
+static size_t argc_params (paralist *params) {
size_t argc = 0;
for (; params != 0; params = params->next)
argc++;
@@ -322,18 +336,18 @@ size_t argc_params (paralist *params) {
}
/* Sentinel for argument list */
-static char argv_sentinel[5] = "step";
+static char argv_sentinel[] = "step";
-/* turn paralist int argc/argv style argument list */
-char **argv_params (paralist *params) {
+/* turn paralist into argc/argv style argument list */
+static char **argv_params (paralist *params, size_t argc) {
char **argv;
- size_t argc = 0;
- argv = pj_calloc (argc_params (params), sizeof (char *));
+ size_t i = 0;
+ argv = pj_calloc (argc, sizeof (char *));
if (0==argv)
return 0;
for (; params != 0; params = params->next)
- argv[argc++] = params->param;
- argv[argc++] = argv_sentinel;
+ argv[i++] = params->param;
+ argv[i++] = argv_sentinel;
return argv;
}
@@ -356,17 +370,15 @@ PJ *PROJECTION(pipeline) {
return 0;
argc = argc_params (P->params);
- argv = argv_params (P->params);
+ P->opaque->argv = argv = argv_params (P->params, argc);
if (0==argv)
return pipeline_freeup (P, ENOMEM);
- /* The elements of current_argv are not used - we just use argv_params */
- /* as allocator for a "large enough" container needed later */
- current_argv = argv_params (P->params);
+ P->opaque->current_argv = current_argv = pj_calloc (argc, sizeof (char *));
if (0==current_argv)
return pipeline_freeup (P, ENOMEM);
- /* Do some syntactic sanity checking */
+ /* Do some syntactical sanity checking */
for (i = 0; i < argc; i++) {
if (0==strcmp ("step", argv[i])) {
if (-1==i_pipeline) {
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 9e1ac4f7..2fcc4e31 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -77,6 +77,8 @@ SET(SRC_LIBPROJ_PJ
PJ_gstmerc.c
PJ_hammer.c
PJ_hatano.c
+ PJ_helmert.c
+ PJ_horner.c
PJ_igh.c
PJ_isea.c
PJ_imw_p.c
diff --git a/src/makefile.vc b/src/makefile.vc
index a2869217..e1eacf04 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -59,7 +59,7 @@ support = \
pj_strtod.obj pj_run_selftests.obj pj_generic_selftest.obj
pipeline = \
- pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj
+ pj_obs_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj
geodesic = geodesic.obj
@@ -147,6 +147,7 @@ install: all
copy *.exe $(INSTDIR)\bin
copy *.dll $(INSTDIR)\bin
copy *.lib $(INSTDIR)\lib
+ copy proj.h $(INSTDIR)\include
copy proj_api.h $(INSTDIR)\include
copy projects.h $(INSTDIR)\include
copy geodesic.h $(INSTDIR)\include
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 469e9a47..598c6dd2 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -6,41 +6,52 @@
XY /* forward projection entry */
pj_fwd(LP lp, PJ *P) {
XY xy;
- double t;
-
- /* check for forward and latitude or longitude overange */
- if ((t = fabs(lp.phi)-M_HALFPI) > EPS || fabs(lp.lam) > 10.) {
- xy.x = xy.y = HUGE_VAL;
- pj_ctx_set_errno( P->ctx, -14);
- } else { /* proceed with projection */
- P->ctx->last_errno = 0;
- pj_errno = 0;
- errno = 0;
-
- if (fabs(t) <= EPS)
- lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
- else if (P->geoc)
- lp.phi = atan(P->rone_es * tan(lp.phi));
- lp.lam -= P->lam0; /* compute del lp.lam */
- if (!P->over)
- lp.lam = adjlon(lp.lam); /* adjust del longitude */
-
- /* Check for NULL pointer */
- if (P->fwd != NULL)
- {
- xy = (*P->fwd)(lp, P); /* project */
- if ( P->ctx->last_errno )
- xy.x = xy.y = HUGE_VAL;
- /* adjust for major axis and easting/northings */
- else {
- xy.x = P->fr_meter * (P->a * xy.x + P->x0);
- xy.y = P->fr_meter * (P->a * xy.y + P->y0);
- }
- }
- else
- {
- xy.x = xy.y = HUGE_VAL;
- }
- }
- return xy;
+ XY err;
+ double t;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.x = err.y = HUGE_VAL;
+
+ if (0==P->fwd)
+ return err;
+
+ P->ctx->last_errno = pj_errno = errno = 0;
+
+ /* Check input coordinates if angular */
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+
+ /* check for forward and latitude or longitude overange */
+ t = fabs(lp.phi)-M_HALFPI;
+ if (t > EPS || fabs(lp.lam) > 10.) {
+ pj_ctx_set_errno( P->ctx, -14);
+ return err;
+ }
+
+ /* Clamp latitude to -90..90 degree range */
+ if (fabs(t) <= EPS)
+ lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
+ else if (P->geoc) /* Maybe redundant and never used. */
+ lp.phi = atan(P->rone_es * tan(lp.phi));
+ lp.lam -= P->lam0; /* compute del lp.lam */
+ if (!P->over)
+ lp.lam = adjlon(lp.lam); /* adjust del longitude */
+ }
+
+ /* Do the transformation */
+ xy = (*P->fwd)(lp, P);
+ if ( P->ctx->last_errno )
+ return err;
+
+ /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xy.x *= P->a;
+ xy.y *= P->a;
+ }
+
+ /* Handle false eastings/northings and non-metric linear units */
+ xy.x = P->fr_meter * (xy.x + P->x0);
+ xy.y = P->fr_meter * (xy.y + P->y0);
+ /* z is not scaled since this is handled by vto_meter outside */
+
+ return xy;
}
diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c
index 63fc529e..0fb2f633 100644
--- a/src/pj_fwd3d.c
+++ b/src/pj_fwd3d.c
@@ -1,47 +1,59 @@
-/* general forward projection */
#define PJ_LIB__
#include <projects.h>
#include <errno.h>
# define EPS 1.0e-12
- XYZ /* forward projection entry */
-pj_fwd3d(LPZ lpz, PJ *P) {
- XYZ xyz;
- double t;
-
- /* check for forward and latitude or longitude overange */
- if ((t = fabs(lpz.phi)-M_HALFPI) > EPS || fabs(lpz.lam) > 10.) {
- xyz.x = xyz.y = xyz.z = HUGE_VAL;
- pj_ctx_set_errno( P->ctx, -14);
- } else { /* proceed with projection */
- P->ctx->last_errno = 0;
- pj_errno = 0;
- errno = 0;
-
- if (fabs(t) <= EPS)
- lpz.phi = lpz.phi < 0. ? -M_HALFPI : M_HALFPI;
- else if (P->geoc) /* Maybe redundant and never used. */
- lpz.phi = atan(P->rone_es * tan(lpz.phi));
- lpz.lam -= P->lam0; /* compute del lp.lam */
- if (!P->over)
- lpz.lam = adjlon(lpz.lam); /* adjust del longitude */
-
- /* Check for NULL pointer */
- if (P->fwd3d != NULL)
- {
- xyz = (*P->fwd3d)(lpz, P); /* project */
- if ( P->ctx->last_errno )
- xyz.x = xyz.y = xyz.z = HUGE_VAL;
- /* adjust for major axis and easting/northings */
- else {
- xyz.x = P->fr_meter * (P->a * xyz.x + P->x0);
- xyz.y = P->fr_meter * (P->a * xyz.y + P->y0);
- /* z is not scaled since this handled by vto_meter outside */
- }
- }
- else
- {
- xyz.x = xyz.y = xyz.z = HUGE_VAL;
- }
- }
- return xyz;
+
+/* 3D Forward transformation */
+
+XYZ pj_fwd3d(LPZ lpz, PJ *P) {
+ XYZ xyz;
+ XYZ err;
+ double t;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.x = err.y = err.z = HUGE_VAL;
+
+ if (0==P->fwd3d)
+ return err;
+
+
+ P->ctx->last_errno = pj_errno = errno = 0;
+
+ /* Check input coordinates if angular */
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+
+ /* check for forward and latitude or longitude overange */
+ t = fabs(lpz.phi)-M_HALFPI;
+ if (t > EPS || fabs(lpz.lam) > 10.) {
+ pj_ctx_set_errno( P->ctx, -14);
+ return err;
+ }
+
+ /* Clamp latitude to -90..90 degree range */
+ if (fabs(t) <= EPS)
+ lpz.phi = lpz.phi < 0. ? -M_HALFPI : M_HALFPI;
+ else if (P->geoc) /* Maybe redundant and never used. */
+ lpz.phi = atan(P->rone_es * tan(lpz.phi));
+ lpz.lam -= P->lam0; /* compute del lp.lam */
+ if (!P->over)
+ lpz.lam = adjlon(lpz.lam); /* adjust del longitude */
+ }
+
+ /* Do the transformation */
+ xyz = (*P->fwd3d)(lpz, P);
+ if ( P->ctx->last_errno )
+ return err;
+
+ /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xyz.x *= P->a;
+ xyz.y *= P->a;
+ }
+
+ /* Handle false eastings/northings and non-metric linear units */
+ xyz.x = P->fr_meter * (xyz.x + P->x0);
+ xyz.y = P->fr_meter * (xyz.y + P->y0);
+ /* z is not scaled since this is handled by vto_meter outside */
+
+ return xyz;
}
diff --git a/src/pj_init.c b/src/pj_init.c
index c468f1d7..b6be52d9 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -35,10 +35,12 @@
#include <errno.h>
#include <ctype.h>
+/* Maximum size of files using the "escape carriage return" feature */
+#define MAX_CR_ESCAPE 65537
typedef struct {
projCtx ctx;
PAFile fid;
- char buffer[8193];
+ char buffer[MAX_CR_ESCAPE];
int buffer_filled;
int at_eof;
} pj_read_state;
@@ -51,6 +53,7 @@ static const char *fill_buffer(pj_read_state *state, const char *last_char)
{
size_t bytes_read;
size_t char_remaining, char_requested;
+ char *r, *w;
/* -------------------------------------------------------------------- */
/* Don't bother trying to read more if we are at eof, or if the */
@@ -85,7 +88,36 @@ static const char *fill_buffer(pj_read_state *state, const char *last_char)
state->buffer[state->buffer_filled + bytes_read] = '\0';
}
- state->buffer_filled += bytes_read;
+/* -------------------------------------------------------------------- */
+/* Line continuations: skip whitespace after escaped newlines */
+/* -------------------------------------------------------------------- */
+ r = state->buffer;
+ w = state->buffer;
+ while (*r) {
+ /* Escaped newline? */
+ while ((r[0]=='\\') && ((r[1]=='\n') || (r[1]=='\r'))) {
+ r += 2;
+ while (isspace (*r))
+ r++;
+ /* we also skip comments immediately after an escaped newline */
+ while (*r=='#') {
+ while( *r && (*r != '\n') )
+ r++;
+ while (isspace (*r))
+ r++;
+ /* Reaching end of buffer while skipping continuation comment is currently an error */
+ if (0==*r) {
+ pj_ctx_set_errno (state->ctx, -2);
+ pj_log (state->ctx, PJ_LOG_ERROR, "init file too big");
+ return 0;
+ }
+ }
+ }
+ *w++ = *r++;
+ }
+ *w = 0;
+ state->buffer_filled += (bytes_read - (r-w));
+
return last_char;
}
@@ -96,11 +128,11 @@ static paralist *
get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
int *found_def) {
pj_read_state *state = (pj_read_state*) calloc(1,sizeof(pj_read_state));
- char sword[302];
+ char sword[MAX_CR_ESCAPE];
+ char *pipeline;
int len;
int in_target = 0;
const char *next_char = NULL;
-
state->fid = fid;
state->ctx = ctx;
next_char = fill_buffer(state, NULL);
@@ -110,6 +142,9 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
len = strlen(name);
*sword = 't';
+ if (0==next_char)
+ return 0;
+
/* loop till we find our target keyword */
while (*next_char)
{
@@ -120,6 +155,8 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
next_char++;
next_char = fill_buffer(state, next_char);
+ if (0==next_char)
+ return 0;
/* for comments, skip past end of line. */
if( *next_char == '#' )
@@ -128,6 +165,8 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
next_char++;
next_char = fill_buffer(state, next_char);
+ if (0==next_char)
+ return 0;
if (*next_char == '\n')
next_char++;
if (*next_char == '\r')
@@ -171,7 +210,7 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
}
/* capture parameter */
- while( *next_char && !isspace(*next_char) )
+ while ( *next_char && !isspace(*next_char) )
{
next_char++;
word_len++;
@@ -180,8 +219,11 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
strncpy(sword+1, start_of_word, word_len);
sword[word_len+1] = '\0';
- /* do not override existing parameter value of same name */
- if (!pj_param(ctx, *start, sword).i) {
+ /* do not override existing parameter value of same name - unless in pipeline definition */
+ pipeline = pj_param(ctx, *start, "sproj").s;
+ if (pipeline && strcmp (pipeline, "pipeline"))
+ pipeline = 0;
+ if (pipeline || !pj_param(ctx, *start, sword).i) {
/* don't default ellipse if datum, ellps or any earth model
information is set. */
if( strncmp(sword+1,"ellps=",6) != 0
@@ -194,8 +236,9 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
{
next = next->next = pj_mkparam(sword+1);
}
+ else
+ next = next->next = pj_mkparam(sword+1);
}
-
}
else
{
@@ -388,6 +431,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
PJ *(*proj)(PJ *);
paralist *curr;
int i;
+ int defer_init_expansion = 0;
PJ *PIN = 0;
ctx->last_errno = 0;
@@ -403,9 +447,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
if (ctx->last_errno) goto bum_call;
/* check if +init present */
- if (pj_param(ctx, start, "tinit").i) {
+ if (pj_param(ctx, start, "tinit").i && ! defer_init_expansion) {
int found_def = 0;
-
+ /* avoid expanding additional inits (as could happen in a pipeline) */
+ defer_init_expansion = 1;
if (!(curr = get_init(ctx,&start, curr,
pj_param(ctx, start, "sinit").s,
&found_def)))
@@ -444,7 +489,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
if (pj_datum_set(ctx, start, PIN)) goto bum_call;
/* set ellipsoid/sphere parameters */
- if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) goto bum_call;
+ if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) {
+ pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
+ goto bum_call;
+ }
PIN->a_orig = PIN->a;
PIN->es_orig = PIN->es;
@@ -464,11 +512,11 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
/* flattening */
PIN->f = 1 - cos (PIN->alpha); /* = 1 - sqrt (1 - PIN->es); */
- PIN->rf = PIN->f? 1/PIN->f: HUGE_VAL;
+ PIN->rf = PIN->f? 1.0/PIN->f: HUGE_VAL;
/* second flattening */
- PIN->f2 = 1/cos (PIN->alpha) - 1;
- PIN->rf = PIN->f2? 1/PIN->f2: HUGE_VAL;
+ PIN->f2 = 1/cos (PIN->alpha) - 1;
+ PIN->rf2 = PIN->f2? 1/PIN->f2: HUGE_VAL;
/* third flattening */
PIN->n = pow (tan (PIN->alpha/2), 2);
diff --git a/src/pj_inv.c b/src/pj_inv.c
index ed05618a..d5393261 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -7,35 +7,49 @@
/* inverse projection entry */
LP pj_inv(XY xy, PJ *P) {
LP lp;
+ LP err;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.lam = err.phi = HUGE_VAL;
+
+ if (0==P->inv)
+ return err;
/* can't do as much preliminary checking as with forward */
if (xy.x == HUGE_VAL || xy.y == HUGE_VAL) {
- lp.lam = lp.phi = HUGE_VAL;
pj_ctx_set_errno( P->ctx, -15);
- return lp;
+ return err;
+ }
+ P->ctx->last_errno = errno = pj_errno = 0;
+
+ /* de-scale and de-offset */
+ xy.x = (xy.x * P->to_meter - P->x0);
+ xy.y = (xy.y * P->to_meter - P->y0);
+
+ /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
+ /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */
+ /* stomps on a and hence depends on this */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xy.x *= P->ra;
+ xy.y *= P->ra;
}
- errno = pj_errno = 0;
- P->ctx->last_errno = 0;
-
- /* descale and de-offset */
- xy.x = (xy.x * P->to_meter - P->x0) * P->ra;
- xy.y = (xy.y * P->to_meter - P->y0) * P->ra;
-
- /* Check for NULL pointer */
- if (P->inv != NULL) {
- lp = (*P->inv)(xy, P); /* inverse project */
- if (P->ctx->last_errno ) {
- lp.lam = lp.phi = HUGE_VAL;
- } else {
- lp.lam += P->lam0; /* reduce from del lp.lam */
- if (!P->over)
- lp.lam = adjlon(lp.lam); /* adjust longitude to CM */
- if (P->geoc && fabs(fabs(lp.phi)-M_HALFPI) > EPS)
- lp.phi = atan(P->one_es * tan(lp.phi));
- }
- } else {
- lp.lam = lp.phi = HUGE_VAL;
+ /* Do inverse transformation */
+ lp = (*P->inv) (xy, P);
+ if (P->ctx->last_errno)
+ return err;
+
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+ /* reduce from del lp.lam */
+ lp.lam += P->lam0;
+
+ /* adjust longitude to central meridian */
+ if (!P->over)
+ lp.lam = adjlon(lp.lam);
+
+ /* This may be redundant and never used */
+ if (P->geoc && fabs(fabs(lp.phi)-M_HALFPI) > EPS)
+ lp.phi = atan(P->one_es * tan(lp.phi));
}
return lp;
diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c
index bf6a9324..c6052ad4 100644
--- a/src/pj_inv3d.c
+++ b/src/pj_inv3d.c
@@ -1,45 +1,57 @@
-/* general inverse projection */
#define PJ_LIB__
#include <projects.h>
#include <errno.h>
# define EPS 1.0e-12
- LPZ /* inverse projection entry */
-pj_inv3d(XYZ xyz, PJ *P) {
- LPZ lpz;
-
- /* can't do as much preliminary checking as with forward */
- if (xyz.x == HUGE_VAL || xyz.y == HUGE_VAL || xyz.z == HUGE_VAL ) {
- lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
- pj_ctx_set_errno( P->ctx, -15);
- return lpz;
- }
-
- errno = pj_errno = 0;
- P->ctx->last_errno = 0;
-
- xyz.x = (xyz.x * P->to_meter - P->x0) * P->ra; /* descale and de-offset */
- xyz.y = (xyz.y * P->to_meter - P->y0) * P->ra;
- /* z is not scaled since that is handled by vto_meter before we get here */
-
- /* Check for NULL pointer */
- if (P->inv3d != NULL)
- {
- lpz = (*P->inv3d)(xyz, P); /* inverse project */
- if (P->ctx->last_errno )
- lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
- else {
- lpz.lam += P->lam0; /* reduce from del lp.lam */
- if (!P->over)
- lpz.lam = adjlon(lpz.lam); /* adjust longitude to CM */
-
- /* This may be redundant and never used */
- if (P->geoc && fabs(fabs(lpz.phi)-M_HALFPI) > EPS)
- lpz.phi = atan(P->one_es * tan(lpz.phi));
- }
- }
- else
- {
- lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
- }
- return lpz;
+
+/* 3D inverse transformation */
+
+LPZ pj_inv3d (XYZ xyz, PJ *P) {
+ LPZ lpz;
+ LPZ err;
+
+ /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
+ err.lam = err.phi = err.z = HUGE_VAL;
+
+ if (0==P->inv3d)
+ return err;
+
+ /* can't do as much preliminary checking as with forward */
+ if (xyz.x == HUGE_VAL || xyz.y == HUGE_VAL || xyz.z == HUGE_VAL ) {
+ pj_ctx_set_errno( P->ctx, -15);
+ return err;
+ }
+
+ P->ctx->last_errno = errno = pj_errno = 0;
+
+ /* de-scale and de-offset */
+ /* z is not de-scaled since that is handled by vto_meter before we get here */
+ xyz.x = (xyz.x * P->to_meter - P->x0);
+ xyz.y = (xyz.y * P->to_meter - P->y0);
+ /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
+ /* Multiplying by ra, rather than dividing by a because the CALCOFI projection */
+ /* stomps on a and hence depends on this */
+ if (P->right==PJ_IO_UNITS_CLASSIC) {
+ xyz.x *= P->ra;
+ xyz.y *= P->ra;
+ }
+
+ /* Do inverse transformation */
+ lpz = (*P->inv3d) (xyz, P);
+ if (P->ctx->last_errno)
+ return err;
+
+ if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+ /* reduce from del lp.lam */
+ lpz.lam += P->lam0;
+
+ /* adjust longitude to central meridian */
+ if (!P->over)
+ lpz.lam = adjlon(lpz.lam);
+
+ /* This may be redundant and never used */
+ if (P->geoc && fabs(fabs(lpz.phi)-M_HALFPI) > EPS)
+ lpz.phi = atan(P->one_es * tan(lpz.phi));
+ }
+
+ return lpz;
}
diff --git a/src/pj_list.h b/src/pj_list.h
index f49fd574..96db13d7 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -52,6 +52,8 @@ PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff")
PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area")
PROJ_HEAD(healpix, "HEALPix")
PROJ_HEAD(rhealpix, "rHEALPix")
+PROJ_HEAD(helmert, "3- and 7-parameter Helmert shift")
+PROJ_HEAD(horner, "Horner polynomial evaluation")
PROJ_HEAD(igh, "Interrupted Goode Homolosine")
PROJ_HEAD(imw_p, "International Map of the World Polyconic")
PROJ_HEAD(isea, "Icosahedral Snyder Equal Area")
diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
index 716c6095..4a6e98a5 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -32,7 +32,7 @@
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
*****************************************************************************/
-#define PJ_OBS_C
+#define PJ_OBS_API_C
#include <proj.h>
#include <projects.h>
#include <geodesic.h>
@@ -40,16 +40,9 @@
#include <math.h>
-/* Used as return value in case of errors */
-const PJ_OBS pj_obs_error = {
- /* Cannot use HUGE_VAL here: MSVC misimplements HUGE_VAL as something that is not compile time constant */
- {{DBL_MAX,DBL_MAX,DBL_MAX,DBL_MAX}},
- {{DBL_MAX,DBL_MAX,DBL_MAX}},
- 0, 0
-};
-
/* Used for zero-initializing new objects */
-const PJ_OBS pj_obs_null = {
+const PJ_COORD pj_coo_null = {{0, 0, 0, 0}};
+const PJ_OBS pj_obs_null = {
{{0, 0, 0, 0}},
{{0, 0, 0}},
0, 0
@@ -77,6 +70,23 @@ double pj_xyz_dist (XYZ a, XYZ b) {
}
+/* Work around non-constness of MSVC HUGE_VAL by providing functions rather than constants */
+PJ_COORD pj_coo_error (void) {
+ PJ_COORD c;
+ c.v[0] = c.v[1] = c.v[2] = c.v[3] = HUGE_VAL;
+ return c;
+}
+
+PJ_OBS pj_obs_error (void) {
+ PJ_OBS obs;
+ obs.coo = pj_coo_error ();
+ obs.anc.v[0] = obs.anc.v[1] = obs.anc.v[2] = HUGE_VAL;
+ obs.id = obs.flags = 0;
+ return obs;
+}
+
+
+
PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
if (0!=P->fwdobs) {
obs = P->fwdobs (obs, P);
@@ -90,8 +100,8 @@ PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
obs.coo.xy = pj_fwd (obs.coo.lp, P);
return obs;
}
- pj_error_set (P, EINVAL);
- return pj_obs_error;
+ pj_err_level (P, EINVAL);
+ return pj_obs_error ();
}
@@ -108,8 +118,8 @@ PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) {
obs.coo.lp = pj_inv (obs.coo.xy, P);
return obs;
}
- pj_error_set (P, EINVAL);
- return pj_obs_error;
+ pj_err_level (P, EINVAL);
+ return pj_obs_error ();
}
@@ -129,8 +139,8 @@ PJ_OBS pj_trans (PJ *P, enum pj_direction direction, PJ_OBS obs) {
break;
}
- pj_error_set (P, EINVAL);
- return pj_obs_error;
+ pj_err_level (P, EINVAL);
+ return pj_obs_error ();
}
@@ -143,7 +153,7 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
return HUGE_VAL;
if (n < 1) {
- pj_error_set (P, EINVAL);
+ pj_err_level (P, EINVAL);
return HUGE_VAL;
}
@@ -152,15 +162,15 @@ double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
for (i = 0; i < n; i++) {
switch (direction) {
case PJ_FWD:
- u.coo.xyz = pj_fwd3d (o.coo.lpz, P);
- o.coo.lpz = pj_inv3d (u.coo.xyz, P);
+ u = pj_fwdobs (o, P);
+ o = pj_invobs (u, P);
break;
case PJ_INV:
- u.coo.lpz = pj_inv3d (o.coo.xyz, P);
- o.coo.xyz = pj_fwd3d (u.coo.lpz, P);
+ u = pj_invobs (o, P);
+ o = pj_fwdobs (u, P);
break;
default:
- pj_error_set (P, EINVAL);
+ pj_err_level (P, EINVAL);
return HUGE_VAL;
}
}
@@ -179,18 +189,24 @@ PJ *pj_create_argv (int argc, char **argv) {
}
-/* From here: Minimum viable support for contexts. The first four functions */
+/* Below: Minimum viable support for contexts. The first four functions */
/* relate to error reporting, debugging, and logging, hence being generically */
/* useful. The remaining is a compact implementation of the more low level */
/* proj_api.h thread contexts, which may or may not be useful */
-int pj_error (PJ *P) {
- return pj_ctx_get_errno (pj_get_ctx(P));
-}
-
-
-void pj_error_set (PJ *P, int err) {
- pj_ctx_set_errno (pj_get_ctx(P), err);
+/* Set error level 0-3, or query current level */
+int pj_err_level (PJ *P, int err_level) {
+ int previous;
+ PJ_CONTEXT *ctx;
+ if (0==P)
+ ctx = pj_get_default_ctx();
+ else
+ ctx = pj_get_ctx (P);
+ previous = pj_ctx_get_errno (ctx);
+ if (PJ_ERR_TELL==err_level)
+ return previous;
+ pj_ctx_set_errno (ctx, err_level);
+ return previous;
}
@@ -211,7 +227,7 @@ enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level) {
/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */
-void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *)) {
+void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *)) {
PJ_CONTEXT *ctx = pj_get_ctx (P);
ctx->app_data = app_data;
if (0!=log)
@@ -223,7 +239,7 @@ void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *))
int pj_context_renew (PJ *P) {
PJ_CONTEXT *ctx = pj_ctx_alloc ();
if (0==ctx) {
- pj_error_set (P, ENOMEM);
+ pj_err_level (P, ENOMEM);
return 1;
}
@@ -231,6 +247,7 @@ int pj_context_renew (PJ *P) {
return 0;
}
+
/* Move daughter to mother's context */
void pj_context_inherit (PJ *mother, PJ *daughter) {
pj_set_ctx (daughter, pj_get_ctx (mother));
diff --git a/src/pj_param.c b/src/pj_param.c
index 434d4535..24c2354e 100644
--- a/src/pj_param.c
+++ b/src/pj_param.c
@@ -23,13 +23,13 @@ pj_mkparam(char *str) {
/* character in `opt' is a parameter type which can take the */
/* values: */
/* */
-/* `t' - test for presence, return TRUE/FALSE in PROJVALUE.i */
-/* `i' - integer value returned in PROJVALUE.i */
-/* `d' - simple valued real input returned in PROJVALUE.f */
+/* `t' - test for presence, return TRUE/FALSE in PROJVALUE.i */
+/* `i' - integer value returned in PROJVALUE.i */
+/* `d' - simple valued real input returned in PROJVALUE.f */
/* `r' - degrees (DMS translation applied), returned as */
-/* radians in PROJVALUE.f */
-/* `s' - string returned in PROJVALUE.s */
-/* `b' - test for t/T/f/F, return in PROJVALUE.i */
+/* radians in PROJVALUE.f */
+/* `s' - string returned in PROJVALUE.s */
+/* `b' - test for t/T/f/F, return in PROJVALUE.i */
/* */
/************************************************************************/
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index 8a2a9c4b..bef42882 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -82,4 +82,5 @@ char *pj_strerrno(int err) {
return note;
}
}
+ return 0;
}
diff --git a/src/proj.def b/src/proj.def
index 3a5e3de9..b01b846c 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -89,23 +89,24 @@ EXPORTS
geod_polygon_testpoint @87
geod_polygon_clear @88
pj_run_selftests @89
- pj_error @90
- pj_create @91
- pj_create_argv @92
- pj_trans @93
- pj_fwdobs @94
- pj_invobs @95
- pj_roundtrip @96
- pj_log_level @97
- pj_error_set @98
- pj_log_set @99
- pj_context_renew @100
- pj_context_inherit @101
- pj_context_free @102
- pj_fileapi_set @103
- pj_log_error @104
- pj_log_debug @105
- pj_log_trace @106
- pj_lp_dist @107
- pj_xy_dist @108
- pj_xyz_dist @109
+ pj_create @90
+ pj_create_argv @91
+ pj_trans @92
+ pj_fwdobs @93
+ pj_invobs @94
+ pj_roundtrip @95
+ pj_log_level @96
+ pj_err_level @97
+ pj_log_func @98
+ pj_context_renew @99
+ pj_context_inherit @100
+ pj_context_free @101
+ pj_fileapi_set @102
+ pj_log_error @103
+ pj_log_debug @104
+ pj_log_trace @105
+ pj_lp_dist @106
+ pj_xy_dist @107
+ pj_xyz_dist @108
+ pj_coo_error @109
+ pj_obs_error @110
diff --git a/src/proj.h b/src/proj.h
index 6c33a916..f6b3bd2a 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -271,10 +271,13 @@ double pj_xy_dist (XY a, XY b);
double pj_xyz_dist (XYZ a, XYZ b);
-#ifndef PJ_OBS_C
-extern const PJ_OBS pj_obs_error;
-extern const PJ_OBS pj_obs_null;
-extern const PJ *pj_shutdown;
+#ifndef PJ_OBS_API_C
+/* Part of MSVC workaround: Make pj_*_null look function-like for symmetry with pj_*_error */
+#define pj_coo_null(x) pj_coo_null
+#define pj_obs_null(x) pj_obs_null
+extern const PJ_COORD pj_coo_null;
+extern const PJ_OBS pj_obs_null;
+extern const PJ *pj_shutdown;
#endif
#ifndef TODEG
@@ -301,17 +304,19 @@ enum pj_log_level {
PJ_LOG_DEBUG_MINOR = 3 /* for proj_api.h compatibility */
};
+/* Set or read error level */
+#define PJ_ERR_TELL -56789
+int pj_err_level (PJ *P, int err_level);
+
/* Set logging level 0-3. Higher number means more debug info. 0 turns it off */
enum pj_log_level pj_log_level (PJ *P, enum pj_log_level log_level);
void pj_log_error (PJ *P, const char *fmt, ...);
void pj_log_debug (PJ *P, const char *fmt, ...);
void pj_log_trace (PJ *P, const char *fmt, ...);
+void pj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));
-void pj_error_set (PJ *P, int err);
-void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *));
-
/* Lower level functionality for handling thread contexts */
int pj_context_renew (PJ *P);
void pj_context_inherit (PJ *mother, PJ *daughter);
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 10:03:27 +0000