[proj] 01/02: Add upstream patches to fix regressions in 5.0.0.
Bas Couwenberg
sebastic at debian.org
Fri Mar 9 06:49:49 UTC 2018
This is an automated email from the git hooks/post-receive script.
sebastic pushed a commit to branch master
in repository proj.
commit c7d93d22ac109dd9b73efd6e64805e15c43f442f
Author: Bas Couwenberg <sebastic at xs4all.nl>
Date: Fri Mar 9 07:34:03 2018 +0100
Add upstream patches to fix regressions in 5.0.0.
- 0001-Revert-fix-to-22.patch
Fixes web mercator projection
- pr845_*-reintroduce-support-for-vertical-scaling.patch
Fixes +vdatum handling in cs2cs
(closes: #892062)
---
debian/changelog | 11 +
debian/patches/0001-Revert-fix-to-22.patch | 86 +++
...-reintroduce-support-for-vertical-scaling.patch | 798 +++++++++++++++++++++
debian/patches/series | 2 +
4 files changed, 897 insertions(+)
diff --git a/debian/changelog b/debian/changelog
index 4591853..aefb8c0 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,14 @@
+proj (5.0.0-3) UNRELEASED; urgency=medium
+
+ * Add upstream patches to fix regressions in 5.0.0:
+ - 0001-Revert-fix-to-22.patch
+ Fixes web mercator projection
+ - pr845_*-reintroduce-support-for-vertical-scaling.patch
+ Fixes +vdatum handling in cs2cs
+ (closes: #892062)
+
+ -- Bas Couwenberg <sebastic at debian.org> Fri, 09 Mar 2018 07:27:23 +0100
+
proj (5.0.0-2) unstable; urgency=medium
* Install docs in all packages.
diff --git a/debian/patches/0001-Revert-fix-to-22.patch b/debian/patches/0001-Revert-fix-to-22.patch
new file mode 100644
index 0000000..32d2c36
--- /dev/null
+++ b/debian/patches/0001-Revert-fix-to-22.patch
@@ -0,0 +1,86 @@
+Description: Revert fix to #22
+ The fix in #22 solved the problem at hand and doing what was expected
+ from the specified parameters. Unfortunately it also removed the slightly
+ hacky "feature" that makes the web mercator work in pj_transform. The
+ web mercator is special since the latitude is computed on the ellipsoid,
+ but behaves as if it was defined on a sphere. Hence it is problematic to
+ change the ellipsoid parameters when using the web mercator, even though
+ that is the geodetically correct thing to do. The web mercator is used in
+ more or less any web mapping application and is thus one of the most
+ frequently used transformations in PROJ. This justifies re-introducing
+ the minor bug reported in #22.
+ .
+ The problem will have to be taken care of properly when pj_transform
+ is removed from the library in favour of the transformation pipelines
+ based API.
+Author: Kristian Evers <kristianevers at gmail.com>
+Origin: https://github.com/OSGeo/proj.4/commit/f41fad3ac2bff401456f31dd3273e163ea7d09af
+
+--- a/nad/proj_outIGNF.dist
++++ b/nad/proj_outIGNF.dist
+@@ -1,16 +1,16 @@
+ +init=./IGNF:NTFG +to +init=./IGNF:RGF93G
+ 3.300866856 43.4477976569 0.0000 3d18'0.915"E 43d26'52.077"N 0.000
+ +init=./IGNF:LAMBE +to +init=./IGNF:LAMB93
+- 600000.0000 2600545.4523 0.0000 652760.737 7033791.244 0.000
++ 600000.0000 2600545.4523 0.0000 652760.737 7033791.243 0.000
+ 135638.3592 2418760.4094 0.0000 187194.062 6855928.882 0.000
+ 998137.3947 2413822.2844 0.0000 1049052.258 6843776.562 0.000
+- 600000.0000 2200000.0000 0.0000 649398.872 6633524.192 0.000
++ 600000.0000 2200000.0000 0.0000 649398.872 6633524.191 0.000
+ 311552.5340 1906457.4840 0.0000 358799.172 6342652.486 0.000
+ 960488.4138 1910172.8812 0.0000 1007068.686 6340907.237 0.000
+ 600000.0000 1699510.8340 0.0000 645204.279 6133556.746 0.000
+-1203792.5981 626873.17210 0.0000 1238875.764 5057405.017 0.000
++1203792.5981 626873.17210 0.0000 1238875.764 5057405.016 0.000
+ +init=./IGNF:LAMBE +to +init=./IGNF:GEOPORTALFXX
+- 600000.0000 2600545.4523 0.0000 179040.148 5610495.276 0.000
++ 600000.0000 2600545.4523 0.0000 179040.148 5610495.275 0.000
+ 135638.3592 2418760.4094 0.0000 -303729.363 5410118.356 0.000
+ 998137.3947 2413822.2844 0.0000 592842.792 5410120.554 0.000
+ 600000.0000 2200000.0000 0.0000 179041.670 5209746.080 0.000
+@@ -37,4 +37,4 @@
+ 2d20'11.7754730" 42d18'00.0824436" 0.0 260109.601 5009175.714 0.000
+ 9d32'12.6680218" 41d24'00.3542556" 0.0 1061637.534 4889066.592 0.000
+ +init=./IGNF:RGR92 +to +init=./IGNF:REUN47
+-3356123.5400 1303218.3090 5247430.6050 3353421.833 1304074.314 5248935.606
++3356123.5400 1303218.3090 5247430.6050 3353421.833 1304074.314 5248935.607
+--- a/src/pj_transform.c
++++ b/src/pj_transform.c
+@@ -748,32 +748,17 @@ int pj_datum_transform( PJ *srcdefn, PJ
+ /* -------------------------------------------------------------------- */
+ if( srcdefn->datum_type == PJD_GRIDSHIFT )
+ {
+- const char* srcnadgrids = pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s;
+-
+ pj_apply_gridshift_2( srcdefn, 0, point_count, point_offset, x, y, z );
+ CHECK_RETURN(srcdefn);
+
+- /* If the gridlist has either "@null" or "null" as its only */
+- /* grid we don't change the ellipsoid parameters, since the */
+- /* datum shift to WGS84 was not performed in practice. */
+- if ( srcnadgrids != NULL &&
+- strcmp("@null", srcnadgrids) && strcmp("null", srcnadgrids) ) {
+- src_a = SRS_WGS84_SEMIMAJOR;
+- src_es = SRS_WGS84_ESQUARED;
+- }
++ src_a = SRS_WGS84_SEMIMAJOR;
++ src_es = SRS_WGS84_ESQUARED;
+ }
+
+ if( dstdefn->datum_type == PJD_GRIDSHIFT )
+ {
+- const char* dstnadgrids = pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s;
+- /* If the gridlist has either "@null" or "null" as its only */
+- /* grid we don't change the ellipsoid parameters, since the */
+- /* datum shift to WGS84 will not be performed. */
+- if ( dstnadgrids != NULL &&
+- strcmp("@null", dstnadgrids) && strcmp("null", dstnadgrids) ) {
+- dst_a = SRS_WGS84_SEMIMAJOR;
+- dst_es = SRS_WGS84_ESQUARED;
+- }
++ dst_a = SRS_WGS84_SEMIMAJOR;
++ dst_es = SRS_WGS84_ESQUARED;
+ }
+
+ /* ==================================================================== */
diff --git a/debian/patches/pr845_Refactor-pj_transform-reintroduce-support-for-vertical-scaling.patch b/debian/patches/pr845_Refactor-pj_transform-reintroduce-support-for-vertical-scaling.patch
new file mode 100644
index 0000000..b09de2d
--- /dev/null
+++ b/debian/patches/pr845_Refactor-pj_transform-reintroduce-support-for-vertical-scaling.patch
@@ -0,0 +1,798 @@
+Subject: Refactor pj_transform, reintroduce support for vertical scaling
+ Change operation step names to reflect their bidirectional nature
+Author: Thomas Knudsen <thokn at sdfe.dk>
+Origin: https://github.com/OSGeo/proj.4/pull/845
+Bug: https://github.com/OSGeo/proj.4/issues/833
+Bug-Debian: https://bugs.debian.org/892062
+
+--- a/src/pj_transform.c
++++ b/src/pj_transform.c
+@@ -27,11 +27,23 @@
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+-#include <projects.h>
++#include "projects.h"
+ #include <string.h>
+ #include <math.h>
+ #include "geocent.h"
+
++
++/* Apply transformation to observation - in forward or inverse direction */
++/* Copied from proj.h */
++enum PJ_DIRECTION {
++ PJ_FWD = 1, /* Forward */
++ PJ_IDENT = 0, /* Do nothing */
++ PJ_INV = -1 /* Inverse */
++};
++typedef enum PJ_DIRECTION PJ_DIRECTION;
++
++
++
+ static int pj_adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
+ long point_count, int point_offset,
+ double *x, double *y, double *z );
+@@ -81,388 +93,468 @@ static const int transient_error[60] = {
+ /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
+ /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0 };
+
+-/************************************************************************/
+-/* pj_transform() */
+-/* */
+-/* Currently this function doesn't recognise if two projections */
+-/* are identical (to short circuit reprojection) because it is */
+-/* difficult to compare PJ structures (since there are some */
+-/* projection specific components). */
+-/************************************************************************/
+-
+-int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
+- double *x, double *y, double *z )
+-
+-{
+- long i;
+
+- srcdefn->ctx->last_errno = 0;
+- dstdefn->ctx->last_errno = 0;
+-
+- if( point_offset == 0 )
+- point_offset = 1;
+
+ /* -------------------------------------------------------------------- */
+ /* Transform unusual input coordinate axis orientation to */
+ /* standard form if needed. */
+ /* -------------------------------------------------------------------- */
+- if( strcmp(srcdefn->axis,"enu") != 0 )
+- {
+- int err;
++static int adjust_axes (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) {
++ /* Nothing to do? */
++ if (0==strcmp(P->axis,"enu"))
++ return 0;
++
++ return pj_adjust_axis( P->ctx, P->axis,
++ dir==PJ_FWD ? 1: 0, n, dist, x, y, z );
++}
+
+- err = pj_adjust_axis( srcdefn->ctx, srcdefn->axis,
+- 0, point_count, point_offset, x, y, z );
+- if( err != 0 )
+- return err;
++
++
++/* ----------------------------------------------------------------------- */
++/* Transform cartesian ("geocentric") source coordinates to lat/long, */
++/* if needed */
++/* ----------------------------------------------------------------------- */
++static int geographic_to_cartesian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) {
++ int res;
++ long i;
++ double fac = P->to_meter;
++
++ /* Nothing to do? */
++ if (!P->is_geocent)
++ return 0;
++
++ if ( z == NULL ) {
++ pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
++ return PJD_ERR_GEOCENTRIC;
++ }
++
++ if (PJ_FWD==dir) {
++ fac = P->fr_meter;
++ res = pj_geodetic_to_geocentric( P->a_orig, P->es_orig, n, dist, x, y, z );
++ if (res)
++ return res;
++ }
++
++ if (fac != 1.0) {
++ for( i = 0; i < n; i++ ) {
++ if( x[dist*i] != HUGE_VAL ) {
++ x[dist*i] *= fac;
++ y[dist*i] *= fac;
++ z[dist*i] *= fac;
++ }
++ }
+ }
+
++ if (PJ_FWD==dir)
++ return 0;
++ return pj_geocentric_to_geodetic(
++ P->a_orig, P->es_orig,
++ n, dist,
++ x, y, z
++ );
++}
++
++
++
++
++
++
++
++
++
++
+ /* -------------------------------------------------------------------- */
+-/* Transform geocentric source coordinates to lat/long. */
++/* Transform destination points to projection coordinates, if */
++/* desired. */
++/* */
++/* Ought to fold this into projected_to_geographic */
+ /* -------------------------------------------------------------------- */
+- if( srcdefn->is_geocent )
++static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *y, double *z) {
++ long i;
++
++ /* Nothing to do? */
++ if (P->is_latlong)
++ return 0;
++ if (P->is_geocent)
++ return 0;
++
++ if(P->fwd3d != NULL)
+ {
+- int err;
+- if( z == NULL )
++ /* Three dimensions must be defined */
++ if ( z == NULL)
+ {
+- pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
++ pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
+ return PJD_ERR_GEOCENTRIC;
+ }
+
+- if( srcdefn->to_meter != 1.0 )
++ for( i = 0; i < n; i++ )
+ {
+- for( i = 0; i < point_count; i++ )
++ XYZ projected_loc;
++ LPZ geodetic_loc;
++
++ geodetic_loc.u = x[dist*i];
++ geodetic_loc.v = y[dist*i];
++ geodetic_loc.w = z[dist*i];
++
++ if (geodetic_loc.u == HUGE_VAL)
++ continue;
++
++ projected_loc = pj_fwd3d( geodetic_loc, P);
++ if( P->ctx->last_errno != 0 )
+ {
+- if( x[point_offset*i] != HUGE_VAL )
++ if( (P->ctx->last_errno != 33 /*EDOM*/
++ && P->ctx->last_errno != 34 /*ERANGE*/ )
++ && (P->ctx->last_errno > 0
++ || P->ctx->last_errno < -44 || n == 1
++ || transient_error[-P->ctx->last_errno] == 0 ) )
++ return P->ctx->last_errno;
++ else
+ {
+- x[point_offset*i] *= srcdefn->to_meter;
+- y[point_offset*i] *= srcdefn->to_meter;
+- z[point_offset*i] *= srcdefn->to_meter;
++ projected_loc.u = HUGE_VAL;
++ projected_loc.v = HUGE_VAL;
++ projected_loc.w = HUGE_VAL;
+ }
+ }
+- }
+
+- err = pj_geocentric_to_geodetic( srcdefn->a_orig, srcdefn->es_orig,
+- point_count, point_offset,
+- x, y, z );
+- if( err != 0 )
+- return err;
++ x[dist*i] = projected_loc.u;
++ y[dist*i] = projected_loc.v;
++ z[dist*i] = projected_loc.w;
++ }
++ return 0;
+ }
+
+-/* -------------------------------------------------------------------- */
+-/* Transform source points to lat/long, if they aren't */
+-/* already. */
+-/* -------------------------------------------------------------------- */
+- else if( !srcdefn->is_latlong )
++ for( i = 0; i <n; i++ )
+ {
++ XY projected_loc;
++ LP geodetic_loc;
+
+- /* Check first if projection is invertible. */
+- if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL))
+- {
+- pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
+- pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR,
+- "pj_transform(): source projection not invertable" );
+- return -17;
+- }
++ geodetic_loc.u = x[dist*i];
++ geodetic_loc.v = y[dist*i];
++
++ if( geodetic_loc.u == HUGE_VAL )
++ continue;
+
+- /* If invertible - First try inv3d if defined */
+- if (srcdefn->inv3d != NULL)
++ projected_loc = pj_fwd( geodetic_loc, P );
++ if( P->ctx->last_errno != 0 )
+ {
+- /* Three dimensions must be defined */
+- if ( z == NULL)
++ if( (P->ctx->last_errno != 33 /*EDOM*/
++ && P->ctx->last_errno != 34 /*ERANGE*/ )
++ && (P->ctx->last_errno > 0
++ || P->ctx->last_errno < -44 || n == 1
++ || transient_error[-P->ctx->last_errno] == 0 ) )
++ return P->ctx->last_errno;
++ else
+ {
+- pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
+- return PJD_ERR_GEOCENTRIC;
++ projected_loc.u = HUGE_VAL;
++ projected_loc.v = HUGE_VAL;
+ }
++ }
+
+- for (i=0; i < point_count; i++)
+- {
+- XYZ projected_loc;
+- XYZ geodetic_loc;
++ x[dist*i] = projected_loc.u;
++ y[dist*i] = projected_loc.v;
++ }
++ return 0;
++}
+
+- projected_loc.u = x[point_offset*i];
+- projected_loc.v = y[point_offset*i];
+- projected_loc.w = z[point_offset*i];
+
+- if (projected_loc.u == HUGE_VAL)
+- continue;
+
+- geodetic_loc = pj_inv3d(projected_loc, srcdefn);
+- if( srcdefn->ctx->last_errno != 0 )
+- {
+- if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
+- && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
+- && (srcdefn->ctx->last_errno > 0
+- || srcdefn->ctx->last_errno < -44 || point_count == 1
+- || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
+- return srcdefn->ctx->last_errno;
+- else
+- {
+- geodetic_loc.u = HUGE_VAL;
+- geodetic_loc.v = HUGE_VAL;
+- geodetic_loc.w = HUGE_VAL;
+- }
+- }
+
+- x[point_offset*i] = geodetic_loc.u;
+- y[point_offset*i] = geodetic_loc.v;
+- z[point_offset*i] = geodetic_loc.w;
+
+- }
++/* ----------------------------------------------------------------------- */
++/* Transform projected source coordinates to lat/long, if needed */
++/* ----------------------------------------------------------------------- */
++static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *y, double *z) {
++ long i;
+
++ /* Nothing to do? */
++ if (P->is_latlong)
++ return 0;
++
++ /* Check first if projection is invertible. */
++ if( (P->inv3d == NULL) && (P->inv == NULL))
++ {
++ pj_ctx_set_errno( pj_get_ctx(P), -17 );
++ pj_log( pj_get_ctx(P), PJ_LOG_ERROR,
++ "pj_transform(): source projection not invertable" );
++ return -17;
++ }
++
++ /* If invertible - First try inv3d if defined */
++ if (P->inv3d != NULL)
++ {
++ /* Three dimensions must be defined */
++ if ( z == NULL)
++ {
++ pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC);
++ return PJD_ERR_GEOCENTRIC;
+ }
+- else
++
++ for (i=0; i < n; i++)
+ {
+- /* Fallback to the original PROJ.4 API 2d inversion - inv */
+- for( i = 0; i < point_count; i++ )
+- {
+- XY projected_loc;
+- LP geodetic_loc;
++ XYZ projected_loc;
++ XYZ geodetic_loc;
+
+- projected_loc.u = x[point_offset*i];
+- projected_loc.v = y[point_offset*i];
++ projected_loc.u = x[dist*i];
++ projected_loc.v = y[dist*i];
++ projected_loc.w = z[dist*i];
+
+- if( projected_loc.u == HUGE_VAL )
+- continue;
++ if (projected_loc.u == HUGE_VAL)
++ continue;
+
+- geodetic_loc = pj_inv( projected_loc, srcdefn );
+- if( srcdefn->ctx->last_errno != 0 )
++ geodetic_loc = pj_inv3d(projected_loc, P);
++ if( P->ctx->last_errno != 0 )
++ {
++ if( (P->ctx->last_errno != 33 /*EDOM*/
++ && P->ctx->last_errno != 34 /*ERANGE*/ )
++ && (P->ctx->last_errno > 0
++ || P->ctx->last_errno < -44 || n == 1
++ || transient_error[-P->ctx->last_errno] == 0 ) )
++ return P->ctx->last_errno;
++ else
+ {
+- if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
+- && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
+- && (srcdefn->ctx->last_errno > 0
+- || srcdefn->ctx->last_errno < -44 || point_count == 1
+- || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
+- return srcdefn->ctx->last_errno;
+- else
+- {
+- geodetic_loc.u = HUGE_VAL;
+- geodetic_loc.v = HUGE_VAL;
+- }
++ geodetic_loc.u = HUGE_VAL;
++ geodetic_loc.v = HUGE_VAL;
++ geodetic_loc.w = HUGE_VAL;
+ }
+-
+- x[point_offset*i] = geodetic_loc.u;
+- y[point_offset*i] = geodetic_loc.v;
+ }
+- }
+- }
+
+-/* -------------------------------------------------------------------- */
+-/* But if they are already lat long, adjust for the prime */
+-/* meridian if there is one in effect. */
+-/* -------------------------------------------------------------------- */
+- if ((srcdefn->is_geocent || srcdefn->is_latlong) && ( srcdefn->from_greenwich != 0.0 ))
+- {
+- for( i = 0; i < point_count; i++ )
+- {
+- if( x[point_offset*i] != HUGE_VAL )
+- x[point_offset*i] += srcdefn->from_greenwich;
++ x[dist*i] = geodetic_loc.u;
++ y[dist*i] = geodetic_loc.v;
++ z[dist*i] = geodetic_loc.w;
++
+ }
++ return 0;
+ }
+
+-/* -------------------------------------------------------------------- */
+-/* Do we need to translate from geoid to ellipsoidal vertical */
+-/* datum? */
+-/* -------------------------------------------------------------------- */
+- if( srcdefn->has_geoid_vgrids && z != NULL )
+- {
+- if( pj_apply_vgridshift( srcdefn, "sgeoidgrids",
+- &(srcdefn->vgridlist_geoid),
+- &(srcdefn->vgridlist_geoid_count),
+- 0, point_count, point_offset, x, y, z ) != 0 )
+- return pj_ctx_get_errno(srcdefn->ctx);
+- }
++ /* Fallback to the original PROJ.4 API 2d inversion - inv */
++ for( i = 0; i < n; i++ ) {
++ XY projected_loc;
++ LP geodetic_loc;
+
+-/* -------------------------------------------------------------------- */
+-/* Convert datums if needed, and possible. */
+-/* -------------------------------------------------------------------- */
+- if( pj_datum_transform( srcdefn, dstdefn, point_count, point_offset,
+- x, y, z ) != 0 )
+- {
+- if( srcdefn->ctx->last_errno != 0 )
+- return srcdefn->ctx->last_errno;
+- else
+- return dstdefn->ctx->last_errno;
+- }
++ projected_loc.u = x[dist*i];
++ projected_loc.v = y[dist*i];
+
+-/* -------------------------------------------------------------------- */
+-/* Do we need to translate from ellipsoidal to geoid vertical */
+-/* datum? */
+-/* -------------------------------------------------------------------- */
+- if( dstdefn->has_geoid_vgrids && z != NULL )
+- {
+- if( pj_apply_vgridshift( dstdefn, "sgeoidgrids",
+- &(dstdefn->vgridlist_geoid),
+- &(dstdefn->vgridlist_geoid_count),
+- 1, point_count, point_offset, x, y, z ) != 0 )
+- return dstdefn->ctx->last_errno;
+- }
++ if( projected_loc.u == HUGE_VAL )
++ continue;
+
+-/* -------------------------------------------------------------------- */
+-/* But if they are staying lat long, adjust for the prime */
+-/* meridian if there is one in effect. */
+-/* -------------------------------------------------------------------- */
+- if ((dstdefn->is_geocent || dstdefn->is_latlong) && ( dstdefn->from_greenwich != 0.0 ))
+- {
+- for( i = 0; i < point_count; i++ )
++ geodetic_loc = pj_inv( projected_loc, P );
++ if( P->ctx->last_errno != 0 )
+ {
+- if( x[point_offset*i] != HUGE_VAL )
+- x[point_offset*i] -= dstdefn->from_greenwich;
++ if( (P->ctx->last_errno != 33 /*EDOM*/
++ && P->ctx->last_errno != 34 /*ERANGE*/ )
++ && (P->ctx->last_errno > 0
++ || P->ctx->last_errno < -44 || n == 1
++ || transient_error[-P->ctx->last_errno] == 0 ) )
++ return P->ctx->last_errno;
++ else
++ {
++ geodetic_loc.u = HUGE_VAL;
++ geodetic_loc.v = HUGE_VAL;
++ }
+ }
++
++ x[dist*i] = geodetic_loc.u;
++ y[dist*i] = geodetic_loc.v;
++ }
++ return 0;
+ }
+
++
++
+ /* -------------------------------------------------------------------- */
+-/* Transform destination latlong to geocentric if required. */
++/* Adjust for the prime meridian if needed. */
+ /* -------------------------------------------------------------------- */
+- if( dstdefn->is_geocent )
+- {
+- if( z == NULL )
+- {
+- pj_ctx_set_errno( dstdefn->ctx, PJD_ERR_GEOCENTRIC );
+- return PJD_ERR_GEOCENTRIC;
+- }
++static int prime_meridian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x) {
++ int i;
++ double pm = P->from_greenwich;
++
++ /* Nothing to do? */
++ if (pm==0.0)
++ return 0;
++ if (!(P->is_geocent || P->is_latlong))
++ return 0;
++
++ if (dir==PJ_FWD)
++ pm = -pm;
++
++ for (i = 0; i < n; i++)
++ if (x[dist*i] != HUGE_VAL)
++ x[dist*i] += pm;
++
++ return 0;
++}
+
+- pj_geodetic_to_geocentric( dstdefn->a_orig, dstdefn->es_orig,
+- point_count, point_offset, x, y, z );
+
+- if( dstdefn->fr_meter != 1.0 )
+- {
+- for( i = 0; i < point_count; i++ )
+- {
+- if( x[point_offset*i] != HUGE_VAL )
+- {
+- x[point_offset*i] *= dstdefn->fr_meter;
+- y[point_offset*i] *= dstdefn->fr_meter;
+- z[point_offset*i] *= srcdefn->fr_meter;
+- }
+- }
+- }
+- }
+
+ /* -------------------------------------------------------------------- */
+-/* Transform destination points to projection coordinates, if */
+-/* desired. */
++/* Adjust for vertical scale factor if needed */
+ /* -------------------------------------------------------------------- */
+- else if( !dstdefn->is_latlong )
+- {
++static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) {
++ int i;
++ double fac = P->vto_meter;
++
++ if (PJ_FWD==dir)
++ fac = P->vfr_meter;
++
++ /* Nothing to do? */
++ if (P->vto_meter==0.0)
++ return 0;
++
++ for (i = 0; i < n; i++)
++ if (z[dist*i] != HUGE_VAL )
++ z[dist*i] *= fac;
+
+- if( dstdefn->fwd3d != NULL)
+- {
+- /* Three dimensions must be defined */
+- if ( z == NULL)
+- {
+- pj_ctx_set_errno( pj_get_ctx(dstdefn), PJD_ERR_GEOCENTRIC);
+- return PJD_ERR_GEOCENTRIC;
+- }
++ return 0;
++}
+
+- for( i = 0; i < point_count; i++ )
+- {
+- XYZ projected_loc;
+- LPZ geodetic_loc;
+
+- geodetic_loc.u = x[point_offset*i];
+- geodetic_loc.v = y[point_offset*i];
+- geodetic_loc.w = z[point_offset*i];
+
+- if (geodetic_loc.u == HUGE_VAL)
+- continue;
++/* -------------------------------------------------------------------- */
++/* Transform to ellipsoidal heights if needed */
++/* -------------------------------------------------------------------- */
++static int geometric_to_orthometric (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) {
++ int err;
++ if (0==P->has_geoid_vgrids)
++ return 0;
++ if (z==0)
++ return PJD_ERR_GEOCENTRIC;
++ err = pj_apply_vgridshift (P, "sgeoidgrids",
++ &(P->vgridlist_geoid),
++ &(P->vgridlist_geoid_count),
++ dir==PJ_FWD ? 1 : 0, n, dist, x, y, z );
++ if (err)
++ return pj_ctx_get_errno(P->ctx);
++ return 0;
++}
+
+- projected_loc = pj_fwd3d( geodetic_loc, dstdefn);
+- if( dstdefn->ctx->last_errno != 0 )
+- {
+- if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
+- && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
+- && (dstdefn->ctx->last_errno > 0
+- || dstdefn->ctx->last_errno < -44 || point_count == 1
+- || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
+- return dstdefn->ctx->last_errno;
+- else
+- {
+- projected_loc.u = HUGE_VAL;
+- projected_loc.v = HUGE_VAL;
+- projected_loc.w = HUGE_VAL;
+- }
+- }
+
+- x[point_offset*i] = projected_loc.u;
+- y[point_offset*i] = projected_loc.v;
+- z[point_offset*i] = projected_loc.w;
+- }
+
+- }
+- else
+- {
+- for( i = 0; i < point_count; i++ )
+- {
+- XY projected_loc;
+- LP geodetic_loc;
++/* -------------------------------------------------------------------- */
++/* Convert datums if needed, and possible. */
++/* -------------------------------------------------------------------- */
++static int datum_transform (PJ *P, PJ *Q, long n, int dist, double *x, double *y, double *z) {
++ if (0==pj_datum_transform (P, Q, n, dist, x, y, z))
++ return 0;
++ if (P->ctx->last_errno)
++ return P->ctx->last_errno;
++ return Q->ctx->last_errno;
++}
+
+- geodetic_loc.u = x[point_offset*i];
+- geodetic_loc.v = y[point_offset*i];
+
+- if( geodetic_loc.u == HUGE_VAL )
+- continue;
+
+- projected_loc = pj_fwd( geodetic_loc, dstdefn );
+- if( dstdefn->ctx->last_errno != 0 )
+- {
+- if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
+- && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
+- && (dstdefn->ctx->last_errno > 0
+- || dstdefn->ctx->last_errno < -44 || point_count == 1
+- || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
+- return dstdefn->ctx->last_errno;
+- else
+- {
+- projected_loc.u = HUGE_VAL;
+- projected_loc.v = HUGE_VAL;
+- }
+- }
+
+- x[point_offset*i] = projected_loc.u;
+- y[point_offset*i] = projected_loc.v;
+- }
+- }
+- }
+
+ /* -------------------------------------------------------------------- */
+ /* If a wrapping center other than 0 is provided, rewrap around */
+ /* the suggested center (for latlong coordinate systems only). */
+ /* -------------------------------------------------------------------- */
+- else if( dstdefn->is_latlong && dstdefn->is_long_wrap_set )
+- {
+- for( i = 0; i < point_count; i++ )
+- {
+- double val = x[point_offset*i];
+- if( val == HUGE_VAL )
+- continue;
++static int long_wrap (PJ *P, long n, int dist, double *x) {
++ long i;
+
+- /* Get fast in ] -2 PI, 2 PI [ range */
+- val = fmod(val, M_TWOPI);
+- while( val < dstdefn->long_wrap_center - M_PI )
+- val += M_TWOPI;
+- while( val > dstdefn->long_wrap_center + M_PI )
+- val -= M_TWOPI;
+- x[point_offset*i] = val;
+- }
++ /* Nothing to do? */
++ if (P->is_geocent)
++ return 0;
++ if (!P->is_long_wrap_set)
++ return 0;
++ if (!P->is_latlong)
++ return 0;
++
++ for (i = 0; i < n; i++ ) {
++ double val = x[dist*i];
++ if (val == HUGE_VAL)
++ continue;
++
++ /* Get fast in ] -2 PI, 2 PI [ range */
++ val = fmod(val, M_TWOPI);
++ while( val < P->long_wrap_center - M_PI )
++ val += M_TWOPI;
++ while( val > P->long_wrap_center + M_PI )
++ val -= M_TWOPI;
++ x[dist*i] = val;
+ }
++ return 0;
++}
+
+-/* -------------------------------------------------------------------- */
+-/* Transform normalized axes into unusual output coordinate axis */
+-/* orientation if needed. */
+-/* -------------------------------------------------------------------- */
+- if( strcmp(dstdefn->axis,"enu") != 0 )
+- {
+- int err;
+
+- err = pj_adjust_axis( dstdefn->ctx, dstdefn->axis,
+- 1, point_count, point_offset, x, y, z );
+- if( err != 0 )
+- return err;
+- }
++
++/************************************************************************/
++/* pj_transform() */
++/* */
++/* Currently this function doesn't recognise if two projections */
++/* are identical (to short circuit reprojection) because it is */
++/* difficult to compare PJ structures (since there are some */
++/* projection specific components). */
++/************************************************************************/
++
++int pj_transform(
++ PJ *srcdefn, PJ *dstdefn,
++ long point_count, int point_offset,
++ double *x, double *y, double *z
++){
++ int err;
++
++ srcdefn->ctx->last_errno = 0;
++ dstdefn->ctx->last_errno = 0;
++
++ if( point_offset == 0 )
++ point_offset = 1;
++
++ /* Bring input to "normal form": longitude, latitude, ellipsoidal height */
++
++ err = adjust_axes (srcdefn, PJ_INV, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++ err = geographic_to_cartesian (srcdefn, PJ_INV, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++ err = projected_to_geographic (srcdefn, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++ err = prime_meridian (srcdefn, PJ_INV, point_count, point_offset, x);
++ if (err)
++ return err;
++ err = height_unit (srcdefn, PJ_INV, point_count, point_offset, z);
++ if (err)
++ return err;
++ err = geometric_to_orthometric (srcdefn, PJ_INV, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++
++ /* At the center of the process we do the datum shift (if needed) */
++
++ err = datum_transform(srcdefn, dstdefn, point_count, point_offset, x, y, z );
++ if (err)
++ return err;
++
++ /* Now get out on the other side: Bring "normal form" to output form */
++
++ err = geometric_to_orthometric (dstdefn, PJ_FWD, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++ err = height_unit (dstdefn, PJ_FWD, point_count, point_offset, z);
++ if (err)
++ return err;
++ err = prime_meridian (dstdefn, PJ_FWD, point_count, point_offset, x);
++ if (err)
++ return err;
++ err = geographic_to_cartesian (dstdefn, PJ_FWD, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++ err = geographic_to_projected (dstdefn, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
++ err = long_wrap (dstdefn, point_count, point_offset, x);
++ if (err)
++ return err;
++ err = adjust_axes (dstdefn, PJ_FWD, point_count, point_offset, x, y, z);
++ if (err)
++ return err;
+
+ return 0;
+ }
+
++
++
+ /************************************************************************/
+ /* pj_geodetic_to_geocentric() */
+ /************************************************************************/
diff --git a/debian/patches/series b/debian/patches/series
new file mode 100644
index 0000000..acc3907
--- /dev/null
+++ b/debian/patches/series
@@ -0,0 +1,2 @@
+0001-Revert-fix-to-22.patch
+pr845_Refactor-pj_transform-reintroduce-support-for-vertical-scaling.patch
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/proj.git
More information about the Pkg-grass-devel
mailing list