[Backport 6.0] Doc: update quickstart with PROJ 6 API + add PROJ 4->6 migration guide (fixes #1403) (#1450)6.0

[Backport 6.0] Doc: update quickstart with PROJ 6 API + add PROJ 4->6 migration guide (fixes #1403)

4 files changed, 287 insertions, 41 deletions

diff --git a/docs/source/development/migration.rst b/docs/source/development/migration.rst
index 99143f20..df622ecb 100644
--- a/docs/source/development/migration.rst
+++ b/docs/source/development/migration.rst
@@ -1,6 +1,151 @@
 .. _API_migration:
 
 ================================================================================
+Version 4 to 6 API Migration
+================================================================================
+
+This is a transition guide for developers wanting to migrate their code to use
+PROJ version 6.
+
+Code example
+###############################################################################
+
+The difference between the old and new API is shown here with a few examples. Below
+we implement the same program with the two different API's. The program reads
+input longitude and latitude from the command line and convert them to
+projected coordinates with the Mercator projection.
+
+We start by writing the program for PROJ 4:
+
+.. code-block:: C
+
+    #include <proj_api.h>
+
+    main(int argc, char **argv) {
+        projPJ pj_merc, pj_latlong;
+        double x, y;
+
+        if (!(pj_longlat = pj_init_plus("+proj=longlat +ellps=clrk66")) )
+            return 1;
+        if (!(pj_merc = pj_init_plus("+proj=merc +datum=clrk66 +lat_ts=33")) )
+            return 1;
+
+        while (scanf("%lf %lf", &x, &y) == 2) {
+            x *= DEG_TO_RAD; /* longitude */
+            y *= DEG_TO_RAD; /* latitude */
+            p = pj_transform(pj_longlat, pj_merc, 1, 1, &x, &y, NULL );
+            printf("%.2f\t%.2f\n", x, y);
+        }
+
+        pj_free(pj_longlat);
+        pj_free(pj_merc);
+
+        return 0;
+    }
+
+The same program implemented using PROJ 6:
+
+.. code-block:: C
+
+    #include <proj.h>
+
+    main(int argc, char **argv) {
+        PJ *P;
+        PJ_COORD c;
+
+        /* NOTE: the use of PROJ strings to describe CRS is strongly discouraged */
+        /* in PROJ 6, as PROJ strings are a poor way of describing a CRS, and */
+        /* more precise its geodetic datum. */
+        /* Use of codes provided by authorities (such as "EPSG:4326", etc...) */
+        /* or WKT strings will bring the full power of the "transformation */
+        /* engine" used by PROJ to determine the best transformation(s) between */
+        /* two CRS. */
+        P = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
+                                   "+proj=longlat +ellps=clrs66",
+                                   "+proj=merc +ellps=clrk66 +lat_ts=33",
+                                   NULL);
+        if (P==0)
+            return 1;
+
+        {
+            /* For that particular use case, this is not needed. */
+            /* proj_normalize_for_visualization() ensures that the coordinate */
+            /* order expected and returned by proj_trans() will be longitude, */
+            /* latitude for geographic CRS, and easting, northing for projected */
+            /* CRS. If instead of using PROJ strings as above, "EPSG:XXXX" codes */
+            /* had been used, this might had been necessary. */
+            PJ* P_for_GIS = proj_normalize_for_visualization(C, P);
+            if( 0 == P_for_GIS )  {
+                proj_destroy(P);
+                return 1;
+            }
+            proj_destroy(P);
+            P = P_for_GIS;
+        }
+
+        while (scanf("%lf %lf", &c.lp.lam, &c.lp.phi) == 2) {
+            /* No need to convert to radian */
+            c = proj_trans(P, PJ_FWD, c);
+            printf("%.2f\t%.2f\n", c.xy.x, c.xy.y);
+        }
+
+        proj_destroy(P);
+    }
+
+
+Function mapping from old to new API
+###############################################################################
+
++---------------------------------------+-------------------------------------------------+
+| Old API functions                     | New API functions                               |
++=======================================+=================================================+
+| pj_fwd                                | :c:func:`proj_trans`                            |
++---------------------------------------+-------------------------------------------------+
+| pj_inv                                | :c:func:`proj_trans`                            |
++---------------------------------------+-------------------------------------------------+
+| pj_fwd3                               | :c:func:`proj_trans`                            |
++---------------------------------------+-------------------------------------------------+
+| pj_inv3                               | :c:func:`proj_trans`                            |
++---------------------------------------+-------------------------------------------------+
+| pj_transform                          | :c:func:`proj_create_crs_to_crs` +              |
+|                                       | (:c:func:`proj_normalize_for_visualization` +)  |
+|                                       | :c:func:`proj_trans`,                           |
+|                                       | :c:func:`proj_trans_array` or                   |
+|                                       | :c:func:`proj_trans_generic`                    |
++---------------------------------------+-------------------------------------------------+
+| pj_init                               | :c:func:`proj_create` /                         |
+|                                       | :c:func:`proj_create_crs_to_crs`                |
++---------------------------------------+-------------------------------------------------+
+| pj_init                               | :c:func:`proj_create` /                         |
+|                                       | :c:func:`proj_create_crs_to_crs`                |
++---------------------------------------+-------------------------------------------------+
+| pj_free                               | :c:func:`proj_destroy`                          |
++---------------------------------------+-------------------------------------------------+
+| pj_is_latlong                         | :c:func:`proj_get_type`                         |
++---------------------------------------+-------------------------------------------------+
+| pj_is_geocent                         | :c:func:`proj_get_type`                         |
++---------------------------------------+-------------------------------------------------+
+| pj_get_def                            | :c:func:`proj_pj_info`                          |
++---------------------------------------+-------------------------------------------------+
+| pj_latlong_from_proj                  | *No direct equivalent*, but can be accomplished |
+|                                       | by chaining :c:func:`proj_create`,              |
+|                                       | :c:func:`proj_crs_get_horizontal_datum` and     |
+|                                       | :c:func:`proj_create_geographic_crs_from_datum` |
++---------------------------------------+-------------------------------------------------+
+| pj_set_finder                         | :c:func:`proj_context_set_file_finder`          |
++---------------------------------------+-------------------------------------------------+
+| pj_set_searchpath                     | :c:func:`proj_context_set_search_paths`         |
++---------------------------------------+-------------------------------------------------+
+| pj_deallocate_grids                   | *No equivalent*                                 |
++---------------------------------------+-------------------------------------------------+
+| pj_strerrno                           | *No equivalent*                                 |
++---------------------------------------+-------------------------------------------------+
+| pj_get_errno_ref                      | :c:func:`proj_errno`                            |
++---------------------------------------+-------------------------------------------------+
+| pj_get_release                        | :c:func:`proj_info`                             |
++---------------------------------------+-------------------------------------------------+
+
+================================================================================
 Version 4 to 5 API Migration
 ================================================================================
 
@@ -66,7 +211,7 @@ Code example
 
 The difference between the old and new API is shown here with a few examples. Below
 we implement the same program with the two different API's. The program reads
-input latitude and longitude from the command line and convert them to
+input longitude and latitude from the command line and convert them to
 projected coordinates with the Mercator projection.
 
 We start by writing the program for PROJ v. 4:
@@ -76,21 +221,24 @@ We start by writing the program for PROJ v. 4:
     #include <proj_api.h>
 
     main(int argc, char **argv) {
-        projPJ pj_merc, pj_latlong;
+        projPJ pj_merc, pj_longlat;
         double x, y;
 
-        if (!(pj_merc = pj_init_plus("+proj=merc +ellps=clrk66 +lat_ts=33")) )
+        if (!(pj_longlat = pj_init_plus("+proj=longlat +ellps=clrk66")) )
             return 1;
-        if (!(pj_latlong = pj_init_plus("+proj=latlong +ellps=clrk66")) )
+        if (!(pj_merc = pj_init_plus("+proj=merc +ellps=clrk66 +lat_ts=33")) )
             return 1;
 
         while (scanf("%lf %lf", &x, &y) == 2) {
-            x *= DEG_TO_RAD;
-            y *= DEG_TO_RAD;
-            p = pj_transform(pj_latlong, pj_merc, 1, 1, &x, &y, NULL );
+            x *= DEG_TO_RAD; /* longitude */
+            y *= DEG_TO_RAD; /* latitude */
+            p = pj_transform(pj_longlat, pj_merc, 1, 1, &x, &y, NULL );
             printf("%.2f\t%.2f\n", x, y);
         }
 
+        pj_free(pj_longlat);
+        pj_free(pj_merc);
+
         return 0;
     }
 
@@ -115,6 +263,7 @@ The same program implemented using PROJ v. 5:
             printf("%.2f\t%.2f\n", c.xy.x, c.xy.y);
         }
 
