Merge efd0752 into e1da8e5

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-55
   :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: 60-66
+  :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
 
 

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.

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);
+
+    /* 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;
+
+    if (0==P) {
+        fprintf(stderr, "Oops\n");
+        return 1;
+    }
 
     /* 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);