From 1583a566a208d2451fb1acc8bcf16fbd8151983e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 29 Sep 2018 10:01:48 +0200 Subject: Add geographic offset transformation method. The Geographic offsets transformation adds an offset to the geographic longitude, latitude coordinates, and an offset to the ellipsoidal height. This method is normally only used when low accuracy is tolerated. It is documented as coordinate operation method code 9619 (for geographic 2D) and 9660 (for geographic 3D) in the EPSG dataset. It can also be used to implement the method Geographic2D with Height Offsets (code 9618) by noting that the input vertical component is a gravity-related height and the output vertical component is the ellispoid height (dh being the geoid undulation). It can also be used to implement the method Vertical offset (code 9616) It is used for example to transform: - from the old Greek geographic 2D CRS to the newer GGRS87 CRS - from Tokyo + JSLD69 height to WGS 84 - from Baltic 1977 height to Black Sea height It is also useful to document the implicit zero-offset transformation we do in pipelines such as +proj=pipeline +step +inv +proj=longlat +ellps=A +step +proj=longlat +ellps=B that can be explicited as +proj=pipeline +step +inv +proj=longlat +ellps=A +step +proj=geogoffset [+dlon=0 +dlat=0 +dh=0] +step +proj=longlat +ellps=B --- .../operations/transformations/geogoffset.rst | 68 +++++++++++++ docs/source/operations/transformations/index.rst | 1 + docs/source/references.bib | 8 ++ src/Makefile.am | 2 +- src/PJ_geogoffset.c | 111 +++++++++++++++++++++ src/lib_proj.cmake | 1 + src/makefile.vc | 2 +- src/pj_list.h | 1 + test/gie/more_builtins.gie | 32 ++++++ 9 files changed, 224 insertions(+), 2 deletions(-) create mode 100644 docs/source/operations/transformations/geogoffset.rst create mode 100644 src/PJ_geogoffset.c diff --git a/docs/source/operations/transformations/geogoffset.rst b/docs/source/operations/transformations/geogoffset.rst new file mode 100644 index 00000000..91c91c83 --- /dev/null +++ b/docs/source/operations/transformations/geogoffset.rst @@ -0,0 +1,68 @@ +.. _geogoffset: + +================================================================================ +Geographic offsets +================================================================================ + +.. versionadded:: 6.0.0 + +The Geographic offsets transformation adds an offset to the geographic longitude, +latitude coordinates, and an offset to the ellipsoidal height. + ++---------------------+----------------------------------------------------------+ +| **Alias** | geogoffset | ++---------------------+----------------------------------------------------------+ +| **Domain** | 3D | ++---------------------+----------------------------------------------------------+ +| **Input type** | Geodetic coordinates (horizontal), meters (vertical) | ++---------------------+----------------------------------------------------------+ +| **output type** | Geodetic coordinates (horizontal), meters (vertical) | ++---------------------+----------------------------------------------------------+ + +This method is normally only used when low accuracy is tolerated. It is documented +as coordinate operation method code 9619 (for geographic 2D) and 9660 (for +geographic 3D) in the EPSG dataset (:cite:`EPSGGuidanceNumber7Part2`) + +It can also be used to implement the method Geographic2D with Height Offsets +(code 9618) by noting that the input vertical component is a gravity-related +height and the output vertical component is the ellispoid height (dh being +the geoid undulation). + +It can also be used to implement the method Vertical offset (code 9616) + +The reverse transformation simply consists in subtracting the offsets. + +Examples +############################################################################### + +Geographic offset from the old Greek geographic 2D CRS to the newer GGRS87 CRS:: + + proj=geogoffset dlon=0.28 dlat=-5.86 + +Conversion from Tokyo + JSLD69 height to WGS 84:: + + proj=geogoffset dlon=-13.97 dlat=7.94 dh=26.9 + +Conversion from Baltic 1977 height to Black Sea height:: + + proj=geogoffset dh=0.4 + + +Parameters +################################################################################ + +Optional +------------------------------------------------------------------------------- + +.. option:: +dlon= + + Offset in longitude, expressed in arc-second, to add. + +.. option:: +dlat= + + Offset in latitude, expressed in arc-second, to add. + +.. option:: +dh= + + Offset in height, expressed in meter, to add. + diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst index 5112ffa2..c2fe6baa 100644 --- a/docs/source/operations/transformations/index.rst +++ b/docs/source/operations/transformations/index.rst @@ -11,6 +11,7 @@ systems are based on different datums. :maxdepth: 1 deformation + geogoffset helmert molodensky hgridshift diff --git a/docs/source/references.bib b/docs/source/references.bib index bd3089ee..0fb077ec 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -82,6 +82,14 @@ Url = {http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf} } +@Manual{EPSGGuidanceNumber7Part2, + Title = {Geomatics Guidance Note 7, part 2 - Coordinate Conversions & Transformations including Formulas}, + Institution = {International Association For Oil And Gas Producers}, + Year = {2018}, + + Url = {https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/} +} + @Manual{Evenden2005, Title = {libproj4: A Comprehensive Library of Cartographic Projection Functions (Preliminary Draft)}, Author = {G. I. Evenden}, diff --git a/src/Makefile.am b/src/Makefile.am index 84846c1c..2320e812 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -89,7 +89,7 @@ libproj_la_SOURCES = \ \ proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \ PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \ - PJ_deformation.c pj_internal.c PJ_axisswap.c + PJ_deformation.c pj_internal.c PJ_axisswap.c PJ_geogoffset.c install-exec-local: rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT) diff --git a/src/PJ_geogoffset.c b/src/PJ_geogoffset.c new file mode 100644 index 00000000..7f917fd7 --- /dev/null +++ b/src/PJ_geogoffset.c @@ -0,0 +1,111 @@ +/************************************************************************ +* Copyright (c) 2018, Even Rouault +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* 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 +#include + +#include "proj.h" +#include "projects.h" + +PROJ_HEAD(geogoffset, "Geographic Offset"); + +struct pj_opaque_geogoffset { + double dlon; + double dlat; + double dh; +}; + + +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { + struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque; + /* offset coordinate */ + obs.lpz.lam += Q->dlon; + obs.lpz.phi += Q->dlat; + obs.lpz.z += Q->dh; + return obs; +} + +static XYZ forward_3d(LPZ lpz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lpz = lpz; + return forward_4d(point, P).xyz; +} + + +static XY forward_2d(LP lp, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lp = lp; + return forward_4d(point, P).xy; +} + + +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { + struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque; + /* offset coordinate */ + obs.lpz.lam -= Q->dlon; + obs.lpz.phi -= Q->dlat; + obs.lpz.z -= Q->dh; + return obs; +} + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xyz = xyz; + return reverse_4d(point, P).lpz; +} + +static LP reverse_2d(XY xy, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xy = xy; + return reverse_4d(point, P).lp; +} + + +/* Arcsecond to radians */ +#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) + + +PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) { + struct pj_opaque_geogoffset *Q = pj_calloc(1, sizeof(struct pj_opaque_geogoffset)); + if (0==Q) + return pj_default_destructor(P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + /* read args */ + Q->dlon = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD; + Q->dlat = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD; + Q->dh = pj_param(P->ctx, P->params, "ddh").f; + + return P; +} diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 0779568a..f1e76ff3 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -72,6 +72,7 @@ SET(SRC_LIBPROJ_PJ PJ_fouc_s.c PJ_gall.c PJ_geoc.c + PJ_geogoffset.c PJ_geos.c PJ_gins8.c PJ_gnom.c diff --git a/src/makefile.vc b/src/makefile.vc index cf9d878c..21db1362 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -61,7 +61,7 @@ support = \ pipeline = \ proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \ - PJ_deformation.obj PJ_axisswap.obj + PJ_deformation.obj PJ_axisswap.obj PJ_geogoffset.obj geodesic = geodesic.obj diff --git a/src/pj_list.h b/src/pj_list.h index dfcfe30a..6d94f4fd 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -46,6 +46,7 @@ PROJ_HEAD(fouc_s, "Foucaut Sinusoidal") PROJ_HEAD(gall, "Gall (Gall Stereographic)") PROJ_HEAD(geoc, "Geocentric Latitude") PROJ_HEAD(geocent, "Geocentric") +PROJ_HEAD(geogoffset, "Geographic Offset") PROJ_HEAD(geos, "Geostationary Satellite View") PROJ_HEAD(gins8, "Ginsburg VIII (TsNIIGAiK)") PROJ_HEAD(gn_sinu, "General Sinusoidal Series") diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index 96212cba..d3fe64f6 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -557,6 +557,38 @@ expect 0 -90 ------------------------------------------------------------------------------- + +------------------------------------------------------------------------------- +Test for PJ_geogoffset +------------------------------------------------------------------------------- +------------------------------------------------------------------------------- +operation +proj=geogoffset +------------------------------------------------------------------------------- +direction forward +tolerance 1mm + +accept 10 20 +expect 10 20 +roundtrip 1 + +------------------------------------------------------------------------------- +operation +proj=geogoffset +dlon=3600 +dlat=-3600 +dh=3 +------------------------------------------------------------------------------- +direction forward +tolerance 1mm + +accept 10 20 +expect 11 19 +roundtrip 1 + +accept 10 20 30 +expect 11 19 33 +roundtrip 1 + +accept 10 20 30 40 +expect 11 19 33 40 +roundtrip 1 + ------------------------------------------------------------------------------- run the few gie-builtin tests, which are currently either awkward or impossible to express in the gie command set -- cgit v1.2.3 From fe29f8acdce81607c11a597f4bffc7ff61fa9c19 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 1 Oct 2018 12:32:52 +0200 Subject: Add a affine transformation method, and make geogoffset as a particular case of it (fixes #535) --- docs/source/operations/transformations/affine.rst | 123 +++++++++++ .../operations/transformations/geogoffset.rst | 2 + docs/source/operations/transformations/index.rst | 1 + src/Makefile.am | 2 +- src/PJ_affine.c | 246 +++++++++++++++++++++ src/PJ_geogoffset.c | 111 ---------- src/lib_proj.cmake | 2 +- src/makefile.vc | 2 +- src/pj_list.h | 1 + test/gie/more_builtins.gie | 46 +++- 10 files changed, 421 insertions(+), 115 deletions(-) create mode 100644 docs/source/operations/transformations/affine.rst create mode 100644 src/PJ_affine.c delete mode 100644 src/PJ_geogoffset.c diff --git a/docs/source/operations/transformations/affine.rst b/docs/source/operations/transformations/affine.rst new file mode 100644 index 00000000..8b1681b3 --- /dev/null +++ b/docs/source/operations/transformations/affine.rst @@ -0,0 +1,123 @@ +.. _affine: + +================================================================================ +Affine transformation +================================================================================ + +.. versionadded:: 6.0.0 + +The affine transformation applies translation and scaling/rotation terms on the +x,y,z coordinates, and translation and scaling on the temporal cordinate. + ++---------------------+----------------------------------------------------------+ +| **Alias** | affine | ++---------------------+----------------------------------------------------------+ +| **Domain** | 4D | ++---------------------+----------------------------------------------------------+ +| **Input type** | XYZT | ++---------------------+----------------------------------------------------------+ +| **output type** | XYZT | ++---------------------+----------------------------------------------------------+ + +By default, the parameters are set for an identity transforms. The transformation +is reversible unless the determinant of the sji matrix is 0, or `tscale` is 0 + + +Parameters +################################################################################ + +Optional +------------------------------------------------------------------------------- + +.. option:: +xoff= + + Offset in X. Default value: 0 + +.. option:: +yoff= + + Offset in Y. Default value: 0 + +.. option:: +zoff= + + Offset in Z. Default value: 0 + +.. option:: +toff= + + Offset in T. Default value: 0 + +.. option:: +s11= + + Rotation/scaling term. Default value: 1 + +.. option:: +s12= + + Rotation/scaling term. Default value: 0 + +.. option:: +s13= + + Rotation/scaling term. Default value: 0 + +.. option:: +s21= + + Rotation/scaling term. Default value: 0 + +.. option:: +s22= + + Rotation/scaling term. Default value: 1 + +.. option:: +s23= + + Rotation/scaling term. Default value: 0 + +.. option:: +s31= + + Rotation/scaling term. Default value: 0 + +.. option:: +s32= + + Rotation/scaling term. Default value: 0 + +.. option:: +s33= + + Rotation/scaling term. Default value: 1 + +.. option:: +tscale= + + Time scaling term. Default value: 1 + + + +Mathematical description ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. math:: + :label: formula + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + T \\ + \end{bmatrix}^{dest} = + \begin{bmatrix} + xoff \\ + yoff \\ + zoff \\ + toff \\ + \end{bmatrix} + + \begin{bmatrix} + s11 & s12 & s13 & 0 \\ + s21 & s22 & s23 & 0 \\ + s31 & s32 & s33 & 0 \\ + 0 & 0 & 0 & tscale \\ + \end{bmatrix} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + T \\ + \end{bmatrix}^{source} + \end{align} + + diff --git a/docs/source/operations/transformations/geogoffset.rst b/docs/source/operations/transformations/geogoffset.rst index 91c91c83..f643485e 100644 --- a/docs/source/operations/transformations/geogoffset.rst +++ b/docs/source/operations/transformations/geogoffset.rst @@ -32,6 +32,8 @@ It can also be used to implement the method Vertical offset (code 9616) The reverse transformation simply consists in subtracting the offsets. +This method is a conveniency wrapper for the more general :ref:`affine`. + Examples ############################################################################### diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst index c2fe6baa..121215f7 100644 --- a/docs/source/operations/transformations/index.rst +++ b/docs/source/operations/transformations/index.rst @@ -10,6 +10,7 @@ systems are based on different datums. .. toctree:: :maxdepth: 1 + affine deformation geogoffset helmert diff --git a/src/Makefile.am b/src/Makefile.am index 2320e812..c1fbc12e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -89,7 +89,7 @@ libproj_la_SOURCES = \ \ proj_4D_api.c PJ_cart.c PJ_pipeline.c PJ_horner.c PJ_helmert.c \ PJ_vgridshift.c PJ_hgridshift.c PJ_unitconvert.c PJ_molodensky.c \ - PJ_deformation.c pj_internal.c PJ_axisswap.c PJ_geogoffset.c + PJ_deformation.c pj_internal.c PJ_axisswap.c PJ_affine.c install-exec-local: rm -f $(DESTDIR)$(bindir)/invproj$(EXEEXT) diff --git a/src/PJ_affine.c b/src/PJ_affine.c new file mode 100644 index 00000000..0d8b6374 --- /dev/null +++ b/src/PJ_affine.c @@ -0,0 +1,246 @@ +/************************************************************************ +* Copyright (c) 2018, Even Rouault +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* 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 +#include + +#include "proj_internal.h" +#include "proj.h" +#include "projects.h" + +PROJ_HEAD(affine, "Affine transformation"); +PROJ_HEAD(geogoffset, "Geographic Offset"); + +struct pj_affine_coeffs { + double s11; + double s12; + double s13; + double s21; + double s22; + double s23; + double s31; + double s32; + double s33; + double tscale; +}; + +struct pj_opaque_affine { + double xoff; + double yoff; + double zoff; + double toff; + struct pj_affine_coeffs forward; + struct pj_affine_coeffs reverse; +}; + + +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { + PJ_COORD newObs; + const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque; + const struct pj_affine_coeffs *C = &(Q->forward); + newObs.xyzt.x = Q->xoff + C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z; + newObs.xyzt.y = Q->yoff + C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z; + newObs.xyzt.z = Q->zoff + C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z; + newObs.xyzt.t = Q->toff + C->tscale * obs.xyzt.t; + return newObs; +} + +static XYZ forward_3d(LPZ lpz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lpz = lpz; + return forward_4d(point, P).xyz; +} + + +static XY forward_2d(LP lp, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.lp = lp; + return forward_4d(point, P).xy; +} + + +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { + PJ_COORD newObs; + const struct pj_opaque_affine *Q = (const struct pj_opaque_affine *) P->opaque; + const struct pj_affine_coeffs *C = &(Q->reverse); + obs.xyzt.x -= Q->xoff; + obs.xyzt.y -= Q->yoff; + obs.xyzt.z -= Q->zoff; + newObs.xyzt.x = C->s11 * obs.xyzt.x + C->s12 * obs.xyzt.y + C->s13 * obs.xyzt.z; + newObs.xyzt.y = C->s21 * obs.xyzt.x + C->s22 * obs.xyzt.y + C->s23 * obs.xyzt.z; + newObs.xyzt.z = C->s31 * obs.xyzt.x + C->s32 * obs.xyzt.y + C->s33 * obs.xyzt.z; + newObs.xyzt.t = C->tscale * (obs.xyzt.t - Q->toff); + return newObs; +} + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xyz = xyz; + return reverse_4d(point, P).lpz; +} + +static LP reverse_2d(XY xy, PJ *P) { + PJ_COORD point = {{0,0,0,0}}; + point.xy = xy; + return reverse_4d(point, P).lp; +} + +static struct pj_opaque_affine * initQ() { + struct pj_opaque_affine *Q = pj_calloc(1, sizeof(struct pj_opaque_affine)); + if (0==Q) + return 0; + + /* default values */ + Q->forward.s11 = 1.0; + Q->forward.s22 = 1.0; + Q->forward.s33 = 1.0; + Q->forward.tscale = 1.0; + + Q->reverse.s11 = 1.0; + Q->reverse.s22 = 1.0; + Q->reverse.s33 = 1.0; + Q->reverse.tscale = 1.0; + + return Q; +} + +static void computeReverseParameters(PJ* P) +{ + struct pj_opaque_affine *Q = (struct pj_opaque_affine *) P->opaque; + + /* cf https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices */ + const double a = Q->forward.s11; + const double b = Q->forward.s12; + const double c = Q->forward.s13; + const double d = Q->forward.s21; + const double e = Q->forward.s22; + const double f = Q->forward.s23; + const double g = Q->forward.s31; + const double h = Q->forward.s32; + const double i = Q->forward.s33; + const double A = e * i - f * h; + const double B = -(d * i - f * g); + const double C = (d * h - e * g); + const double D = -(b * i - c * h); + const double E = (a * i - c * g); + const double F = -(a * h - b * g); + const double G = b * f - c * e; + const double H = -(a * f - c * d); + const double I = a * e - b * d; + const double det = a * A + b * B + c * C; + if( det == 0.0 || Q->forward.tscale == 0.0 ) { + if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) { + proj_log_debug(P, "Affine: matrix non invertible"); + } + P->inv4d = NULL; + P->inv3d = NULL; + P->inv = NULL; + } else { + Q->reverse.s11 = A / det; + Q->reverse.s12 = D / det; + Q->reverse.s13 = G / det; + Q->reverse.s21 = B / det; + Q->reverse.s22 = E / det; + Q->reverse.s23 = H / det; + Q->reverse.s31 = C / det; + Q->reverse.s32 = F / det; + Q->reverse.s33 = I / det; + Q->reverse.tscale = 1.0 / Q->forward.tscale; + } +} + +PJ *TRANSFORMATION(affine,0 /* no need for ellipsoid */) { + struct pj_opaque_affine *Q = initQ(); + if (0==Q) + return pj_default_destructor(P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_WHATEVER; + P->right = PJ_IO_UNITS_WHATEVER; + + /* read args */ + Q->xoff = pj_param(P->ctx, P->params, "dxoff").f; + Q->yoff = pj_param(P->ctx, P->params, "dyoff").f; + Q->zoff = pj_param(P->ctx, P->params, "dzoff").f; + Q->toff = pj_param(P->ctx, P->params, "dtoff").f; + + if(pj_param (P->ctx, P->params, "ts11").i) { + Q->forward.s11 = pj_param(P->ctx, P->params, "ds11").f; + } + Q->forward.s12 = pj_param(P->ctx, P->params, "ds12").f; + Q->forward.s13 = pj_param(P->ctx, P->params, "ds13").f; + Q->forward.s21 = pj_param(P->ctx, P->params, "ds21").f; + if(pj_param (P->ctx, P->params, "ts22").i) { + Q->forward.s22 = pj_param(P->ctx, P->params, "ds22").f; + } + Q->forward.s23 = pj_param(P->ctx, P->params, "ds23").f; + Q->forward.s31 = pj_param(P->ctx, P->params, "ds31").f; + Q->forward.s32 = pj_param(P->ctx, P->params, "ds32").f; + if(pj_param (P->ctx, P->params, "ts33").i) { + Q->forward.s33 = pj_param(P->ctx, P->params, "ds33").f; + } + if(pj_param (P->ctx, P->params, "ttscale").i) { + Q->forward.tscale = pj_param(P->ctx, P->params, "dtscale").f; + } + + computeReverseParameters(P); + + return P; +} + + +/* Arcsecond to radians */ +#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) + + +PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) { + struct pj_opaque_affine *Q = initQ(); + if (0==Q) + return pj_default_destructor(P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + /* read args */ + Q->xoff = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD; + Q->yoff = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD; + Q->zoff = pj_param(P->ctx, P->params, "ddh").f; + + return P; +} diff --git a/src/PJ_geogoffset.c b/src/PJ_geogoffset.c deleted file mode 100644 index 7f917fd7..00000000 --- a/src/PJ_geogoffset.c +++ /dev/null @@ -1,111 +0,0 @@ -/************************************************************************ -* Copyright (c) 2018, Even Rouault -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* 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 -#include - -#include "proj.h" -#include "projects.h" - -PROJ_HEAD(geogoffset, "Geographic Offset"); - -struct pj_opaque_geogoffset { - double dlon; - double dlat; - double dh; -}; - - -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque; - /* offset coordinate */ - obs.lpz.lam += Q->dlon; - obs.lpz.phi += Q->dlat; - obs.lpz.z += Q->dh; - return obs; -} - -static XYZ forward_3d(LPZ lpz, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.lpz = lpz; - return forward_4d(point, P).xyz; -} - - -static XY forward_2d(LP lp, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.lp = lp; - return forward_4d(point, P).xy; -} - - -static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { - struct pj_opaque_geogoffset *Q = (struct pj_opaque_geogoffset *) P->opaque; - /* offset coordinate */ - obs.lpz.lam -= Q->dlon; - obs.lpz.phi -= Q->dlat; - obs.lpz.z -= Q->dh; - return obs; -} - -static LPZ reverse_3d(XYZ xyz, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.xyz = xyz; - return reverse_4d(point, P).lpz; -} - -static LP reverse_2d(XY xy, PJ *P) { - PJ_COORD point = {{0,0,0,0}}; - point.xy = xy; - return reverse_4d(point, P).lp; -} - - -/* Arcsecond to radians */ -#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0) - - -PJ *TRANSFORMATION(geogoffset,0 /* no need for ellipsoid */) { - struct pj_opaque_geogoffset *Q = pj_calloc(1, sizeof(struct pj_opaque_geogoffset)); - if (0==Q) - return pj_default_destructor(P, ENOMEM); - P->opaque = (void *) Q; - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = forward_2d; - P->inv = reverse_2d; - - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - - /* read args */ - Q->dlon = pj_param(P->ctx, P->params, "ddlon").f * ARCSEC_TO_RAD; - Q->dlat = pj_param(P->ctx, P->params, "ddlat").f * ARCSEC_TO_RAD; - Q->dh = pj_param(P->ctx, P->params, "ddh").f; - - return P; -} diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index f1e76ff3..a359369e 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -40,6 +40,7 @@ SET(SRC_LIBPROJ_PJ nad_init.c PJ_aea.c PJ_aeqd.c + PJ_affine.c PJ_airy.c PJ_aitoff.c PJ_august.c @@ -72,7 +73,6 @@ SET(SRC_LIBPROJ_PJ PJ_fouc_s.c PJ_gall.c PJ_geoc.c - PJ_geogoffset.c PJ_geos.c PJ_gins8.c PJ_gnom.c diff --git a/src/makefile.vc b/src/makefile.vc index 21db1362..73c68c0f 100644 --- a/src/makefile.vc +++ b/src/makefile.vc @@ -61,7 +61,7 @@ support = \ pipeline = \ proj_4D_api.obj PJ_cart.obj PJ_pipeline.obj PJ_horner.obj PJ_helmert.obj \ PJ_vgridshift.obj PJ_hgridshift.obj PJ_unitconvert.obj PJ_molodensky.obj \ - PJ_deformation.obj PJ_axisswap.obj PJ_geogoffset.obj + PJ_deformation.obj PJ_axisswap.obj PJ_affine.obj geodesic = geodesic.obj diff --git a/src/pj_list.h b/src/pj_list.h index 6d94f4fd..72a9af7f 100644 --- a/src/pj_list.h +++ b/src/pj_list.h @@ -7,6 +7,7 @@ static const char PJ_LIST_H_ID[] = "@(#)pj_list.h 4.5 95/08/09 GIE REL"; */ PROJ_HEAD(aea, "Albers Equal Area") PROJ_HEAD(aeqd, "Azimuthal Equidistant") +PROJ_HEAD(affine, "Affine transformation") PROJ_HEAD(airy, "Airy") PROJ_HEAD(aitoff, "Aitoff") PROJ_HEAD(alsk, "Mod. Stererographics of Alaska") diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index d3fe64f6..10feea70 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -559,7 +559,7 @@ expect 0 -90 ------------------------------------------------------------------------------- -Test for PJ_geogoffset +Test for PJ_affine ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- operation +proj=geogoffset @@ -589,6 +589,50 @@ accept 10 20 30 40 expect 11 19 33 40 roundtrip 1 +------------------------------------------------------------------------------- +operation +proj=affine +------------------------------------------------------------------------------- +direction forward +tolerance 1mm + +accept 10 20 30 40 +expect 10 20 30 40 +roundtrip 1 + +------------------------------------------------------------------------------- +operation +proj=affine +xoff=1 +yoff=2 +zoff=3 +toff=4 +s11=11 +s12=12 +s13=13 +s21=21 +s22=22 +s23=23 +s31=-31 +s32=32 +s33=33 +tscale=34 +------------------------------------------------------------------------------- +direction forward +tolerance 1mm + +accept 2 49 10 100 +expect 741.0000 1352.0000 1839.0000 3404.0000 +roundtrip 1 + +accept 2 49 10 +expect 741.0000 1352.0000 1839.0000 +roundtrip 1 + +accept 2 49 +expect 611.0000 1122.0000 +roundtrip 1 + +------------------------------------------------------------------------------- +# Non invertible +operation +proj=affine +s11=0 +s22=0 +s23=0 +------------------------------------------------------------------------------- +direction reverse +accept 0 0 0 0 +expect failure + +------------------------------------------------------------------------------- +# Non invertible +operation +proj=affine +tscale=0 +------------------------------------------------------------------------------- +direction reverse +accept 0 0 0 0 +expect failure + ------------------------------------------------------------------------------- run the few gie-builtin tests, which are currently either awkward or impossible to express in the gie command set -- cgit v1.2.3