+        proj_destroy(P);
     }
 
 Looking at the two different programs, there's a few immediate
@@ -155,27 +304,27 @@ Function mapping from old to new API
 +---------------------------------------+---------------------------------------+
 | Old API functions                     | New API functions                     |
 +=======================================+=======================================+
-| pj_fwd                                | proj_trans                            |
+| pj_fwd                                | :c:func:`proj_trans`                  |
 +---------------------------------------+---------------------------------------+
-| pj_inv                                | proj_trans                            |
+| pj_inv                                | :c:func:`proj_trans`                  |
 +---------------------------------------+---------------------------------------+
-| pj_fwd3                               | proj_trans                            |
+| pj_fwd3                               | :c:func:`proj_trans`                  |
 +---------------------------------------+---------------------------------------+
-| pj_inv3                               | proj_trans                            |
+| pj_inv3                               | :c:func:`proj_trans`                  |
 +---------------------------------------+---------------------------------------+
 | pj_transform                          | proj_trans_array or proj_trans_generic|
 +---------------------------------------+---------------------------------------+
-| pj_init                               | proj_create                           |
+| pj_init                               | :c:func:`proj_create`                 |
 +---------------------------------------+---------------------------------------+
-| pj_init_plus                          | proj_create                           |
+| pj_init_plus                          | :c:func:`proj_create`                 |
 +---------------------------------------+---------------------------------------+
