diff options
| -rw-r--r-- | docs/source/operations/transformations/geogoffset.rst | 68 | ||||
| -rw-r--r-- | docs/source/operations/transformations/index.rst | 1 | ||||
| -rw-r--r-- | docs/source/references.bib | 8 | ||||
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/PJ_geogoffset.c | 111 | ||||
| -rw-r--r-- | src/lib_proj.cmake | 1 | ||||
| -rw-r--r-- | src/makefile.vc | 2 | ||||
| -rw-r--r-- | src/pj_list.h | 1 | ||||
| -rw-r--r-- | test/gie/more_builtins.gie | 32 |
9 files changed, 224 insertions, 2 deletions
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=<value> + + Offset in longitude, expressed in arc-second, to add. + +.. option:: +dlat=<value> + + Offset in latitude, expressed in arc-second, to add. + +.. option:: +dh=<value> + + 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 <even.rouault at spatialys.com> +* +* 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 <errno.h> +#include <math.h> + +#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 |