Merge pull request #1385 from kbevers/doc-axis-ordering

Doc axis ordering

3 files changed, 89 insertions, 9 deletions

diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index 4052ff82..0066c3a2 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -221,10 +221,16 @@ Coordinate transformation
         3. of length one, i.e. a constant, which will be treated as a fully
            populated array of that constant value
 
+    .. note:: Even though he coordinate components are named :c:data:`x`, :c:data:`y`,
+              :c:data:`z` and :c:data:`t`, axis ordering of the to and from CRS
+              is respected. Transformations exhibit the same behaviour
+              as if they were gathered in a :c:type:`PJ_COORD` struct.
+
+
     The strides, :c:data:`sx`, :c:data:`sy`, :c:data:`sz`, :c:data:`st`,
     represent the step length, in bytes, between
     consecutive elements of the corresponding array. This makes it possible for
-    :c:func:`proj_transform` to handle transformation of a large class of application
+    :c:func:`proj_trans_generic` to handle transformation of a large class of application
     specific data structures, without necessarily understanding the data structure
     format, as in:
 
@@ -250,21 +256,22 @@ Coordinate transformation
             0, 0                          /*  and the time is the constant 0.00 s */
         );
 
-    This is similar to the inner workings of the deprecated pj_transform function, but the
-    stride functionality has been generalized to work for any size of basic unit,
-    not just a fixed number of doubles.
+    This is similar to the inner workings of the deprecated :c:func:`pj_transform`
+    function, but the stride functionality has been generalized to work for any
+    size of basic unit, not just a fixed number of doubles.
 
     In most cases, the stride will be identical for x, y, z, and t, since they will
-    typically be either individual arrays (stride = sizeof(double)), or strided
-    views into an array of application specific data structures (stride = sizeof (...)).
+    typically be either individual arrays (``stride = sizeof(double)``), or strided
+    views into an array of application specific data structures (``stride = sizeof (...)``).
 
     But in order to support cases where :c:data:`x`, :c:data:`y`, :c:data:`z`,
     and :c:data:`t` come from heterogeneous sources, individual strides,
     :c:data:`sx`, :c:data:`sy`, :c:data:`sz`, :c:data:`st`, are used.
 
-    .. note:: Since :c:func:`proj_transform` does its work *in place*, this means that even the
-              supposedly constants (i.e. length 1 arrays) will return from the call in altered
-              state. Hence, remember to reinitialize between repeated calls.
+    .. note:: Since :c:func:`proj_trans_generic` does its work *in place*,
+              this means that even the supposedly constants (i.e. length 1 arrays)
+              will return from the call in altered state. Hence, remember to
+              reinitialize between repeated calls.
 
     :param PJ* P: Transformation object
     :param `direction`: Transformation direction
diff --git a/docs/source/faq.rst b/docs/source/faq.rst
index d6f118ce..e6a2e54a 100644
--- a/docs/source/faq.rst
+++ b/docs/source/faq.rst
@@ -96,3 +96,63 @@ erroneous transformations.
 
 For compatibility reasons PROJ supports several WKT dialects
 (see :option:`projinfo -o`). If possible WKT2 should be used.
+
+Why is the axis ordering in PROJ not consistent?
+--------------------------------------------------------------------------------
+
+PROJ respects the axis ordering as it was defined by the authority in charge of
+a given coordinate reference system. This is in accordance to the ISO19111 standard
+:cite:`ISO19111`. Unfortunately most GIS software on the market doesn't follow this
+standard. Before version 6, PROJ did not respect the standard either. This causes
+some problems while the rest of the industry conforms to the standard. PROJ intends
+to spearhead this effort, hopefully setting a good example for the rest of the
+geospatial industry.
+
+Customarily in GIS the first component in a coordinate tuple has been aligned with
+the east/west direction and the second component with the north/south direction.
+For many coordinate reference systems this is also what is defined by the authority.
+There are however exceptions, especially when dealing with coordinate systems that
+don't align with the cardinal directions of a compass. For example it is not
+obvious which coordinate component aligns to which axis in a skewed coordinate
+system with a 45 degrees angle against the north direction. Similarly, a geocentric
+cartesain coordinate system usually has the z-component aligned with the rotational
+axis of the earth and hence the axis points towards north. Both cases are
+incompatible with the convention of always having the x-component be the east/west
+axis, the y-component the north/south axis and the z-component the up/down axis.
+
+In most cases coordinate reference systems with geodetic coordinates expect the
+input ordered as latitude/longitude  (typically with the EPSG dataset), however,
+internally PROJ expects an longitude/latitude ordering for all projections. This
+is generally hidden for users but in a few cases it is exposed at the surface
+level of PROJ, most prominently in the :program:`proj` utility which expects
+longitude/latitude ordering of input date (unless :program:`proj -r` is used).
+
+In case of doubt about the axis order of a specific CRS :program:`projinfo` is
+able to provide an answer. Simply look up the CRS and examine the axis specification
+of the Well-Known Text output:
+
+.. code-block:: console
+
+    projinfo EPSG:4326
+    PROJ.4 string:
+    +proj=longlat +datum=WGS84 +no_defs +type=crs
+
+    WKT2_2018 string:
+    GEOGCRS["WGS 84",
+        DATUM["World Geodetic System 1984",
+            ELLIPSOID["WGS 84",6378137,298.257223563,
+                LENGTHUNIT["metre",1]]],
+        PRIMEM["Greenwich",0,
+            ANGLEUNIT["degree",0.0174532925199433]],
+        CS[ellipsoidal,2],
+            AXIS["geodetic latitude (Lat)",north,
+                ORDER[1],
+                ANGLEUNIT["degree",0.0174532925199433]],
+            AXIS["geodetic longitude (Lon)",east,
+                ORDER[2],
+                ANGLEUNIT["degree",0.0174532925199433]],
+        USAGE[
+            SCOPE["unknown"],
+            AREA["World"],
+            BBOX[-90,-180,90,180]],
+        ID["EPSG",4326]]
diff --git a/docs/source/references.bib b/docs/source/references.bib
index 43b187b2..9829495f 100644
--- a/docs/source/references.bib
+++ b/docs/source/references.bib
@@ -160,6 +160,19 @@
   Url                      = {https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/}
 }
 
+@TechReport{ISO19111,
+  Title                    = {{Geographic information -- Referencing by coordinates}},
+  Author                   = {ISO},
+  Year                     = {2019},
+  Type                     = {Standard},
+  Key                      = {ISO 19111:2019},
+  Month                    = Jan,
+  Volume                   = {2019},
+  Address                  = {Geneva, CH},
+  Institution              = {International Organization for Standardization},
+  Url                      = {http://docs.opengeospatial.org/as/18-005r4/18-005r4.html}
+}
+
 @Article{Jenny2015,
   Title                    = {A compromise aspect-adaptive cylindrical projection for world maps},
   Author                   = {Bernhard Jenny and Bojan Šavrič and Tom Patterson},