-| pj_free                               | proj_destroy                          |
+| pj_free                               | :c:func:`proj_destroy`                |
 +---------------------------------------+---------------------------------------+
-| pj_is_latlong                         | proj_angular_output                   |
+| pj_is_latlong                         | :c:func:`proj_angular_output`         |
 +---------------------------------------+---------------------------------------+
-| pj_is_geocent                         | proj_angular_outout                   |
+| pj_is_geocent                         | :c:func:`proj_angular_output`         |
 +---------------------------------------+---------------------------------------+
-| pj_get_def                            | proj_pj_info                          |
+| pj_get_def                            | :c:func:`proj_pj_info`                |
 +---------------------------------------+---------------------------------------+
 | pj_latlong_from_proj                  | *No equivalent*                       |
 +---------------------------------------+---------------------------------------+
@@ -187,7 +336,7 @@ Function mapping from old to new API
 +---------------------------------------+---------------------------------------+
 | pj_strerrno                           | *No equivalent*                       |
 +---------------------------------------+---------------------------------------+
-| pj_get_errno_ref                      | proj_errno                            |
+| pj_get_errno_ref                      | :c:func:`proj_errno`                  |
 +---------------------------------------+---------------------------------------+
-| pj_get_release                        | proj_info                             |
+| pj_get_release                        | :c:func:`proj_info`                   |
 +---------------------------------------+---------------------------------------+
diff --git a/docs/source/development/quickstart.rst b/docs/source/development/quickstart.rst
index b6c5d17c..d53d98fd 100644
--- a/docs/source/development/quickstart.rst
+++ b/docs/source/development/quickstart.rst
@@ -25,7 +25,7 @@ See the :doc:`reference for more info on data types <reference/datatypes>`.
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 43-45
+  :lines: 43-46
   :dedent: 4
 
 For use in multi-threaded programs the :c:type:`PJ_CONTEXT` threading-context is used.
