diff options
Diffstat (limited to '_sources/operations/transformations')
| -rw-r--r-- | _sources/operations/transformations/affine.rst.txt | 123 | ||||
| -rw-r--r-- | _sources/operations/transformations/defmodel.rst.txt | 63 | ||||
| -rw-r--r-- | _sources/operations/transformations/deformation.rst.txt | 229 | ||||
| -rw-r--r-- | _sources/operations/transformations/geogoffset.rst.txt | 70 | ||||
| -rw-r--r-- | _sources/operations/transformations/helmert.rst.txt | 482 | ||||
| -rw-r--r-- | _sources/operations/transformations/hgridshift.rst.txt | 118 | ||||
| -rw-r--r-- | _sources/operations/transformations/horner.rst.txt | 213 | ||||
| -rw-r--r-- | _sources/operations/transformations/index.rst.txt | 24 | ||||
| -rw-r--r-- | _sources/operations/transformations/molobadekas.rst.txt | 147 | ||||
| -rw-r--r-- | _sources/operations/transformations/molodensky.rst.txt | 96 | ||||
| -rw-r--r-- | _sources/operations/transformations/tinshift.rst.txt | 215 | ||||
| -rw-r--r-- | _sources/operations/transformations/vgridshift.rst.txt | 133 | ||||
| -rw-r--r-- | _sources/operations/transformations/xyzgridshift.rst.txt | 93 |
13 files changed, 2006 insertions, 0 deletions
diff --git a/_sources/operations/transformations/affine.rst.txt b/_sources/operations/transformations/affine.rst.txt new file mode 100644 index 00000000..d2a4e26b --- /dev/null +++ b/_sources/operations/transformations/affine.rst.txt @@ -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 coordinate. + ++---------------------+----------------------------------------------------------+ +| **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=<value> + + Offset in X. Default value: 0 + +.. option:: +yoff=<value> + + Offset in Y. Default value: 0 + +.. option:: +zoff=<value> + + Offset in Z. Default value: 0 + +.. option:: +toff=<value> + + Offset in T. Default value: 0 + +.. option:: +s11=<value> + + Rotation/scaling term. Default value: 1 + +.. option:: +s12=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s13=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s21=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s22=<value> + + Rotation/scaling term. Default value: 1 + +.. option:: +s23=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s31=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s32=<value> + + Rotation/scaling term. Default value: 0 + +.. option:: +s33=<value> + + Rotation/scaling term. Default value: 1 + +.. option:: +tscale=<value> + + 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/_sources/operations/transformations/defmodel.rst.txt b/_sources/operations/transformations/defmodel.rst.txt new file mode 100644 index 00000000..77906924 --- /dev/null +++ b/_sources/operations/transformations/defmodel.rst.txt @@ -0,0 +1,63 @@ +.. _defmodel: + +================================================================================ +Multi-component time-based deformation model +================================================================================ + +.. versionadded:: 7.1.0 + ++---------------------+--------------------------------------------------------------------+ +| **Alias** | defmodel | ++---------------------+--------------------------------------------------------------------+ +| **Input type** | Geodetic or projected coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++---------------------+--------------------------------------------------------------------+ +| **Output type** | Geodetic or projected coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++---------------------+--------------------------------------------------------------------+ +| **Domain** | 4D | ++---------------------+--------------------------------------------------------------------+ +| **Available forms** | Forward and inverse | ++---------------------+--------------------------------------------------------------------+ + +The defmodel transformation can be used to represent most deformation models +currently in use, in particular for areas subject to complex deformation, including +large scale secular crustal deformation near plate boundaries and vertical deformation +due to Glacial Isostatic Adjustment (GIA). These can often be represented by a constant +velocity model. Additionally, many areas suffer episodic deformation events such as +earthquakes and aseismic slow slip event. + +The transformation relies on a "master" JSON file, describing general metadata on +the deformation model, its spatial and temporal extent, and listing spatial +components whose values are stored in :ref:`Geodetic TIFF grids (GTG) <geodetictiffgrids>`. +The valuation of each component is modulated by a time function (constant, step, +reverse step, velocity, piecewise, exponential). + +All details on the content of this JSON file are given in the `Proposal for encoding +of a Deformation Model <https://github.com/linz/deformation-model-format/blob/master/doc/JsonGeotiffDeformationModelFormat_20200501.pdf>`__ + +If input coordinates are given in the geographic domain (resp. projected domain), +the output will also be in the geographic domain (resp. projected domain). +The domain should be the corresponding to the source_crs metadata of the model. + +This transformation is a generalization of the :ref:`Kinematic datum shifting utilizing a deformation model <deformation>` transformation. + +Parameters +------------------------------------------------------------------------------- + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +model=<filename> + + Filename to the JSON master file for the deformation model. + + +Example +------------------------------------------------------------------------------- + +Transforming a point with the LINZ NZGD2000 deformation model: + +:: + + echo 166.7133850980 -44.5105886020 293.3700 2007.689 | cct +proj=defmodel +model=nzgd2000-20180701.json diff --git a/_sources/operations/transformations/deformation.rst.txt b/_sources/operations/transformations/deformation.rst.txt new file mode 100644 index 00000000..84b87949 --- /dev/null +++ b/_sources/operations/transformations/deformation.rst.txt @@ -0,0 +1,229 @@ +.. _deformation: + +================================================================================ +Kinematic datum shifting utilizing a deformation model +================================================================================ + +.. versionadded:: 5.0.0 + +Perform datum shifts means of a deformation/velocity model. + ++-----------------+--------------------------------------------------------------------+ +| **Alias** | deformation | ++-----------------+--------------------------------------------------------------------+ +| **Input type** | Cartesian coordinates (spatial), decimalyears (temporal). | ++-----------------+--------------------------------------------------------------------+ +| **Output type** | Cartesian coordinates (spatial), decimalyears (temporal). | ++-----------------+--------------------------------------------------------------------+ +| **Domain** | 4D | ++-----------------+--------------------------------------------------------------------+ +| **Input type** | Geodetic coordinates | ++-----------------+--------------------------------------------------------------------+ +| **Output type** | Geodetic coordinates | ++-----------------+--------------------------------------------------------------------+ + + +The deformation operation is used to adjust coordinates for intraplate deformations. +Usually the transformation parameters for regional plate-fixed reference frames such as +the ETRS89 does not take intraplate deformation into account. It is assumed that +tectonic plate of the region is rigid. Often times this is true, but near the plate +boundary and in areas with post-glacial uplift the assumption breaks. Intraplate +deformations can be modelled and then applied to the coordinates so that +they represent the physical world better. In PROJ this is done with the deformation +operation. + +The horizontal grid is stored in CTable2 format and the vertical grid is stored in the +GTX format. Both grids are expected to contain grid-values in units of +mm/year. GDAL both reads and writes both file formats. Using GDAL for +construction of new grids is recommended. + +Starting with PROJ 7.0, use of a GeoTIFF format is recommended to store both +the horizontal and vertical velocities. + +More complex deformations can be done with the :ref:`Multi-component time-based deformation model <defmodel>` transformation. + +Example +------------------------------------------------------------------------------- + +In :cite:`Hakli2016` coordinate transformation including a deformation model is described. +The paper describes how coordinates from the global ITRFxx frames are transformed to the +local Nordic realisations of ETRS89. Scandinavia is an area with significant post-glacial +rebound. The deformations from the post-glacial uplift is not accounted for in the +official ETRS89 transformations so in order to get accurate transformations in the Nordic +countries it is necessary to apply the deformation model. The transformation from ITRF2008 +to the Danish realisation of ETRS89 is in PROJ described as:: + + + proj = pipeline ellps = GRS80 + # ITRF2008@t_obs -> ITRF2000@t_obs + step init = ITRF2008:ITRF2000 + # ITRF2000@t_obs -> ETRF2000@t_obs + step proj=helmert t_epoch = 2000.0 convention=position_vector + x = 0.054 rx = 0.000891 drx = 8.1e-05 + y = 0.051 ry = 0.00539 dry = 0.00049 + z = -0.048 rz = -0.008712 drz = -0.000792 + # ETRF2000@t_obs -> NKG_ETRF00@2000.0 + step proj = deformation t_epoch = 2000.0 + grids = ./eur_nkg_nkgrf03vel_realigned.tif + inv + # NKG_ETRF@2000.0 -> ETRF92@2000.0 + step proj=helmert convention=position_vector s = -0.009420e + x = 0.03863 rx = 0.00617753 + y = 0.147 ry = 5.064e-05 + z = 0.02776 rz = 4.729e-05 + # ETRF92@2000.0 -> ETRF92@1994.704 + step proj = deformation dt = -5.296 + grids = ./eur_nkg_nkgrf03vel_realigned.tif + +From this we can see that the transformation from ITRF2008 to the Danish realisation of +ETRS89 is a combination of Helmert transformations and adjustments with a deformation +model. The first use of the deformation operation is:: + + proj = deformation t_epoch = 2000.0 grids = ./eur_nkg_nkgrf03vel_realigned.tif + +Here we set the central epoch of the transformation, 2000.0. The observation epoch +is expected as part of the input coordinate tuple. The deformation model is +described by two grids, specified with :option:`+xy_grids` and :option:`+z_grids`. +The first is the horizontal part of the model and the second is the vertical +component. + +Parameters +------------------------------------------------------------------------------- + +.. option:: +xy_grids=<list> + + Comma-separated list of grids to load. If a grid is prefixed by an ``@`` the + grid is considered optional and PROJ will the not complain if the grid is + not available. + + Grids for the horizontal component of a deformation model is expected to be + in CTable2 format. + + .. note:: :option:`+xy_grids` is mutually exclusive with :option:`+grids` + +.. option:: +z_grids=<list> + + Comma-separated list of grids to load. If a grid is prefixed by an `@` the + grid is considered optional and PROJ will the not complain if the grid is + not available. + + Grids for the vertical component of a deformation model is expected to be + in either GTX format. + + .. note:: :option:`+z_grids` is mutually exclusive with :option:`+grids` + +.. option:: +grids=<list> + + .. versionadded:: 7.0.0 + + Comma-separated list of grids to load. If a grid is prefixed by an `@` the + grid is considered optional and PROJ will the not complain if the grid is + not available. + + Grids should be in GeoTIFF format with the first 3 components being + respectively the easting, northing and up velocities in mm/year. + Setting the Description and Unit Type GDAL band metadata items is strongly + recommended, so that gdalinfo reports: + + :: + + Band 1 Block=... Type=Float32, ColorInterp=Gray + Description = east_velocity + Unit Type: mm/year + Band 2 Block=... Type=Float32, ColorInterp=Undefined + Description = north_velocity + Unit Type: mm/year + Band 3 Block=... Type=Float32, ColorInterp=Undefined + Description = up_velocity + Unit Type: mm/year + + .. note:: :option:`+grids` is mutually exclusive with :option:`+xy_grids` + and :option:`+z_grids` + +.. option:: +t_epoch=<value> + + Central epoch of transformation given in decimalyears. Will be used in + conjunction with the observation time from the input coordinate to + determine :math:`dt` as used in eq. :eq:`apply_velocity` below. + + .. note:: :option:`+t_epoch` is mutually exclusive with :option:`+dt` + +.. option:: +dt=<value> + + .. versionadded:: 6.0.0 + + :math:`dt` as used in eq. :eq:`apply_velocity` below. Is useful when + no observation time is available in the input coordinate or when + a deformation for a specific timespan needs to be applied in a + transformation. :math:`dt` is given in units of decimalyears. + + .. note:: :option:`+dt` is mutually exclusive with :option:`+t_epoch` + +Mathematical description +------------------------------------------------------------------------------- + +Mathematically speaking, application of a deformation model is simple. The deformation model is +represented as a grid of velocities in three dimensions. Coordinate corrections are +applied in cartesian space. For a given coordinate, :math:`(X, Y, Z)`, velocities +:math:`(V_X, V_Y, V_Z)` can be interpolated from the gridded model. The time span +between :math:`t_{obs}` and :math:`t_c` determine the magnitude of the coordinate +correction as seen in eq. :eq:`apply_velocity` below. + +.. math:: + :label: apply_velocity + + \begin{align} + \begin{pmatrix} + X \\ + Y \\ + Z \\ + \end{pmatrix}_B = + \begin{pmatrix} + X \\ + Y \\ + Z \\ + \end{pmatrix}_A + + (t_{obs} - t_c) + \begin{pmatrix} + V_X \\ + V_Y \\ + V_Z \\ + \end{pmatrix} + \end{align} + +Corrections are done in cartesian space. + +Coordinates of the gridded model are in ENU (east, north, up) space because it +would otherwise require an enormous 3 dimensional grid to handle the corrections +in cartesian space. Keeping the correction in lat/long space reduces the +complexity of the grid significantly. Consequently though, the input coordinates +needs to be converted to lat/long space when searching for corrections in the +grid. This is done with the :ref:`cart<cart>` operation. The converted grid +corrections can then be applied to the input coordinates in cartesian space. The +conversion from ENU space to cartesian space is done in the following way: + +.. math:: + :label: enu2xyz + + \begin{align} + \begin{pmatrix} + X \\ + Y \\ + Z \\ + \end{pmatrix} = + \begin{pmatrix} + -\sin\phi \cos\lambda N - \sin\lambda E + \cos\phi \cos\lambda U \\ + -\sin\phi \sin\lambda N + \sin\lambda E + \cos\phi \sin\lambda U \\ + \cos\phi N + \sin\phi U \\ + \end{pmatrix} + \end{align} + +where :math:`\phi` and :math:`\lambda` are the latitude and longitude of the coordinate +that is searched for in the grid. :math:`(E, N, U)` are the grid values in ENU-space and +:math:`(X, Y, Z)` are the corrections converted to cartesian space. + + +See also +----------------------------------------------------------------------------- + +#. :ref:`Behavioural changes from version 5 to 6<differences_deformation>` diff --git a/_sources/operations/transformations/geogoffset.rst.txt b/_sources/operations/transformations/geogoffset.rst.txt new file mode 100644 index 00000000..6ca4f64c --- /dev/null +++ b/_sources/operations/transformations/geogoffset.rst.txt @@ -0,0 +1,70 @@ +.. _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:`IOGP2018`) + +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 ellipsoid 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. + +This method is a conveniency wrapper for the more general :ref:`affine`. + +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/_sources/operations/transformations/helmert.rst.txt b/_sources/operations/transformations/helmert.rst.txt new file mode 100644 index 00000000..c0f6934f --- /dev/null +++ b/_sources/operations/transformations/helmert.rst.txt @@ -0,0 +1,482 @@ +.. _helmert: + +================================================================================ +Helmert transform +================================================================================ + +.. versionadded:: 5.0.0 + +The Helmert transformation changes coordinates from one reference frame to +another by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and +14-parameter kinematic counterparts. + + ++-----------------+-------------------------------------------------------------------+ +| **Alias** | helmert | ++-----------------+-------------------------------------------------------------------+ +| **Domain** | 2D, 3D and 4D | ++-----------------+-------------------------------------------------------------------+ +| **Input type** | Cartesian coordinates (spatial), decimalyears (temporal). | ++-----------------+-------------------------------------------------------------------+ +| **Output type** | Cartesian coordinates (spatial), decimalyears (temporal). | ++-----------------+-------------------------------------------------------------------+ +| **Input type** | Cartesian coordinates | ++-----------------+-------------------------------------------------------------------+ +| **Output type** | Cartesian coordinates | ++-----------------+-------------------------------------------------------------------+ + +The Helmert transform, in all its various incarnations, is used to perform reference +frame shifts. The transformation operates in cartesian space. It can be used to transform +planar coordinates from one datum to another, transform 3D cartesian +coordinates from one static reference frame to another or it can be used to do fully +kinematic transformations from global reference frames to local static frames. + +All of the parameters described in the table above are marked as optional. This is true +as long as at least one parameter is defined in the setup of the transformation. +The behavior of the transformation depends on which parameters are used in the setup. +For instance, if a rate of change parameter is specified a kinematic version of the +transformation is used. + +The kinematic transformations require an observation time of the coordinate, as well +as a central epoch for the transformation. The latter is usually documented +alongside the rest of the transformation parameters for a given transformation. +The central epoch is controlled with the parameter `t_epoch`. The observation +time is given as part of the coordinate when using PROJ's 4D-functionality. + +Examples ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Transforming coordinates from NAD72 to NAD83 using the 4 parameter 2D Helmert: + +:: + + proj=helmert convention=coordinate_frame x=-9597.3572 y=.6112 s=0.304794780637 theta=-1.244048 + +Simplified transformations from ITRF2008/IGS08 to ETRS89 using 7 parameters: + +:: + + proj=helmert convention=coordinate_frame x=0.67678 y=0.65495 z=-0.52827 + rx=-0.022742 ry=0.012667 rz=0.022704 s=-0.01070 + +Transformation from `ITRF2000` to `ITRF93` using 15 parameters: + +:: + + proj=helmert convention=position_vector + x=0.0127 y=0.0065 z=-0.0209 s=0.00195 + dx=-0.0029 dy=-0.0002 dz=-0.0006 ds=0.00001 + rx=-0.00039 ry=0.00080 rz=-0.00114 + drx=-0.00011 dry=-0.00019 drz=0.00007 + t_epoch=1988.0 + +Parameters ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. note:: + + All parameters are optional but at least one should be used, otherwise the + operation will return the coordinates unchanged. + +.. option:: +convention=coordinate_frame/position_vector + + .. versionadded:: 5.2.0 + + Indicates the convention to express the rotational terms when a 3D-Helmert / + 7-parameter more transform is involved. As soon as a rotational parameter + is specified (one of ``rx``, ``ry``, ``rz``, ``drx``, ``dry``, ``drz``), + ``convention`` is required. + + The two conventions are equally popular and a frequent source of confusion. + The coordinate frame convention is also described as an clockwise + rotation of the coordinate frame. It corresponds to EPSG method code + 1032 (in the geocentric domain) or 9607 (in the geographic domain) + The position vector convention is also described as an anticlockwise + (counter-clockwise) rotation of the coordinate frame. + It corresponds to as EPSG method code 1033 (in the geocentric domain) or + 9606 (in the geographic domain). + + This parameter is ignored when only a 3-parameter + (translation terms only: ``x``, ``y``, ``z``) , 4-parameter (3-parameter + and ``theta``) or 6-parameter (3-parameter and their derivative terms) + is used. + + The result obtained with parameters specified in a given convention + can be obtained in the other convention by negating the rotational parameters + (``rx``, ``ry``, ``rz``, ``drx``, ``dry``, ``drz``) + + .. note:: This parameter obsoletes ``transpose`` which was present in + PROJ 5.0 and 5.1, and is forbidden starting with PROJ 5.2 + +.. option:: +x=<value> + + Translation of the x-axis given in meters. + +.. option:: +y=<value> + + Translation of the y-axis given in meters. + +.. option:: +z=<value> + + Translation of the z-axis given in meters. + +.. option:: +s=<value> + + Scale factor given in ppm. + +.. option:: +rx=<value> + + X-axis rotation in the 3D Helmert given arc seconds. + + +.. option:: +ry=<value> + + Y-axis rotation in the 3D Helmert given in arc seconds. + +.. option:: +rz=<value> + + Z-axis rotation in the 3D Helmert given in arc seconds. + + +.. option:: +theta=<value> + + Rotation angle in the 2D Helmert given in arc seconds. + +.. option:: +dx=<value> + + Translation rate of the x-axis given in m/year. + +.. option:: +dy=<value> + + Translation rate of the y-axis given in m/year. + +.. option:: +dz=<value> + + Translation rate of the z-axis given in m/year. + +.. option:: +ds=<value> + + Scale rate factor given in ppm/year. + +.. option:: +drx=<value> + + Rotation rate of the x-axis given in arc seconds/year. + +.. option:: +dry=<value> + + Rotation rate of the y-axis given in arc seconds/year. + +.. option:: +drz=<value> + + Rotation rate of the y-axis given in arc seconds/year. + +.. option:: +t_epoch=<value> + + Central epoch of transformation given in decimalyear. Only used + spatiotemporal transformations. + +.. option:: +exact + + Use exact transformation equations. + + See :eq:`rot_exact` + +.. option:: +transpose + + .. deprecated:: 5.2.0 (removed) + + Transpose rotation matrix and follow the **Position Vector** rotation + convention. If :option:`+transpose` is not added the **Coordinate Frame** + rotation convention is used. + + + +Mathematical description ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +In the notation used below, :math:`\hat{P}` is the rate of change of a given transformation +parameter :math:`P`. :math:`\dot{P}` is the kinematically adjusted version of :math:`P`, +described by + +.. math:: + :label: propagation + + \dot{P}= P + \hat{P}\left(t - t_{central}\right) + +where :math:`t` is the observation time of the coordinate and :math:`t_{central}` is +the central epoch of the transformation. Equation :eq:`propagation` can be used to +propagate all transformation parameters in time. + +Superscripts of vectors denote the reference frame the coordinates in the vector belong to. + + +2D Helmert +------------------------------------------------------------------------------- + +The simplest version of the Helmert transform is the 2D case. In the 2-dimensional +case only the horizontal coordinates are changed. The coordinates can be +translated, rotated and scale. Translation is controlled with the `x` and `y` +parameters. The rotation is determined by `theta` and the scale is controlled with +the `s` parameters. + +.. note:: + + The scaling parameter `s` is unitless for the 2D Helmert, as opposed to the + 3D version where the scaling parameter is given in units of ppm. + +Mathematically the 2D Helmert is described as: + +.. math:: + :label: 4param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + \end{bmatrix}^B = + \begin{bmatrix} + T_x \\ + T_y \\ + \end{bmatrix} + + s + \begin{bmatrix} + \hphantom{-}\cos \theta & \sin \theta \\ + -\sin \theta & \cos \theta \\ + \end{bmatrix} + \begin{bmatrix} + X \\ + Y \\ + \end{bmatrix}^A + \end{align} + + +:eq:`4param` can be extended to a time-varying kinematic version by +adjusting the parameters with :eq:`propagation` to :eq:`4param`, which yields +the kinematic 2D Helmert transform: + +.. math:: + :label: 8param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + \end{bmatrix}^B = + \begin{bmatrix} + \dot{T_x} \\ + \dot{T_y} \\ + \end{bmatrix} + + s(t) + \begin{bmatrix} + \hphantom{-}\cos \dot{\theta} & \sin \dot{\theta} \\ + -\sin\ \dot{\theta} & \cos \dot{\theta} \\ + \end{bmatrix} + \begin{bmatrix} + X \\ + Y \\ + \end{bmatrix}^A + \end{align} + +All parameters in :eq:`8param` are determined by the use of :eq:`propagation`, +which applies the rate of change to each individual parameter for a given +timespan between :math:`t` and :math:`t_{central}`. + + +3D Helmert +------------------------------------------------------------------------------- + +The general form of the 3D Helmert is + +.. math:: + :label: general-helmert + + + \begin{align} + V^B = T + \left(1 + s \times 10^{-6}\right) \mathbf{R} V^A + \end{align} + +Where :math:`T` is a vector consisting of the three translation parameters, :math:`s` +is the scaling factor and :math:`\mathbf{R}` is a rotation matrix. :math:`V^A` and +:math:`V^B` are coordinate vectors, with :math:`V^A` being the input coordinate and +:math:`V^B` is the output coordinate. + +In the *Position Vector* convention, we define :math:`R_x = radians \left( rx \right)`, +:math:`R_z = radians \left( ry \right)` and :math:`R_z = radians \left( rz \right)` + +In the *Coordinate Frame* convention, :math:`R_x = - radians \left( rx \right)`, +:math:`R_z = - radians \left( ry \right)` and :math:`R_z = - radians \left( rz \right)` + +The rotation matrix is composed of three rotation matrices, one for each axis. + +.. math:: + + \begin{align} + \mathbf{R}_X &= \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos R_x & -\sin R_x \\ 0 & \sin R_x & \cos R_x \end{bmatrix} + \end{align} + +.. math:: + + \begin{align} + \mathbf{R}_Y &= \begin{bmatrix} \cos R_y & 0 & \sin R_y\\ 0 & 1 & 0\\ -\sin R_y & 0 & \cos R_y \end{bmatrix} + \end{align} + +.. math:: + + \begin{align} + \mathbf{R}_Z &= \begin{bmatrix} \cos R_z & -\sin R_z & 0\\ \sin R_z & \cos R_z & 0\\ 0 & 0 & 1 \end{bmatrix} + \end{align} + +The three rotation matrices can be combined in one: + +.. math:: + + \begin{align} + \mathbf{R} = \mathbf{R_X} \mathbf{R_Y} \mathbf{R_Y} + \end{align} + + +For :math:`\mathbf{R}`, this yields: + +.. math:: + :label: rot_exact + + \begin{bmatrix} + \cos R_y \cos R_z & -\cos R_x \sin R_z + & \sin R_x \sin R_z + \\ + & \sin R_x \sin R_y \cos R_z & \cos R_x \sin R_y \cos R_z \\ + \cos R_y\sin R_z & \cos R_x \cos R_z + & - \sin R_x \cos R_z + \\ + & \sin R_x \sin R_y \sin R_z & \cos R_x \sin R_y \sin R_z \\ + -\sin R_y & \sin R_x \cos R_y & \cos R_x \cos R_y \\ + \end{bmatrix} + + +Using the small angle approximation the rotation matrix can be simplified to + +.. math:: + :label: rot_approx + + \begin{align} \mathbf{R} = + \begin{bmatrix} + 1 & -R_z & R_y \\ + Rz & 1 & -R_x \\ + -Ry & R_x & 1 \\ + \end{bmatrix} + \end{align} + +Which allow us to express the most common version of the Helmert transform, +using the approximated rotation matrix: + + +.. math:: + :label: 7param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^B = + \begin{bmatrix} + T_x \\ + T_y \\ + T_z \\ + \end{bmatrix} + + \left(1 + s \times 10^{-6}\right) + \begin{bmatrix} + 1 & -R_z & R_y \\ + Rz & 1 & -R_x \\ + -Ry & R_x & 1 \\ + \end{bmatrix} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^A + \end{align} + +If the rotation matrix is transposed, or the sign of the rotation terms negated, +the rotational part of the transformation is effectively reversed. +This is what happens when switching between the 2 conventions ``position_vector`` +and ``coordinate_frame`` + +Applying :eq:`propagation` we get the kinematic version of the approximated +3D Helmert: + +.. math:: + :label: 14param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^B = + \begin{bmatrix} + \dot{T_x} \\ + \dot{T_y} \\ + \dot{T_z} \\ + \end{bmatrix} + + \left(1 + \dot{s} \times 10^{-6}\right) + \begin{bmatrix} + 1 & -\dot{R_z} & \dot{R_y} \\ + \dot{R_z} & 1 & -\dot{R_x} \\ + -\dot{R_y} & \dot{R_x} & 1 \\ + \end{bmatrix} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^A + \end{align} + + + + +The Helmert transformation can be applied without using the rotation parameters, +in which case it becomes a simple translation of the origin of the coordinate +system. When using the Helmert in this version equation :eq:`general-helmert` +simplifies to: + +.. math:: + :label: 3param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^B = + \begin{bmatrix} + T_x \\ + T_y \\ + T_z \\ + \end{bmatrix} + + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^A + \end{align} + +That after application of :eq:`propagation` has the following kinematic +counterpart: + +.. math:: + :label: 6param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^B = + \begin{bmatrix} + \dot{T_x} \\ + \dot{T_y} \\ + \dot{T_z} \\ + \end{bmatrix} + + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^A + \end{align} diff --git a/_sources/operations/transformations/hgridshift.rst.txt b/_sources/operations/transformations/hgridshift.rst.txt new file mode 100644 index 00000000..d27059e9 --- /dev/null +++ b/_sources/operations/transformations/hgridshift.rst.txt @@ -0,0 +1,118 @@ +.. _hgridshift: + +================================================================================ +Horizontal grid shift +================================================================================ + +.. versionadded:: 5.0.0 + +Change of horizontal datum by grid shift. + ++-----------------+--------------------------------------------------------------------+ +| **Alias** | hgridshift | ++-----------------+--------------------------------------------------------------------+ +| **Domain** | 2D, 3D and 4D | ++-----------------+--------------------------------------------------------------------+ +| **Input type** | Geodetic coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++-----------------+--------------------------------------------------------------------+ +| **Output type** | Geodetic coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++-----------------+--------------------------------------------------------------------+ + +The horizontal grid shift is done by offsetting the planar input coordinates by +a specific amount determined by the loaded grids. The simplest use case of the +horizontal grid shift is applying a single grid:: + + +proj=hgridshift +grids=nzgr2kgrid0005.gsb + + +More than one grid can be loaded at the same time, for instance in case the +dataset needs to be transformed spans several countries. In this example grids +of the continental US, Alaska and Canada is loaded at the same time:: + + +proj=hgridshift +grids=@conus,@alaska,@ntv2_0.gsb,@ntv_can.dat + +The ``@`` in the above example states that the grid is optional, in case the grid +is not found in the PROJ search path. The list of grids is prioritized so that +grids in the start of the list takes precedence over the grids in the back of the +list. + +PROJ supports CTable2, NTv1 and NTv2 files for horizontal grid corrections. Details +about all three formats can be found in the GDAL documentation and/or driver source +code. GDAL reads and writes all three formats. Using GDAL for construction of +new grids is recommended. + + +Temporal gridshifting +################################################################################ +.. versionadded:: 5.1.0 + +By initializing the horizontal gridshift operation with a central epoch, it can be +used as a step function applying the grid offsets only if a coordinate is +transformed from an epoch before grids central epoch to an epoch after. This is +handy in transformations where it is necessary to handle deformations caused by +seismic activity. + +The central epoch of the grid is controlled with :option:`+t_epoch` and the final +epoch of the coordinate is set with :option:`+t_final`. The observation epoch of +the coordinate is part of the coordinate tuple. + +Suppose we want to model the deformation of the 2008 earthquake in Iceland in +a transformation of data from 2005 to 2009:: + + echo 63.992 -21.014 10.0 2005.0 | cct +proj=hgridshift +grids=iceland2008.gsb +t_epoch=2008.4071 +t_final=2009.0 + 63.9920021 -21.0140013 10.0 2005.0 + +.. note:: + The timestamp of the resulting coordinate is still 2005.0. The observation + time is always kept unchanged as it would otherwise be impossible to do the + inverse transformation. + +Temporal gridshifting is especially powerful in transformation pipelines where +several gridshifts can be chained together, effectively acting as a series of +step functions that can be applied to a coordinate that is propagated through +time. In the following example we establish a pipeline that allows transformation +of coordinates from any given epoch up until the current date, applying only +those gridshifts that have central epochs between the observation epoch and +the final epoch:: + + +proj=pipeline +t_final=now + +step +proj=hgridshift +grids=earthquake_1.gsb +t_epoch=2010.421 + +step +proj=hgridshift +grids=earthquake_2.gsb +t_epoch=2013.853 + +step +proj=hgridshift +grids=earthquake_3.gsb +t_epoch=2017.713 + +.. note:: + The special epoch *now* is used when specifying the final epoch with + :option:`+t_final`. This results in coordinates being transformed to the + current date. Additionally, :option:`+t_final` is used as a + :ref:`global pipeline parameter<global-pipeline-parameter>`, which means + that it is applied to all the steps in the pipeline. + +In the above transformation, a coordinate with observation epoch 2009.32 would +be subject to all three gridshift steps in the pipeline. A coordinate with +observation epoch 2014.12 would only by offset by the last step in the pipeline. + + +Parameters +################################################################################ + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +grids=<list> + + Comma-separated list of grids to load. If a grid is prefixed by an `@` the + grid is considered optional and PROJ will the not complain if the grid is + not available. + + Grids are expected to be in CTable2, NTv1 or NTv2 format. + +Optional ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. include:: ../options/t_epoch.rst +.. versionadded:: 5.1.0 +.. include:: ../options/t_final.rst +.. versionadded:: 5.1.0 + diff --git a/_sources/operations/transformations/horner.rst.txt b/_sources/operations/transformations/horner.rst.txt new file mode 100644 index 00000000..ad09e8ab --- /dev/null +++ b/_sources/operations/transformations/horner.rst.txt @@ -0,0 +1,213 @@ +.. _horner: + +================================================================================ +Horner polynomial evaluation +================================================================================ + +.. versionadded:: 5.0.0 + ++-----------------+-------------------------------------------------------------------+ +| **Alias** | horner | ++-----------------+-------------------------------------------------------------------+ +| **Domain** | 2D and 3D | ++-----------------+-------------------------------------------------------------------+ +| **Input type** | Geodetic and projected coordinates | ++-----------------+-------------------------------------------------------------------+ +| **Output type** | Geodetic and projected coordinates | ++-----------------+-------------------------------------------------------------------+ + +The Horner polynomial evaluation scheme is used for transformations between reference +frames where one or both are inhomogeneous or internally distorted. This will typically +be reference frames created before the introduction of space geodetic techniques such +as GPS. + +Horner polynomials, or Multiple Regression Equations as they are also known as, have +their strength in being able to create complicated mappings between coordinate +reference frames while still being lightweight in both computational cost and disk +space used. + +PROJ implements two versions of the Horner evaluation scheme: Real and complex +polynomial evaluation. Below both are briefly described. For more details consult +:cite:`Ruffhead2016` and :cite:`IOGP2018`. + +The polynomial evaluation in real number space is defined by the following equations: + +.. math:: + :label: real_poly + + \Delta X = \sum_{i,j} u_{i,j} U^i V^j + + \Delta Y = \sum_{i,j} v_{i,j} U^i V^j + +where + +.. math:: + :label: UV + + U = X_{in} - X_{origin} + + V = Y_{in} - Y_{origin} + +and :math:`u_{i,j}` and :math:`v_{i,j}` are coefficients that make up +the polynomial. + +The final coordinates are determined as + +.. math:: + :label: xy_out + + X_{out} = X_{in} + \Delta X + + Y_{out} = Y_{in} + \Delta Y + +The inverse transform is the same as the above but requires a different set of +coefficients. + +Evaluation of the complex polynomials are defined by the following equations: + +.. math:: + :label: complex_poly + + \Delta X + i \Delta Y = \sum_{j=1}^n (c_{2j-1} + i c_{2j}) (U + i V)^j + +Where :math:`n` is the degree of the polynomial. :math:`U` and :math:`V` are +defined as in :eq:`UV` and the resulting coordinates are again determined +by :eq:`xy_out`. + +Examples +################################################################################ + +Mapping between Danish TC32 and ETRS89/UTM zone 32 using polynomials in real +number space:: + + +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 + + +Mapping between Danish System Storebælt and ETRS89/UTM zone 32 using complex +polynomials:: + + +proj=horner + +ellps=intl + +range=500000 + +fwd_origin=4.94690026817276e+05,6.13342113183056e+06 + +inv_origin=6.19480258923588e+05,6.13258568148837e+06 + +deg=3 + +fwd_c=6.13258562111350e+06,6.19480105709997e+05,9.99378966275206e-01,-2.82153291753490e-02,-2.27089979140026e-10,-1.77019590701470e-09,1.08522286274070e-14,2.11430298751604e-15 + +inv_c=6.13342118787027e+06,4.94690181709311e+05,9.99824464710368e-01,2.82279070814774e-02,7.66123542220864e-11,1.78425334628927e-09,-1.05584823306400e-14,-3.32554258683744e-15 + +Parameters +################################################################################ + +Setting up Horner polynomials requires many coefficients being explicitly +written, even for polynomials of low degree. For this reason it is recommended +to store the polynomial definitions in an :ref:`init file <init_files>` for +easier writing and reuse. + +Required +------------------------------------------------------------------------------- + +Below is a list of required parameters that can be set for the Horner polynomial +transformation. As stated above, the transformation takes to forms, either using +real or complex polynomials. These are divided into separate sections below. +Parameters from the two sections are mutually exclusive, that is parameters +describing real and complex polynomials can't be mixed. + +.. include:: ../options/ellps.rst + +.. option:: +deg=<value> + + Degree of polynomial + +.. option:: +fwd_origin=<northing,easting> + + Coordinate of origin for the forward mapping + +.. option:: +inv_origin=<northing,easting> + + Coordinate of origin for the inverse mapping + +Real polynomials +.............................................................................. + +The following parameters has to be set if the transformation consists of +polynomials in real space. Each parameter takes a comma-separated list of +coefficients. The number of coefficients is governed by the degree, :math:`d`, +of the polynomial: + +.. math:: + + N = \frac{(d + 1)(d + 2)}{2} + +.. option:: +fwd_u=<u_11,u_12,...,u_ij,..,u_mn> + + Coefficients for the forward transformation i.e. latitude to northing + as described in :eq:`real_poly`. + +.. option:: +fwd_v=<v_11,v_12,...,v_ij,..,v_mn> + + Coefficients for the forward transformation i.e. longitude to easting + as described in :eq:`real_poly`. + +.. option:: +inv_u=<u_11,u_12,...,u_ij,..,u_mn> + + Coefficients for the inverse transformation i.e. latitude to northing + as described in :eq:`real_poly`. + +.. option:: +inv_v=<v_11,v_12,...,v_ij,..,v_mn> + + Coefficients for the inverse transformation i.e. longitude to easting + as described in :eq:`real_poly`. + + + +Complex polynomials +.............................................................................. + +The following parameters has to be set if the transformation consists of +polynomials in complex space. Each parameter takes a comma-separated list of +coefficients. The number of coefficients is governed by the degree, :math:`d`, +of the polynomial: + +.. math:: + + N = 2d + 2 + +.. option:: +fwd_c=<c_1,c_2,...,c_N> + + Coefficients for the complex forward transformation + as described in :eq:`complex_poly`. + +.. option:: +inv_c=<c_1,c_2,...,c_N> + + Coefficients for the complex inverse transformation + as described in :eq:`complex_poly`. + +Optional +------------------------------------------------------------------------------- + +.. option:: +range=<value> + + Radius of the region of validity. + +.. option:: +uneg + + Express latitude as southing. Only applies for complex polynomials. + +.. option:: +vneg + + Express longitude as westing. Only applies for complex polynomials. + + +Further reading +################################################################################ + +#. `Wikipedia <https://en.wikipedia.org/wiki/Horner%27s_method>`_ diff --git a/_sources/operations/transformations/index.rst.txt b/_sources/operations/transformations/index.rst.txt new file mode 100644 index 00000000..a413b8f6 --- /dev/null +++ b/_sources/operations/transformations/index.rst.txt @@ -0,0 +1,24 @@ +.. _transformation_list: + +================================================================================ +Transformations +================================================================================ + +Transformations coordinate operation in which the two coordinate reference +systems are based on different datums. + +.. toctree:: + :maxdepth: 1 + + affine + defmodel + deformation + geogoffset + helmert + horner + molodensky + molobadekas + hgridshift + tinshift + vgridshift + xyzgridshift diff --git a/_sources/operations/transformations/molobadekas.rst.txt b/_sources/operations/transformations/molobadekas.rst.txt new file mode 100644 index 00000000..b7d638bf --- /dev/null +++ b/_sources/operations/transformations/molobadekas.rst.txt @@ -0,0 +1,147 @@ +.. _molobadekas: + +================================================================================ +Molodensky-Badekas transform +================================================================================ + +.. versionadded:: 6.0.0 + +The Molodensky-Badekas transformation changes coordinates from one reference frame to +another by means of a 10-parameter shift. + +.. note:: + + It should not be confused with the :ref:`Molodensky` transform which + operates directly in the geodetic coordinates. Molodensky-Badekas can rather + be seen as a variation of :ref:`Helmert` + ++-----------------+-------------------------------------------------------------------+ +| **Alias** | molobadekas | ++-----------------+-------------------------------------------------------------------+ +| **Domain** | 3D | ++-----------------+-------------------------------------------------------------------+ +| **Input type** | Cartesian coordinates | ++-----------------+-------------------------------------------------------------------+ +| **Output type** | Cartesian coordinates | ++-----------------+-------------------------------------------------------------------+ + +The Molodensky-Badekas transformation is a variation of the :ref:`Helmert` where +the rotational terms are not directly applied to the ECEF coordinates, but on +cartesian coordinates relative to a reference point (usually close to Earth surface, +and to the area of use of the transformation). When ``px`` = ``py`` = ``pz`` = 0, +this is equivalent to a 7-parameter Helmert transformation. + +Example ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Transforming coordinates from La Canoa to REGVEN: + +:: + + proj=molobadekas convention=coordinate_frame + x=-270.933 y=115.599 z=-360.226 rx=-5.266 ry=-1.238 rz=2.381 + s=-5.109 px=2464351.59 py=-5783466.61 pz=974809.81 + + +Parameters ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. note:: + + All parameters (except convention) are optional but at least one should be + used, otherwise the operation will return the coordinates unchanged. + +.. option:: +convention=coordinate_frame/position_vector + + Indicates the convention to express the rotational terms when a 3D-Helmert / + 7-parameter more transform is involved. + + The two conventions are equally popular and a frequent source of confusion. + The coordinate frame convention is also described as an clockwise + rotation of the coordinate frame. It corresponds to EPSG method code + 1034 (in the geocentric domain) or 9636 (in the geographic domain) + The position vector convention is also described as an anticlockwise + (counter-clockwise) rotation of the coordinate frame. + It corresponds to as EPSG method code 1061 (in the geocentric domain) or + 1063 (in the geographic domain). + + The result obtained with parameters specified in a given convention + can be obtained in the other convention by negating the rotational parameters + (``rx``, ``ry``, ``rz``) + +.. option:: +x=<value> + + Translation of the x-axis given in meters. + +.. option:: +y=<value> + + Translation of the y-axis given in meters. + +.. option:: +z=<value> + + Translation of the z-axis given in meters. + +.. option:: +s=<value> + + Scale factor given in ppm. + +.. option:: +rx=<value> + + X-axis rotation given arc seconds. + +.. option:: +ry=<value> + + Y-axis rotation given in arc seconds. + +.. option:: +rz=<value> + + Z-axis rotation given in arc seconds. + +.. option:: +px=<value> + + Coordinate along the x-axis of the reference point given in meters. + +.. option:: +py=<value> + + Coordinate along the y-axis of the reference point given in meters. + +.. option:: +pz=<value> + + Coordinate along the z-axis of the reference point given in meters. + +Mathematical description ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + +In the *Position Vector* convention, we define :math:`R_x = radians \left( rx \right)`, +:math:`R_z = radians \left( ry \right)` and :math:`R_z = radians \left( rz \right)` + +In the *Coordinate Frame* convention, :math:`R_x = - radians \left( rx \right)`, +:math:`R_z = - radians \left( ry \right)` and :math:`R_z = - radians \left( rz \right)` + +.. math:: + :label: 10param + + \begin{align} + \begin{bmatrix} + X \\ + Y \\ + Z \\ + \end{bmatrix}^{output} = + \begin{bmatrix} + T_x + P_x\\ + T_y + P_y\\ + T_z + P_z\\ + \end{bmatrix} + + \left(1 + s \times 10^{-6}\right) + \begin{bmatrix} + 1 & -R_z & R_y \\ + Rz & 1 & -R_x \\ + -Ry & R_x & 1 \\ + \end{bmatrix} + \begin{bmatrix} + X^{input} - P_x\\ + Y^{input} - P_y\\ + Z^{input} - P_z\\ + \end{bmatrix} + \end{align} diff --git a/_sources/operations/transformations/molodensky.rst.txt b/_sources/operations/transformations/molodensky.rst.txt new file mode 100644 index 00000000..df0b00a2 --- /dev/null +++ b/_sources/operations/transformations/molodensky.rst.txt @@ -0,0 +1,96 @@ +.. _molodensky: + +================================================================================ +Molodensky transform +================================================================================ + +.. versionadded:: 5.0.0 + +The Molodensky transformation resembles a :ref:`Helmert` with zero +rotations and a scale of unity, but converts directly from geodetic coordinates to +geodetic coordinates, without the intermediate shifts to and from cartesian +geocentric coordinates, associated with the Helmert transformation. +The Molodensky transformation is simple to implement and to parametrize, requiring +only the 3 shifts between the input and output frame, and the corresponding +differences between the semimajor axes and flattening parameters of the reference +ellipsoids. Due to its algorithmic simplicity, it was popular prior to the +ubiquity of digital computers. Today, it is mostly interesting for historical +reasons, but nevertheless indispensable due to the large amount of data that has +already been transformed that way :cite:`EversKnudsen2017`. + ++---------------------+----------------------------------------------------------+ +| **Alias** | molodensky | ++---------------------+----------------------------------------------------------+ +| **Domain** | 3D | ++---------------------+----------------------------------------------------------+ +| **Input type** | Geodetic coordinates (horizontal), meters (vertical) | ++---------------------+----------------------------------------------------------+ +| **output type** | Geodetic coordinates (horizontal), meters (vertical) | ++---------------------+----------------------------------------------------------+ + +The Molodensky transform can be used to perform a datum shift from coordinate +:math:`(\phi_1, \lambda_1, h_1)` to :math:`(\phi_2, \lambda_2, h_2)` where the two +coordinates are referenced to different ellipsoids. This is based on three +assumptions: + + 1. The cartesian axes, :math:`X, Y, Z`, of the two ellipsoids are parallel. + 2. The offset, :math:`\delta X, \delta Y, \delta Z`, between the two ellipsoid + are known. + 3. The characteristics of the two ellipsoids, expressed as the difference in + semimajor axis (:math:`\delta a`) and flattening (:math:`\delta f`), are known. + +The Molodensky transform is mostly used for transforming between old systems +dating back to the time before computers. The advantage of the Molodensky transform +is that it is fairly simple to compute by hand. The ease of computation come at the +cost of limited accuracy. + +A derivation of the mathematical formulas for the Molodensky transform can be found +in :cite:`Deakin2004`. + + +Examples +############################################################################### + +The abridged Molodensky:: + + proj=molodensky a=6378160 rf=298.25 da=-23 df=-8.120449e-8 dx=-134 dy=-48 dz=149 abridged + +The same transformation using the standard Molodensky:: + + proj=molodensky a=6378160 rf=298.25 da=-23 df=-8.120449e-8 dx=-134 dy=-48 dz=149 + + +Parameters +################################################################################ + +Required +------------------------------------------------------------------------------- + +.. option:: +da=<value> + + Difference in semimajor axis of the defining ellipsoids. + +.. option:: +df=<value> + + Difference in flattening of the defining ellipsoids. + +.. option:: +dx=<value> + + Offset of the X-axes of the defining ellipsoids. + +.. option:: +dy=<value> + + Offset of the Y-axes of the defining ellipsoids. + +.. option:: +dz=<value> + + Offset of the Z-axes of the defining ellipsoids. + +.. include:: ../options/ellps.rst + +Optional +------------------------------------------------------------------------------- + +.. option:: +abridged + + Use the abridged version of the Molodensky transform. diff --git a/_sources/operations/transformations/tinshift.rst.txt b/_sources/operations/transformations/tinshift.rst.txt new file mode 100644 index 00000000..2135266b --- /dev/null +++ b/_sources/operations/transformations/tinshift.rst.txt @@ -0,0 +1,215 @@ +.. _tinshift: + +================================================================================ +Triangulated Irregular Network based transformation +================================================================================ + +.. versionadded:: 7.2.0 + ++---------------------+--------------------------------------------------------------------+ +| **Alias** | tinshift | ++---------------------+--------------------------------------------------------------------+ +| **Input type** | Projected or geographic coordinates (horizontal), meters (vertical)| ++---------------------+--------------------------------------------------------------------+ +| **Output type** | Projected or geographic coordinates (horizontal), meters (vertical)| ++---------------------+--------------------------------------------------------------------+ +| **Domain** | 2D or 3D | ++---------------------+--------------------------------------------------------------------+ +| **Available forms** | Forward and inverse | ++---------------------+--------------------------------------------------------------------+ + + +The ``tinshift`` transformation takes one mandatory +argument, ``file``, that points to a JSON file, which contains the +triangulation and associated metadata. Input and output coordinates must be +geographic or projected coordinates. +Depending on the content of the JSON file, horizontal, vertical or both +components of the coordinates may be transformed. + +The transformation is invertible, with the same computational complexity than +the forward transformation. + +Parameters +------------------------------------------------------------------------------- + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +file=<filename> + + Filename to the JSON file for the TIN. + + +Example +------------------------------------------------------------------------------- + +Transforming a point with the Finland EPSG:2393 ("KKJ / Finland Uniform +Coordinate System") projected CRS to EPSG:3067 ("ETRS89 / TM35FIN(E,N)") + +:: + + $ echo 3210000.0000 6700000.0000 0 2020 | cct +proj=tinshift +file=./triangulation_kkj.json + + 209948.3217 6697187.0009 0.0000 2020 + + +Algorithm ++++++++++ + +Internally, ``tinshift`` ingest the whole file into memory. It is considered that +triangulation should be small enough for that. + +When a point is transformed, one must find the triangle into which it falls into. +Instead of iterating over all triangles, we build a in-memory quadtree to speed-up +the identification of candidates triangles. + +To determine if a point falls into a triangle, one computes its 3 +`barycentric coordinates <https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Conversion_between_barycentric_and_Cartesian_coordinates>`_ +from its projected coordinates, :math:`\lambda_i` for :math:`i=1,2,3`. +They are real values (in the [0,1] range for a point inside the triangle), +giving the weight of each of the 3 vertices of the triangles. + +Once those weights are known, interpolating the target horizontal +coordinate is a matter of doing the linear combination of those weights with +the target horizontal coordinates at the 3 vertices of the triangle (:math:`Xt_i` and :math:`Yt_i`): + +.. math:: + + X_{target} = Xt_1 * \lambda_1 + Xt_2 * \lambda_2 + Xt_3 * \lambda_3 + + Y_{target} = Yt_1 * \lambda_1 + Yt_2 * \lambda_2 + Yt_3 * \lambda_3 + +This interpolation is exact at the vertices of the triangulation, and has linear properties +inside each triangle. It is completely equivalent to other formulations of +triangular interpolation, such as + +.. math:: + + X_{target} = A + X_{source} * B + Y_{source} * C + + Y_{target} = D + X_{source} * E + Y_{source} * F + +where the A, B, C, D, E, F constants (for a given triangle) are found by solving +the 2 systems of 3 linear equations, constraint by the source and target coordinate pairs +of the 3 vertices of the triangle: + +.. math:: + + Xt_i = A + Xs_i * B + Ys_i * C + + Yt_i = D + Xs_i * E + Ys_i * F + +Similarly for a vertical coordinate transformation, where :math:`Zoff_i` is the vertical +offset at each vertex of the triangle: + +.. math:: + + Z_{target} = Z_{source} + Zoff_1 * \lambda_1 + Zoff_2 * \lambda_2 + Zoff_3 * \lambda_3 + +Constraints on the triangulation +++++++++++++++++++++++++++++++++ + +No check is done on the consistence of the triangulation. It is highly +recommended that triangles do not overlap each other (when considering the +source coordinates or the forward transformation, or the target coordinates for +the inverse transformation), otherwise which triangle will be selected is +unspecified. Besides that, the triangulation does not need to have particular +properties (like being a Delaunay triangulation) + +File format ++++++++++++ + +The triangulation is stored in a text-based format, using JSON as a serialization. + +Below a minimal example, from the KKJ to ETRS89 transformation, with just a +single triangle: + +.. literalinclude:: ../../../../data/tests/tinshift_crs_implicit.json + :language: json + + +So after the generic metadata, we define the input and output CRS (informative +only), and that the transformation affects horizontal components of +coordinates. We name the columns of the ``vertices`` and ``triangles`` arrays. +We defined the source and target coordinates of each vertex, and define a +triangle by referring to the index of its vertices in the ``vertices`` array. + +More formally, the specific items for the triangulation file are: + +input_crs + String identifying the CRS of source coordinates + in the vertices. Typically ``EPSG:XXXX``. If the transformation is for vertical + component, this should be the code for a compound CRS (can be EPSG:XXXX+YYYY + where XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). + For example, for the KKJ->ETRS89 transformation, this is EPSG:2393 + (``KKJ / Finland Uniform Coordinate System``). The input coordinates are assumed + to be passed in the "normalized for visualisation" / "GIS friendly" order, + that is longitude, latitude for geographic coordinates and + easting, northing for projected coordinates. + + +output_crs + String identifying the CRS of target coordinates in the vertices. + Typically ``EPSG:XXXX``. If the transformation is for vertical component, + this should be the code for a compound CRS (can be EPSG:XXXX+YYYY where + XXXX is the code of the horizontal CRS and YYYY the code of the vertical CRS). + For example, for the KKJ->ETRS89 transformation, this is EPSG:3067 + (\"ETRS89 / TM35FIN(E,N)\"). The output coordinates will be returned in + the "normalized for visualisation" / "GIS friendly" order, + that is longitude, latitude for geographic coordinates and + easting, northing for projected coordinates. + + +transformed_components + Array which may contain one or two strings: "horizontal" when horizontal + components of the coordinates are transformed and/or "vertical" when the + vertical component is transformed. + + +fallback_strategy + String identifying how to treat points that do not fall into any of the + specified triangles. This item is available for ``format_version`` >= 1.1. + Possible values are ``none``, ``nearest_side`` and ``nearest_centroid``. The + default is ``none`` and signifies, that points that fall outside the + specified triangles are not transformed. This is also the behaviour for + ``format_version`` before 1.1. If ``fallback_strategy`` is set to + ``nearest_side``, then points that do not fall into any triangle are + transformed according to the triangle closest to them by euclidean distance. + If ``fallback_strategy`` is set to ``nearest_centroid``, then points that do + not fall into any triangle are transformed according to the triangle with the + closest centroid (intersection of its medians). + + +vertices_columns + Specify the name of the columns of the rows in the ``vertices`` + array. There must be exactly as many elements in ``vertices_columns`` as in a + row of ``vertices``. The following names have a special meaning: ``source_x``, + ``source_y``, ``target_x``, ``target_y``, ``source_z``, ``target_z`` and + ``offset_z``. ``source_x`` and ``source_y`` are compulsory. + ``source_x`` is for the source longitude (in degree) or easting. + ``source_y`` is for the source latitude (in degree) or northing. + ``target_x`` and ``target_y`` are compulsory when ``horizontal`` is specified + in ``transformed_components``. (``source_z`` and ``target_z``) or + ``offset_z`` are compulsory when ``vertical`` is specified in ``transformed_components`` + + +triangles_columns + Specify the name of the columns of the rows in the + ``triangles`` array. There must be exactly as many elements in ``triangles_columns`` + as in a row of ``triangles``. The following names have a special meaning: + ``idx_vertex1``, ``idx_vertex2``, ``idx_vertex3``. They are compulsory. + + +vertices + An array whose items are themselves arrays with as many columns as + described in ``vertices_columns``. + + +triangles + An array whose items are themselves arrays with as many columns as + described in ``triangles_columns``. + The value of the ``idx_vertexN`` columns must be indices + (between 0 and len(``vertices``-1) of items of the ``vertices`` array. + +A `JSON schema <https://proj.org/schemas/triangulation.schema.json>`_ is available +for this file format. diff --git a/_sources/operations/transformations/vgridshift.rst.txt b/_sources/operations/transformations/vgridshift.rst.txt new file mode 100644 index 00000000..5a31ae8f --- /dev/null +++ b/_sources/operations/transformations/vgridshift.rst.txt @@ -0,0 +1,133 @@ +.. _vgridshift: + +================================================================================ +Vertical grid shift +================================================================================ + +.. versionadded:: 5.0.0 + +Change Vertical datum change by grid shift + ++-----------------+--------------------------------------------------------------------+ +| **Alias** | vgridshift | ++-----------------+--------------------------------------------------------------------+ +| **Domain** | 3D and 4D | ++-----------------+--------------------------------------------------------------------+ +| **Input type** | Geodetic coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++-----------------+--------------------------------------------------------------------+ +| **Output type** | Geodetic coordinates (horizontal), meters (vertical), | +| | decimalyear (temporal) | ++-----------------+--------------------------------------------------------------------+ + +The vertical grid shift is done by offsetting the vertical input coordinates by +a specific amount determined by the loaded grids. The simplest use case of the +horizontal grid shift is applying a single grid. Here we change the vertical +reference from the ellipsoid to the global geoid model, EGM96:: + + +proj=vgridshift +grids=egm96_15.gtx + + +More than one grid can be loaded at the same time, for instance in the case where +a better geoid model than the global is available for a certain area. Here the +gridshift is set up so that the local DVR90 geoid model takes precedence over +the global model:: + + +proj=vgridshift +grids=@dvr90.gtx,egm96_15.gtx + +The ``@`` in the above example states that the grid is optional, in case the grid +is not found in the PROJ search path. The list of grids is prioritized so that +grids in the start of the list takes precedence over the grids in the back of the +list. + +PROJ supports the GTX file format for vertical grid corrections. Details +about all the format can be found in the GDAL documentation. GDAL both reads and +writes the format. Using GDAL for construction of new grids is recommended. + +Temporal gridshifting +################################################################################ +.. versionadded:: 5.1.0 + +By initializing the vertical gridshift operation with a central epoch, it can be +used as a step function applying the grid offsets only if a coordinate is +transformed from an epoch before grids central epoch to an epoch after. This is +handy in transformations where it is necessary to handle deformations caused by +seismic activity. + +The central epoch of the grid is controlled with :option:`+t_epoch` and the final +epoch of the coordinate is set with :option:`+t_final`. The observation epoch of +the coordinate is part of the coordinate tuple. + +Suppose we want to model the deformation of the 2008 earthquake in Iceland in +a transformation of data from 2005 to 2009:: + + echo 63.992 -21.014 10.0 2005.0 | cct +proj=vgridshift +grids=iceland2008.gtx +t_epoch=2008.4071 +t_final=2009.0 + 63.992 -21.014 10.11 2005.0 + +.. note:: + The timestamp of the resulting coordinate is still 2005.0. The observation + time is always kept unchanged as it would otherwise be impossible to do the + inverse transformation. + +Temporal gridshifting is especially powerful in transformation pipelines where +several gridshifts can be chained together, effectively acting as a series of +step functions that can be applied to a coordinate that is propagated through +time. In the following example we establish a pipeline that allows transformation +of coordinates from any given epoch up until the current date, applying only +those gridshifts that have central epochs between the observation epoch and +the final epoch:: + + +proj=pipeline +t_final=now + +step +proj=vgridshift +grids=earthquake_1.gtx +t_epoch=2010.421 + +step +proj=vgridshift +grids=earthquake_2.gtx +t_epoch=2013.853 + +step +proj=vgridshift +grids=earthquake_3.gtx +t_epoch=2017.713 + +.. note:: + The special epoch *now* is used when specifying the final epoch with + :option:`+t_final`. This results in coordinates being transformed to the + current date. Additionally, :option:`+t_final` is used as a + :ref:`global pipeline parameter<global-pipeline-parameter>`, which means + that it is applied to all the steps in the pipeline. + +In the above transformation, a coordinate with observation epoch 2009.32 would +be subject to all three gridshift steps in the pipeline. A coordinate with +observation epoch 2014.12 would only by offset by the last step in the pipeline. + +Parameters +################################################################################ + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +grids=<list> + + Comma-separated list of grids to load. If a grid is prefixed by an `@` the + grid is considered optional and PROJ will the not complain if the grid is + not available. + + Grids are expected to be in GTX format. + +Optional ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. include:: ../options/t_epoch.rst +.. versionadded:: 5.1.0 +.. include:: ../options/t_final.rst +.. versionadded:: 5.1.0 + +.. option:: +multiplier=<value> + + Specify the multiplier to apply to the grid value in the forward transformation + direction, such that: + + .. math:: + :label: multiplier_formula + + Z_{target} = Z_{source} + multiplier \times gridvalue + + The multiplier can be used to control whether the gridvalue should be added + or subtracted, and if unit conversion must be done (the multiplied gridvalue + must be expressed in metre). + + Note that the default is `-1.0` for historical reasons. +.. versionadded:: 5.2.0 diff --git a/_sources/operations/transformations/xyzgridshift.rst.txt b/_sources/operations/transformations/xyzgridshift.rst.txt new file mode 100644 index 00000000..e76b3a63 --- /dev/null +++ b/_sources/operations/transformations/xyzgridshift.rst.txt @@ -0,0 +1,93 @@ +.. _xyzgridshift: + +================================================================================ +Geocentric grid shift +================================================================================ + +.. versionadded:: 7.0.0 + +Geocentric translation using a grid shift + ++-----------------+-------------------------------------------------------------------+ +| **Alias** | xyzgridshift | ++-----------------+-------------------------------------------------------------------+ +| **Domain** | 3D | ++-----------------+-------------------------------------------------------------------+ +| **Input type** | Cartesian coordinates | ++-----------------+-------------------------------------------------------------------+ +| **Output type** | Cartesian coordinates | ++-----------------+-------------------------------------------------------------------+ + +Perform a geocentric translation by bilinear interpolation of dx, dy, dz +translation values from a grid. The grid is referenced against either the +2D geographic CRS corresponding to the input (or sometimes output) CRS. + +This method is described (in French) in :cite:`NTF_88` +and as EPSG operation method code 9655 in :cite:`IOGP2018` (§2.4.4.1.1 +France geocentric interpolation). + +The translation in the grids are added to the input coordinates in the forward direction, +and subtracted in the reverse direction. +By default (if ``grid_ref=input_crs``), in the forward direction, the input coordinates +are converted to their geographic equivalent to directly read and interpolate from +the grid. In the reverse direction, an iterative method is used to be able to find +the grid locations to read. +If ``grid_ref=output_crs`` is used, then the reverse strategy is applied: iterative +method in the forward direction, and direct read in the reverse direction. + +Example +------------------------------------------------------------------------------- + +NTF to RGF93 transformation using :file:`gr3df97a.tif` grid + +:: + + +proj=pipeline + +step +proj=push +v_3 + +step +proj=cart +ellps=clrk80ign + +step +proj=xyzgridshift +grids=gr3df97a.tif +grid_ref=output_crs +ellps=GRS80 + +step +proj=cart +ellps=GRS80 +inv + +step +proj=pop +v_3 + +Parameters +################################################################################ + +The ellipsoid parameters should be the ones consistent with ``grid_ref``. +They are used to perform a geocentric to geographic conversion to find the +translation parameters. + + +Required ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. include:: ../options/ellps.rst + +.. option:: +grids=<list> + + Comma-separated list of grids to load. If a grid is prefixed by an ``@`` the + grid is considered optional and PROJ will the not complain if the grid is + not available. + + Grids are expected to be in GeoTIFF format. If no metadata is provided, + the first, second and third samples are assumed to be the geocentric + translation along X, Y and Z axis respectively, in metres. + +Optional ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +.. option:: +grid_ref=input_crs/output_crs + + Specify in which CRS the grid is referenced to. The default value is + input_crs. That is the grid is referenced in the geographic CRS corresponding + to the input geocentric CRS. + + If output_crs is specified, the grid is referenced in the geographic CRS corresponding + to the output geocentric CRS. This is for example the case for the French + :file:`gr3df97a.tif` grid converting from NTF to RGF93, but referenced against RGF93. + Thus in the forward direction (NTF->RGF93), an iterative method is used to find + the appropriate shift. + +.. option:: +multiplier=<value> + + Specify the multiplier to apply to the grid values. Defaults to 1.0 + |