From 7ea472f54849f3230d5865884ef75c4edcda3fa8 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 28 Aug 2017 08:07:48 +0200 Subject: Major restructure and rewrite of docs. A section on "Using PROJ.4" has been introduced giving a introduction on how to use PROJ.4. The Parameters section has been removed and the content reworked into other sections. Additionally the order of chapters has been changed to provide a more complete and readable experience for users of the documentation. --- docs/source/usage/apps/cs2cs.rst | 143 ++++++++++++++++ docs/source/usage/apps/geod.rst | 170 +++++++++++++++++++ docs/source/usage/apps/index.rst | 20 +++ docs/source/usage/apps/proj.rst | 164 ++++++++++++++++++ docs/source/usage/index.rst | 20 +++ docs/source/usage/projections.rst | 163 ++++++++++++++++++ docs/source/usage/quickstart.rst | 68 ++++++++ docs/source/usage/resource_files.rst | 34 ++++ docs/source/usage/transformation.rst | 314 +++++++++++++++++++++++++++++++++++ 9 files changed, 1096 insertions(+) create mode 100644 docs/source/usage/apps/cs2cs.rst create mode 100644 docs/source/usage/apps/geod.rst create mode 100644 docs/source/usage/apps/index.rst create mode 100644 docs/source/usage/apps/proj.rst create mode 100644 docs/source/usage/index.rst create mode 100644 docs/source/usage/projections.rst create mode 100644 docs/source/usage/quickstart.rst create mode 100644 docs/source/usage/resource_files.rst create mode 100644 docs/source/usage/transformation.rst (limited to 'docs/source/usage') diff --git a/docs/source/usage/apps/cs2cs.rst b/docs/source/usage/apps/cs2cs.rst new file mode 100644 index 00000000..675379ef --- /dev/null +++ b/docs/source/usage/apps/cs2cs.rst @@ -0,0 +1,143 @@ +.. _cs2cs: + +================================================================================ +cs2cs +================================================================================ + + +``cs2cs`` performs transformation between the source and destination cartographic coordinate system +on a set of input points. The coordinate system transformation can include translation between +projected and geographic coordinates as well as the application of datum shifts. + + +Synopsis +******** + +:: + + cs2cs [ -eEfIlrstvwW [ args ] ] [ +opts[=arg] ] [ +to [+opts[=arg]] ] file[s] + +Description +*********** +The following control parameters can appear in any order: + +:: + + -I method to specify inverse translation, convert from +to coordinate system to the + primary coordinate system defined. + + -ta A specifies a character employed as the first character to denote a control line + to be passed through without processing. This option applicable to ascii input + only. (# is the default value). + + -e string + String is an arbitrary string to be output if an error is detected during data + transformations. The default value is: *\t*. Note that if the -b, -i or -o + options are employed, an error is returned as HUGE_VAL value for both + return values. + + -E causes the input coordinates to be copied to the output line prior to printing the + converted values. + + -l[p|P|=|e|u|d]id + List projection identifiers with -l, -lp or -lP (expanded) that can be selected + with +proj. -l=id gives expanded description of projection id. List + ellipsoid identifiers with -le, that can be selected with +ellps,-lu list of + cartesian to meter conversion factors that can be selected with +units or -ld + list of datums that can be selected with +datum. + + -r This options reverses the order of the expected input from longitude-latitude or + x-y to latitude-longitude or y-x. + + -s This options reverses the order of the output from x-y or longitude-latitude to + y-x or latitude-longitude. + + -f format + Format is a printf format string to control the form of the output values. For + inverse projections, the output will be in degrees when this option is employed. + If a format is specified for inverse projection the output data will be in deci- + mal degrees. The default format is "%.2f" for forward projection and DMS for + inverse. + + -[w|W]n + N is the number of significant fractional digits to employ for seconds output (when + the option is not specified, -w3 is assumed). When -W is employed the fields + will be constant width and with leading zeroes. + + -v causes a listing of cartographic control parameters tested for and used by the + program to be printed prior to input data. + +The ``+args`` run-line arguments are associated with cartographic parameters. Usage varies with +projection and for a complete description consult the `projection pages <../projections/index.html>`_ + +The cs2cs program requires two coordinate system +definitions. The first (or primary is defined +based on all projection parameters not appearing +after the +to argument. All projection parameters +appearing after the ``+to`` argument are considered +the definition of the second coordinate +system. If there is no second coordinate system +defined, a geographic coordinate system based on +the datum and ellipsoid of the source coordinate +system is assumed. Note that the source and +destination coordinate system can both be projections, +both be geographic, or one of each and +may have the same or different datums. + +Additional projection control parameters may be +contained in two auxiliary control files: the +first is optionally referenced with the +``+init=file:id`` and the second is always processed +after the name of the projection has been established +from either the run-line or the contents +of ``+init`` file. The environment parameter +PROJ_LIB establishes the default directory for a +file reference without an absolute path. This +is also used for supporting files like datum +shift files. + +One or more files (processed in left to right +order) specify the source of data to be transformed. +A - will specify the location of processing standard input. +If no files are specified, the input is assumed to be from stdin. +For input data the two data values must be in +the first two white space separated fields and +when both input and output are ASCII all trailing +portions of the input line are appended to +the output line. + +Input geographic data (longitude and latitude) +must be in DMS or decimal degrees format and +input cartesian data must be in units consistent +with the ellipsoid major axis or sphere radius +units. Output geographic coordinates will normally +be in DMS format (use ``-f %.12f`` for decimal +degrees with 12 decimal places), while projected +(cartesian) coordinates will be in linear +(meter, feet) units. + + +Example +******* +The following script + +:: + + cs2cs +proj=latlong +datum=NAD83 + +to +proj=utm +zone=10 +datum=NAD27 + -r < [ -afFIlptwW [ args ] ] [ +args ] file[s] + + invgeod +ellps= [ -afFIlptwW [ args ] ] [ +args ] file[s] + +Description +*********** +The following command-line options can appear in any order: + +:: + + -I Specifies that the inverse geodesic computation is to be + performed. May be used with execution of geod as an + alternative to invgeod execution. + + -a Latitude and longitudes of the initial and terminal + points, forward and back azimuths and distance are output. + + -ta A specifies a character employed as the first character + to denote a control line to be passed through without + processing. + + -le Gives a listing of all the ellipsoids that may be + selected with the +ellps= option. + + -lu Gives a listing of all the units that may be selected + with the +units= option. + + -[f|F] format + Format is a printf format string to control the output + form of the geographic coordinate values (f) or distance + value (F). The default mode is DMS for geographic + coordinates and "%.3f" for distance. + + -[w|W]n + N is the number of significant fractional digits to + employ for seconds output (when the option is not + specified, -w3 is assumed). When -W is employed the fields + will be constant width with leading zeroes. + + -p This option causes the azimuthal values to be output as + unsigned DMS numbers between 0 and 360 degrees. Also + note -f. + +The ``+args`` command-line options are associated with geodetic +parameters for specifying the ellipsoidal or sphere to use. +See ``proj`` `documentation `_ for full list of these parameters and +controls. The options are processed in left to right order +from the command line. Reentry of an option is ignored with +the first occurrence assumed to be the desired value. + +One or more files (processed in left to right order) specify +the source of data to be transformed. A - will specify the +location of processing standard input. If no files are specified, +the input is assumed to be from stdin. + +For direct determinations input data must be in latitude, longitude, +azimuth and distance order and output will be latitude, +longitude and back azimuth of the terminus point. Latitude, +longitude of the initial and terminus point are input for the +inverse mode and respective forward and back azimuth from the +initial and terminus points are output along with the distance +between the points. + +Input geographic coordinates (latitude and longitude) and +azimuthal data must be in decimal degrees or DMS format and +input distance data must be in units consistent with the ellipsoid +major axis or sphere radius units. The latitude must lie +in the range [-90d,90d]. Output geographic coordinates will be +in DMS (if the -f switch is not employed) to 0.001" with trailing, +zero-valued minute-second fields deleted. Output distance +data will be in the same units as the ellipsoid or sphere +radius. + +The Earth's ellipsoidal figure may be selected in the same manner +as program ``proj`` by using ``+ellps=``, ``+a=``, ``+es=``, etc. + +Geod may also be used to determine intermediate points along +either a geodesic line between two points or along an arc of +specified distance from a geographic point. In both cases an +initial point must be specified with ``+lat_1=lat`` and ``+lon_1=lon`` +parameters and either a terminus point ``+lat_2=lat`` and +``+lon_2=lon`` or a distance and azimuth from the initial point +with ``+S=distance`` and ``+A=azimuth`` must be specified. + +If points along a geodesic are to be determined then either +``+n_S=integer`` specifying the number of intermediate points +and/or ``+del_S=distance`` specifying the incremental distance +between points must be specified. + +To determine points along an arc equidistant from the initial +point both ``+del_A=angle`` and ``+n_A=integer`` must be specified +which determine the respective angular increments and number of +points to be determined. + +Examples +******** + +The following script determines the geodesic azimuths and distance in U.S. +statute miles from Boston, MA, to Portland, OR: + +:: + + geod +ellps=clrk66 <`_ + +#. `C. F. F. Karney, Algorithms for Geodesics, J. Geodesy 87, 43-55 (2013) `_. + `Addendum `_ + +#. `The online geodesic bibliography `_ + + diff --git a/docs/source/usage/apps/index.rst b/docs/source/usage/apps/index.rst new file mode 100644 index 00000000..73c5c63c --- /dev/null +++ b/docs/source/usage/apps/index.rst @@ -0,0 +1,20 @@ +.. _apps: + +================================================================================ +Applications +================================================================================ + +Bundled with PROJ.4 comes a set of small command line utilities. The +``proj`` program is limited to converting between geographic and projection +coordinates within one datum. The ``cs2cs`` program operates similarly, but +allows translation between any pair of definable coordinate systems, including +support for basic datum translation. The ``geod`` program provides the ability +to do geodesic (great circle) computations. + +.. toctree:: + :maxdepth: 1 + + proj + cs2cs + geod + diff --git a/docs/source/usage/apps/proj.rst b/docs/source/usage/apps/proj.rst new file mode 100644 index 00000000..424a0fef --- /dev/null +++ b/docs/source/usage/apps/proj.rst @@ -0,0 +1,164 @@ +.. _proj: + +================================================================================ +proj +================================================================================ + + +.. Index:: proj + +``proj`` and ``invproj`` perform respective forward and inverse transformation of cartographic data to +or from cartesian data with a wide range of selectable projection functions. + + +Synopsis +******** +:: + + proj [ -bcCeEfiIlmorsStTvVwW [ args ] ] [ +args ] file[s] + invproj [ -bcCeEfiIlmorsStTwW [ args ] ] [ +args ] file[s] + + +Description +*********** +The following control parameters can appear in any order + +:: + + -b Special option for binary coordinate data input and output through standard input + and standard output. Data is assumed to be in system type double floating point + words. This option is to be used when proj is a son process and allows bypassing + formatting operations. + + -i Selects binary input only (see -b option). + + -C Check. Invoke all built in self tests and report. Get more verbose report by + preceding with the -V option). + + -I alternate method to specify inverse projection. Redundant when used with invproj. + + -o Selects binary output only (see -b option). + + -ta A specifies a character employed as the first character to denote a control line + to be passed through without processing. This option applicable to ascii input + only. (# is the default value). + + -e string + String is an arbitrary string to be output if an error is detected during data + transformations. The default value is: *\t*. Note that if the -b, -i or -o + options are employed, an error is returned as HUGE_VAL value for both return values. + + -E causes the input coordinates to be copied to the output line prior to printing the + converted values. + + + -l[p|P|=|e|u|d]id + List projection identifiers with -l, -lp or -lP (expanded) that can be selected + with +proj. -l=id gives expanded description of projection id. List + ellipsoid identifiers with -le, that can be selected with +ellps, -lu list of + cartesian to meter conversion factors that can be selected with +units or -ld + list of datums that can be selected with +datum. + + -r This options reverses the order of the expected input from longitude-latitude or + x-y to latitude-longitude or y-x. + + -s This options reverses the order of the output from x-y or longitude-latitude to + y-x or latitude-longitude. + + -S Causes estimation of meridinal and parallel scale factors, area scale factor and + angular distortion, and maximum and minimum scale factors to be listed between <> + for each input point. For conformal projections meridinal and parallel scales + factors will be equal and angular distortion zero. Equal area projections will + have an area factor of 1. + + -m mult + The cartesian data may be scaled by the mult parameter. When processing data in + a forward projection mode the cartesian output values are multiplied by mult + otherwise the input cartesian values are divided by mult before inverse projection. + If the first two characters of mult are 1/ or 1: then the reciprocal value of mult + is employed. + + -f format + Format is a printf format string to control the form of the output values. For + inverse projections, the output will be in degrees when this option is employed. + The default format is "%.2f" for forward projection and DMS for inverse. + + -[w|W]n + N is the number of significant fractional digits to employ for seconds output (when + the option is not specified, -w3 is assumed). When -W is employed the fields + will be constant width and with leading zeroes. + + -v causes a listing of cartographic control parameters tested for and used by the + program to be printed prior to input data. Should not be used with the -T + option. + + -V This option causes an expanded annotated listing of the characteristics of the + projected point. -v is implied with this option. + + -T ulow,uhi,vlow,vhi,res[,umax,vmax] + This option creates a set of bivariate Chebyshev polynomial coefficients that + approximate the selected cartographic projection on stdout. The values low and + hi denote the range of the input where the u or v prefixes apply to respective + longitude-x or latitude-y depending upon whether a forward or inverse projection + is selected. Res is an integer number specifying the power of 10 precision of + the approximation. For example, a res of -3 specifies an approximation with an + accuracy better than .001. Umax, and vmax specify maximum degree of the polynomials + (default: 15). + + +The ``+args`` run-line arguments are associated with cartographic parameters. Usage varies with +projection and for a complete description consult the `projection pages <../projections/index.html>`_ + +Additional projection control parameters may be +contained in two auxiliary control files: the +first is optionally referenced with the +``+init=file:id`` and the second is always processed +after the name of the projection has been established +from either the run-line or the contents +of +init file. The environment parameter +``PROJ_LIB`` establishes the default directory for a +file reference without an absolute path. This +is also used for supporting files like datum +shift files. + +One or more files (processed in left to right +order) specify the source of data to be trans- +formed. A - will specify the location of pro- +cessing standard input. If no files are speci- +fied, the input is assumed to be from stdin. +For ASCII input data the two data values must be +in the first two white space separated fields +and when both input and output are ASCII all +trailing portions of the input line are appended +to the output line. + +Input geographic data (longitude and latitude) +must be in DMS format and input cartesian data +must be in units consistent with the ellipsoid +major axis or sphere radius units. Output geo- +graphic coordinates will be in DMS (if the -w +switch is not employed) and precise to 0.001" +with trailing, zero-valued minute-second fields +deleted. + +Example +******* +The following script + +:: + + proj +proj=utm +lon_0=112w +ellps=clrk66 + -r <` available in the library. This section +is devoted to the generic parameters that can be used on any projection in the +PROJ.4 library. + +Below is a list of PROJ.4 parameters which can be applied to most coordinate +system definitions. This table does not attempt to describe the parameters +particular to particular projection types. These can be found on the pages +documenting the individual :doc:`projections<../projections/index>`. + + ========== ================================================================ + Parameter Description + ========== ================================================================ + +a Semimajor radius of the ellipsoid axis + +axis Axis orientation + +b Semiminor radius of the ellipsoid axis + +ellps Ellipsoid name (see ``proj -le``) + +k Scaling factor (deprecated) + +k_0 Scaling factor + +lat_0 Latitude of origin + +lon_0 Central meridian + +lon_wrap Center longitude to use for wrapping (see below) + +no_defs Don't use the /usr/share/proj/proj_def.dat defaults file + +over Allow longitude output outside -180 to 180 range, disables + wrapping (see below) + +pm Alternate prime meridian (typically a city name, see below) + +proj Projection name (see ``proj -l``) + +units meters, US survey feet, etc. + +vunits vertical units. + +x_0 False easting + +y_0 False northing + ========== ================================================================ + +In the sections below most of the parameters are explained in details. + +Units ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Horizontal units can be specified using the ``+units`` keyword with a symbolic +name for a unit (ie. ``us-ft``). Alternatively the translation to meters can be +specified with the ``+to_meter`` keyword (ie. 0.304800609601219 for US feet). The +``-lu`` argument to ``cs2cs`` or ``proj`` can be used to list symbolic unit names. +The default unit for projected coordinates is the meter. +A few special projections deviate from this behaviour, most notably the +latlong pseudo-projection that returns degrees. + +Vertical (Z) units can be specified using the ``+vunits`` keyword with a +symbolic name for a unit (ie. ``us-ft``). Alternatively the translation to +meters can be specified with the ``+vto_meter`` keyword (ie. 0.304800609601219 +for US feet). The ``-lu`` argument to ``cs2cs`` or ``proj`` can be used to list +symbolic unit names. If no vertical units are specified, the vertical units will +default to be the same as the horizontal coordinates. + +.. note:: + ``proj`` do not handle vertical units at all and hence the ``+vto_meter`` + argument will be ignored. + +Scaling of output units can be done by applying the ``+k_0`` argument. The +returned coordinates are scaled by the value assigned with the ``+k_0`` +parameter. + +False Easting/Northing ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Virtually all coordinate systems allow for the presence of a false easting +(``+x_0``) and northing (``+y_0``). Note that these values are always expressed in +meters even if the coordinate system is some other units. Some coordinate +systems (such as UTM) have implicit false easting and northing values. + +Longitude Wrapping ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +By default PROJ.4 wraps output longitudes in the range -180 to 180. The ``+over`` +switch can be used to disable the default wrapping which is done at a low level +in ``pj_inv()``. This is particularly useful with projections like the +:doc:`equidistant cylindrical<../projections/eqc>` where it would be desirable for +X values past -20000000 (roughly) to continue past -180 instead of wrapping to +180. + +The ``+lon_wrap`` option can be used to provide an alternative means of doing +longitude wrapping within ``pj_transform()``. The argument to this option is a +center longitude. So ``+lon_wrap=180`` means wrap longitudes in the range 0 to +360. Note that ``+over`` does **not** disable ``+lon_wrap``. + +Prime Meridian ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +A prime meridian may be declared indicating the offset between the prime +meridian of the declared coordinate system and that of greenwich. A prime +meridian is clared using the "pm" parameter, and may be assigned a symbolic +name, or the longitude of the alternative prime meridian relative to greenwich. + +Currently prime meridian declarations are only utilized by the +``pj_transform()`` API call, not the ``pj_inv()`` and ``pj_fwd()`` calls. +Consequently the user utility ``cs2cs`` does honour prime meridians but the +``proj`` user utility ignores them. + +The following predeclared prime meridian names are supported. These can be +listed using with ``cs2cs -lm``. + + =========== ================ + Meridian Longitude + =========== ================ + greenwich 0dE + lisbon 9d07'54.862"W + paris 2d20'14.025"E + bogota 74d04'51.3"E + madrid 3d41'16.48"W + rome 12d27'8.4"E + bern 7d26'22.5"E + jakarta 106d48'27.79"E + ferro 17d40'W + brussels 4d22'4.71"E + stockholm 18d3'29.8"E + athens 23d42'58.815"E + oslo 10d43'22.5"E + =========== ================ + +Example of use. The location ``long=0``, ``lat=0`` in the greenwich based lat/long +coordinates is translated to lat/long coordinates with Madrid as the prime +meridian. + +:: + + cs2cs +proj=latlong +datum=WGS84 +to +proj=latlong +datum=WGS84 +pm=madrid + 0 0 (input) + 3d41'16.48"E 0dN 0.000 (output) + + +Axis orientation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +Starting in PROJ 4.8.0, the +axis argument can be used to control the axis +orientation of the coordinate system. The default orientation is "easting, +northing, up" but directions can be flipped, or axes flipped using combinations +of the axes in the +axis switch. The values are: + +* "e" - Easting +* "w" - Westing +* "n" - Northing +* "s" - Southing +* "u" - Up +* "d" - Down + +They can be combined in +axis in forms like: + +* ``+axis=enu`` - the default easting, northing, elevation. +* ``+axis=neu`` - northing, easting, up - useful for "lat/long" geographic + coordinates, or south orientated transverse mercator. +* ``+axis=wnu`` - westing, northing, up - some planetary coordinate systems + have "west positive" coordinate systems + +.. note:: + + The ``+axis`` argument does not work with the ``proj`` command line + utility. + + diff --git a/docs/source/usage/quickstart.rst b/docs/source/usage/quickstart.rst new file mode 100644 index 00000000..d162ace9 --- /dev/null +++ b/docs/source/usage/quickstart.rst @@ -0,0 +1,68 @@ +.. _quickstart: + +================================================================================ +Quick start +================================================================================ + +Coordinate transformations are defined by, what in PROJ.4 terminology is +known as, "proj-strings". A proj-string describes any transformation regardless of +how simple or complicated it might be. The simplest case is projection of geodetic +coordinates. This section focuses on the simpler cases and introduces the basic +anatomy of the proj-string. The complex cases are discussed in +:doc:`transformation`. + +A proj-strings holds the parameters of a given coordinate transformation, e.g. + +:: + + +proj=merc +lat_ts=56.5 +ellps=GRS80 + +I.e. a proj-string consists of a projection specifier, ``+proj``, a number of +parameters that applies to the projection and, if needed, a description of a +datum shift. In the example above geodetic coordinates are transformed to +projected space with the :doc:`Mercator projection<../projections/merc>` with +the latitude of true scale at 56.5 degrees north on the GRS80 ellipsoid. Every +projection in PROJ.4 is identified by a shorthand such as ``merc`` in the above +example. + +By using the above projection definition as parameters for the command line +utility ``proj`` we can convert the geodetic coordinates to projected space: + +:: + + $ proj +proj=merc +lat_ts=56.5 +ellps=GRS80 + +If called as above ``proj`` will be in interactive mode, letting you type the +input data manually and getting a responce presented on screen. ``proj`` +works as any UNIX filter though, which means that you can also pipe data to +the utility, for instance by using the ``echo`` command: + +:: + + $ echo 55.2 12.2 | proj +proj=merc +lat_ts=56.5 +ellps=GRS80 + 3399483.80 752085.60 + + +PROJ.4 also comes bundled with the ``cs2cs`` utility which is used to transform +from onecoordinate reference system to another. Say we want to convert +the above Mercator coordinates to UTM, we can do that with ``cs2cs``: + +:: + + $ echo 3399483.80 752085.60 | cs2cs +proj=merc +lat_ts=56.5 +ellps=GRS80 +to +proj=utm +zone=32 + 6103992.36 1924052.47 0.00 + +Notice the ``+to`` parameter that seperates the source and destination +projection definitions. + +If you happen to know the EPSG identifiers for the two cordinates reference +systems you are transforming between you can use those with ``cs2cs``: + +:: + + $ echo 56 12 | cs2cs +init=epsg:4326 +to +init=epsg:25832 + 231950.54 1920310.71 0.00 + +In the above example we transform geodetic coordinates in the WGS84 reference +frame to UTM zone 32N coordinates in the ETRS89 reference frame. +UTM coordinates diff --git a/docs/source/usage/resource_files.rst b/docs/source/usage/resource_files.rst new file mode 100644 index 00000000..0881cb7c --- /dev/null +++ b/docs/source/usage/resource_files.rst @@ -0,0 +1,34 @@ +.. _resource_files: + +================================================================================ +Resource files +================================================================================ + +Init files +------------------------------------------------------------------------------- + +epsg ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +nad27 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +nad83 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +ITRF ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +IGNF ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +GL27 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +esri ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + +Defaults file +------------------------------------------------------------------------------- + diff --git a/docs/source/usage/transformation.rst b/docs/source/usage/transformation.rst new file mode 100644 index 00000000..98e607a7 --- /dev/null +++ b/docs/source/usage/transformation.rst @@ -0,0 +1,314 @@ +.. _transformation: + +================================================================================ +Geodetic transformation +================================================================================ + +PROJ.4 can do everything from the most simple projection to very complex +transformations across many reference frames. While originally developed as a +tool for cartographic projections, PROJ.4 has over time evolved into a powerfull +generic coordinate transformation engine that makes it possible to do both +large scale cartographic projections as well as coordinate transformation at a +geodetic high precision level. This chapter delves into the details of how +geodetec transformations of varying complexity can be performed. + +In PROJ.4, two frameworks for geodetic transformations exists, the *cs2cs* +framework and the *transformation pipelines* framework. The first is the original, +and limited, framework for doing geodetic transforms in PROJ.4 The latter is a +newer addition that aims to be a more complete transformation framework. Both are +described in the sections below. Large portions of the text are based on +[EversKnudsen2017]_. + +Before describing the details of the two frameworks, let us first remark that +most cases of geodetic transformations can be expressed as a series of elementary +operations, the output of one operation being the input of the next. E.g. when +going from UTM zone 32, datum ED50, to UTM zone 32, datum ETRS89, one must, in the +simplest case, go through 5 steps: + +1. Back-project the UTM coordinates to geographic coordinates +2. Convert the geographic coordinates to 3D cartesian geocentric coordinates +3. Apply a Helmert transformation from ED50 to ETRS89 +4. Convert back from cartesian to geographic coordinates +5. Finally project the geographic coordinates to UTM zone 32 planar coordinates. + + +Transformation pipelines ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +The homology between the above steps and a Unix shell style pipeline is evident. +It is there the main architectural inspiration behind the transformation pipeline +framework. The pipeline framework is realized by utilizing a special "projection", +that takes as its user supplied arguments, a series of elementary operations, +which it strings together in order to implement the full transformation needed. +Additionally, a number of elementary geodetic operations, including Helmert +transformations, general high order polynomial shifts and the Molodensky +transformation are available as part of the pipeline framework. +In anticipation of upcoming support for full time-varying transformations, we +also introduce a 4D spatiotemporal data type, and a programming interface +(API) for handling this. + +The Molodensky transformation converts directly from geodetic coordinates +in one datum, to geodetic coordinates in another datum, while the (typically more +accurate) Helmert transformation converts from 3D cartesian to 3D cartesian +coordinates. So when using the Helmert transformation one typically needs to do an +initial conversion from geodetic to cartesian coordinates, and a final conversion +the other way round, to arrive at the desired result. Fortunately, this three-step +compound transformation has the attractive characteristic that each step depends +only on the output of the immediately preceding step. Hence, we can build a +geodetic-to-geodetic Helmert transformation by tying together the outputs and inputs +of 3 steps (geodetic-to-cartesian → Helmert → cartesian-to-geodetic), pipeline style. +The pipeline driver, makes this kind of chained transformations possible. +The implementation is compact, consisting of just one pseudo-projection, called +``pipeline``, which takes as its arguments strings of elementary projections +(note: "projection" is the, slightly misleading, PROJ.4 term used for any kind of +transformation). +The pipeline pseudo projection is supplemented by a number of elementary +transformations, all in all providing a framework for building high accuracy +solutions for a wide spectrum of geodetic tasks. + + +As a first example, let us take a look at the iconic +*geodetic → Cartesian → Helmert → geodetic* case (steps 2 to 4 in the example in +the itroduction). In PROJ.4 it can be implemented as + +:: + + proj=pipeline + step proj=cart ellps=intl + step proj=helmert + x=-81.0703 y=-89.3603 z=-115.7526 + rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645 + step proj=cart inv ellps=GRS80 + +The pipeline can be expanded at both ends to accommodate whatever coordinate type +is needed for input and output: In the example below, we transform from the +deprecated Danish System 45, a 2D system with some tension in the original defining +network, to UTM zone 33, ETRS89. The tension is reduced using a polynomial +transformation (the init=./s45b... step, s45b.pol is a file containing the +polynomial coefficients), taking the S45 coordinates to a technical coordinate +system (TC32), defined to represent "UTM zone 32 coordinates, as they would look if +the Helmert transformation between ED50 and ETRS89 was perfect". The TC32 +coordinates are then converted back to geodetic(ED50) coordinates, using an +inverse UTM projection, further to cartesian(ED50), then to cartesian(ETRS89), +using the relevant Helmert transformation, and back to geodetic(ETRS89), before +finally being projected onto the UTM zone 33, ETRS89 system. All in all a 6 step +pipeline, implementing a transformation with centimeter level accuracy from a +deprecated system with decimeter level tensions. + +:: + + proj=pipeline + step init=./s45b.pol:s45b_tc32 + step proj=utm inv ellps=intl zone=32 + step proj=cart ellps=intl + step proj=helmert + x=-81.0703 y=-89.3603 z=-115.7526 + rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645 + step proj=cart inv ellps=GRS80 + step proj=utm ellps=GRS80 zone=33 + +With the pipeline framework spatiotemporal transformation is possible. This is +possible by leveraging the time dimension in PROJ.4 that enables 4D coordinates +(three spatial components and one temporal component) to be passed through a +transformation pipeline. In the example below a transformation from ITRF93 to +ITRF2000 is defined. The temporal component is given as GPS weeks in the input +data, but the 14-parameter Helmert transform expects temporal units in decimalyears. +Hence the first step in the pipeline is the unitconvert pseudo-projection that makes +sure the correct units are passed along to the Helmert transform. +Most parameters of the Helmert transform are taken from [AltamimiEtAl2002]_, +except the epoch which is the epoch of the transformation. The default setting is to +use “coordinate frame” convention of the Helmert transform, but “position vector” +convention can also be used. The last step in the pipeline is converting the +coordinate timestamps back to GPS weeks. + +:: + + proj=pipeline + step proj=unitconvert t_in=gps_week t_out=decimalyear + step proj=helmert + x=0.0127 y=0.0065 z=-0.0209 s=0.00195 + rx=0.00039 ry=-0.00080 rz=0.00114 + dx=-0.0029 dy=-0.0002 dz=-0.0006 ds=0.00001 + drx=0.00011 dry=0.00019 drz=-0.00007 + epoch=1988.0 + step proj=unitconvert t_in=decimalyear t_out=gps_week + + +cs2cs paradigm ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + ============ ============================================================== + Parameter Description + ============ ============================================================== + +datum Datum name (see ``proj -ld``) + +geoidgrids Filename of GTX grid file to use for vertical datum transforms + +nadgrids Filename of NTv2 grid file to use for datum transforms + +towgs84 3 or 7 term datum transform parameters + +to_meter Multiplier to convert map units to 1.0m + +vto_meter Vertical conversion to meters + ============ ============================================================== + +The *cs2cs* framework delivers a subset of the geodetic transformations available +with the *pipeline* framework. Coordinate transformations done in this framework +are transformed in a two-step process with WGS84 as a pivot datum That is, the +input coordinates are transformed to WGS84 geodetic coordinates and then transformed +from WGS84 coordinates to the specified output coordinate reference system, by +utilizing either the Helmert transform, datum shift grids or a combination of both. +Datum shifts can be described in a proj-string with the parameters ``+towgs84``, +``+nadgrids`` and ``+geoidgrids``. +An inverse transform exists for all three and is applied if +specified in the input proj-string. The most common is ``+towgs84``, which is used to +define a 3- or 7-parameter Helmert shift from the input reference frame to WGS84. +Exactly which realization of WGS84 is not specified, hence a fair amount of +uncertainty is introduced in this step of the transformation. With the +nadgrids +parameter a non-lineaer planar correction derived from interpolation in a +correction grid can be applied. Originally this was implemented as a means to +transform coordinates between the american datums NAD27 and NAD83, but corrections +can be applied for any datum for which a correction grid exists. The inverse +transform for the horizontal grid shift is "dumb", in the sense that the +correction grid is applied verbatim without taking into account that the inverse +operation is non-linear. Similar to the horizontal grid correction, ``+geoidgrids`` +can be used to perform grid corrections in the vertical component. +Both grid correction methods allow inclusion of more than one grid in the same +transformation + +In contrast to the *transformation pipeline* framework, transformations with the +*cs2cs* framework are expressed as two separate proj-strings. One proj-string *to* +WGS84 and one *from* WGS84. Together they form the mapping from the source +coordinate reference system to the destination coordinate reference system. +When used with the ``cs2cs`` the source and destination CRS's are separated by the +special ``+to`` parameter. + +The following example demonstrates converting from the Greek GGRS87 datum +to WGS84 with the ``+towgs84`` parameter. + +:: + + cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62 + +to +proj=latlong +datum=WGS84 + 20 35 + 20d0'5.467"E 35d0'9.575"N 8.570 + +The EPSG database provides this example for transforming from WGS72 to WGS84 +using an approximated 7 parameter transformation. + +:: + + cs2cs +proj=latlong +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.219 \ + +to +proj=latlong +datum=WGS84 + 4 55 + 4d0'0.554"E 55d0'0.09"N 3.223 + + +Grid Based Datum Adjustments ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +In many places (notably North America and Austrialia) national geodetic +organizations provide grid shift files for converting between different datums, +such as NAD27 to NAD83. These grid shift files include a shift to be applied +at each grid location. Actually grid shifts are normally computed based on an +interpolation between the containing four grid points. + +PROJ.4 supports use of grid files for shifting between various reference frames. +The grid shift table formats are ctable (the binary format produced by the PROJ.4 +``nad2bin`` program), NTv1 (the old Canadian format), and NTv2 (``.gsb`` - the new +Canadian and Australian format). + +The text in this section is based on the *cs2cs* framework. Gridshifting is off +course also possible with the *pipeline* framework. The major difference between the +two is that the *cs2cs* framework is limited to grid mappings to WGS84, whereas with +*transformation pipelines* it is possible to perform grid shifts between any two +reference frames, as long as a grid exists. + +Use of grid shifts with ``cs2cs`` is specified using the ``+nadgrids`` +keyword in a coordinate system definition. For example: + +:: + + % cs2cs +proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat \ + +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF + -111 50 + EOF + 111d0'2.952"W 50d0'0.111"N 0.000 + +In this case the ``/usr/local/share/proj/ntv1_can.dat`` grid shift file was +loaded, and used to get a grid shift value for the selected point. + +It is possible to list multiple grid shift files, in which case each will be +tried in turn till one is found that contains the point being transformed. + +:: + + cs2cs +proj=latlong +ellps=clrk66 \ + +nadgrids=conus,alaska,hawaii,stgeorge,stlrnc,stpaul \ + +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF + -111 44 + EOF + 111d0'2.788"W 43d59'59.725"N 0.000 + + +Skipping Missing Grids +................................................................................ + +The special prefix ``@`` may be prefixed to a grid to make it optional. If it +not found, the search will continue to the next grid. Normally any grid not +found will cause an error. For instance, the following would use the +``ntv2_0.gsb`` file if available (see :ref:`nonfreegrids`), otherwise it would +fallback to using the ``ntv1_can.dat`` file. + +:: + + cs2cs +proj=latlong +ellps=clrk66 +nadgrids=@ntv2_0.gsb,ntv1_can.dat \ + +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF + -111 50 + EOF + 111d0'3.006"W 50d0'0.103"N 0.000 + +The null Grid +................................................................................ + +A special ``null`` grid shift file is shift with releases after 4.4.6 (not +inclusive). This file provides a zero shift for the whole world. It may be +listed at the end of a nadgrids file list if you want a zero shift to be +applied to points outside the valid region of all the other grids. Normally if +no grid is found that contains the point to be transformed an error will occur. + +:: + + cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null \ + +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF + -111 45 + EOF + 111d0'3.006"W 50d0'0.103"N 0.000 + + cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null \ + +to +proj=latlong +ellps=GRS80 +datum=NAD83 << EOF + -111 44 + -111 55 + EOF + 111d0'2.788"W 43d59'59.725"N 0.000 + 111dW 55dN 0.000 + +For more information see the chapter on :ref:`Grids`. + +Caveats +................................................................................ + +* Where grids overlap (such as conus and ntv1_can.dat for instance) the first + found for a point will be used regardless of whether it is appropriate or + not. So, for instance, ``+nadgrids=ntv1_can.dat``,conus would result in + the Canadian data being used for some areas in the northern United States + even though the conus data is the approved data to use for the area. + Careful selection of files and file order is necessary. In some cases + border spanning datasets may need to be pre-segmented into Canadian and + American points so they can be properly grid shifted +* There are additional grids for shifting between NAD83 and various HPGN + versions of the NAD83 datum. Use of these haven't been tried recently so + you may encounter problems. The FL.lla, WO.lla, MD.lla, TN.lla and WI.lla + are examples of high precision grid shifts. Take care! +* Additional detail on the grid shift being applied can be found by setting + the PROJ_DEBUG environment variable to a value. This will result in output + to stderr on what grid is used to shift points, the bounds of the various + grids loaded and so forth +* The *cs2cs* framework always assumes that grids contain a shift **to** NAD83 (essentially + WGS84). Other types of grids can be used with the *pipeline* framework. -- cgit v1.2.3