@@ -35,13 +35,43 @@ this in detail.
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 49
+  :lines: 50
   :dedent: 4
 
 Next we create the :c:type:`PJ` transformation object ``P`` with the function
-:c:func:`proj_create`. :c:func:`proj_create` takes the threading context ``C``
-created above, and a proj-string that defines the desired transformation.
-Here we transform from geodetic coordinate to UTM zone 32N.
+:c:func:`proj_create_crs_to_crs`. :c:func:`proj_create_crs_to_crs` takes the threading context ``C``
+created above, a string that describes the source coordinate reference system (CRS),
+a string that describes the target CRS and an optional description of the area of
+use.
+The strings for the source or target CRS may be PROJ strings (``+proj=longlat +datum=WGS84``),
+CRS identified by their code (``EPSG:4326`` or ``urn:ogc:def:crs:EPSG::4326``) or
+by a well-known text (WKT) 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]]
+
+).
+The use of PROJ strings to describe a CRS is considered as legacy (one of the
+main weakness of PROJ strings is their inability to describe a geodetic datum,
+other than the few ones hardcoded in the ``+datum`` parameter).
+Here we transform from geographic coordinates to UTM zone 32N.
 It is recommended to create one threading-context per thread used by the program.
 This ensures that all :c:type:`PJ` objects created in the same context will be
 sharing resources such as error-numbers and loaded grids.
@@ -51,43 +81,72 @@ details.
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 51-53
+  :lines: 52-60
   :dedent: 4
 
-PROJ uses it's own data structures for handling coordinates. Here we use a
+:c:func:`proj_create_crs_to_crs` creates a transformation object, which accepts
+coordinates expressed in the units and axis order of the definition of the
+source CRS, and return transformed coordinates in the units and axis order of
+the definition of the target CRS.
+For almost most geographic CRS, the units will be in most cases degrees (in
+rare cases, such as EPSG:4807 / NTF (Paris), this can be grads). For geographic
+CRS defined by the EPSG authority, the order of coordinates is latitude first,
+longitude second. When using a PROJ string, on contrary the order will be
+longitude first, latitude second.
+For projected CRS, the units may vary (metre, us-foot, etc..). For projected
+CRS defined by the EPSG authority, and with EAST / NORTH directions, the order
+might be easting first, northing second, or the reverse.  When using a PROJ string,
+the order will be easting first, northing second, except if the ``+axis``
+parameter modifies it.
+
+If for the needs of your software, you want
+a uniform axis order (and thus do not care about axis order mandated by the
+authority defining the CRS), the :c:func:`proj_normalize_for_visualization`
+function can be used to modify the PJ* object returned by
+:c:func:`proj_create_crs_to_crs` so that it accepts as input and returns as
+output coordinates using the traditional GIS order, that is longitude, latitude
+(followed by elevation, time) for geographic CRS and easting, northing for most
+projected CRS.
+
+.. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
+  :language: c
+  :lines: 65-71
+  :dedent: 4
+
+PROJ uses its own data structures for handling coordinates. Here we use a
 :c:type:`PJ_COORD` which is easily assigned with the function :c:func:`proj_coord`.
-Note that the input values are converted to radians with :c:func:`proj_torad`.
-This is necessary since PROJ is using radians internally. See :doc:`transformations`
-for further details.
+When using +proj=longlat, the order of coordinates is longitude, latitude,
+and values are expressed in degrees. If you used instead a EPSG geographic CRS,
+like EPSG:4326 (WGS84), it would be latitude, longitude.
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 57
+  :lines: 76
   :dedent: 4
 
 The coordinate defined above is transformed with :c:func:`proj_trans`. For this
 a :c:type:`PJ` object, a transformation direction (either forward or inverse)
 and the coordinate is needed. The transformed coordinate is returned in ``b``.
