diff options
| -rw-r--r-- | docs/source/development/quickstart.rst | 91 | ||||
| -rw-r--r-- | docs/source/development/reference/functions.rst | 19 | ||||
| -rw-r--r-- | examples/pj_obs_api_mini_demo.c | 31 |
3 files changed, 119 insertions, 22 deletions
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); |