-Here the forward (:c:type:`PJ_FWD`) transformation from geodetic to UTM is made.
+Here the forward (:c:type:`PJ_FWD`) transformation from geographic to UTM is made.
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 60-61
+  :lines: 79-80
   :dedent: 4
 
-The inverse transformation (UTM to geodetic) is done similar to above,
+The inverse transformation (UTM to geographic) is done similar to above,
 this time using :c:type:`PJ_INV` as the direction.
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 62-63
+  :lines: 81-82
   :dedent: 4
 
 Before ending the program the allocated memory needs to be released again:
 
 .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c
   :language: c
-  :lines: 66-67
+  :lines: 85-86
   :dedent: 4
 
 
diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index 00653ad9..319face3 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -144,6 +144,25 @@ paragraph for more details.
     :type `area`: PJ_AREA
     :returns: :c:type:`PJ*`
 
+.. c:function:: PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ* obj)
+
+    .. versionadded:: 6.1.0
+
+    Returns a PJ* object whose axis order is the one expected for
+    visualization purposes.
+
+    The input object must be a coordinate operation, that has been created with
+    proj_create_crs_to_crs().
+    If the axis order of its source or target CRS is northing,easting, then an
+    axis swap operation will be inserted.
+
+    The returned :c:type:`PJ`-pointer should be deallocated with :c:func:`proj_destroy`.
+
+    :param PJ_CONTEXT* ctx: Threading context.
+    :param `obj`: Object of type CoordinateOperation
+    :returns: :c:type:`PJ*`
+
+
 .. c:function:: PJ* proj_destroy(PJ *P)
 
     Deallocate a :c:type:`PJ` transformation object.
diff --git a/examples/pj_obs_api_mini_demo.c b/examples/pj_obs_api_mini_demo.c
index 94520490..3df94e2d 100644
--- a/examples/pj_obs_api_mini_demo.c
+++ b/examples/pj_obs_api_mini_demo.c
@@ -42,23 +42,42 @@
 int main (void) {
     PJ_CONTEXT *C;
     PJ *P;
+    PJ* P_for_GIS;
     PJ_COORD a, b;
 
     /* or you may set C=PJ_DEFAULT_CTX if you are sure you will     */
     /* use PJ objects from only one thread                          */
     C = proj_context_create();
 
-    P = proj_create (C, "+proj=utm +zone=32 +ellps=GRS80");
-    if (0==P)
-        return puts ("Oops"), 0;
+    P = proj_create_crs_to_crs (C,
+                                "EPSG:4326",
+                                "+proj=utm +zone=32 +datum=WGS84", /* or EPSG:32632 */
+                                NULL);
+
+    if (0==P) {
+        fprintf(stderr, "Oops\n");
+        return 1;
+    }
+
+    /* This will ensure that the order of coordinates for the input CRS */
+    /* will be longitude, latitude, whereas EPSG:4326 mandates latitude, */
+    /* longitude */
+    P_for_GIS = proj_normalize_for_visualization(C, P);
+    if( 0 == P_for_GIS )  {
+        fprintf(stderr, "Oops\n");
+        return 1;
+    }
+    proj_destroy(P);
+    P = P_for_GIS;
 
     /* a coordinate union representing Copenhagen: 55d N, 12d E    */
-    /* note: PROJ.4 works in radians, hence the proj_torad() calls */
-    a = proj_coord (proj_torad(12), proj_torad(55), 0, 0);
+    /* Given that we have used proj_normalize_for_visualization(), the order of
+    /* coordinates is longitude, latitude, and values are expressed in degrees. */
+    a = proj_coord (12, 55, 0, 0);
 
     /* transform to UTM zone 32, then back to geographical */
     b = proj_trans (P, PJ_FWD, a);
-    printf ("easting: %g, northing: %g\n", b.enu.e, b.enu.n);
+    printf ("easting: %.3f, northing: %.3f\n", b.enu.e, b.enu.n);
     b = proj_trans (P, PJ_INV, b);
     printf ("longitude: %g, latitude: %g\n", b.lp.lam, b.lp.phi);