Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add C++ API to transform coordinate and a C++ quickstart #3705

Merged
merged 6 commits into from
Apr 19, 2023
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Doc: add a C++ quickstart (fixes #3593)
rouault committed Apr 9, 2023
commit 66dbb7f58b0dc4557114e03ec688b990f1afca48
1 change: 1 addition & 0 deletions docs/source/development/index.rst
Original file line number Diff line number Diff line change
@@ -13,6 +13,7 @@ PROJ project or using the library in their own software.

dev_practices
quickstart
quickstart_cpp
errorhandling
reference/index
cmake
4 changes: 2 additions & 2 deletions docs/source/development/quickstart.rst
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
.. _dev_quickstart:

================================================================================
Quick start
Quick start for C API usage
================================================================================

This is a short introduction to the PROJ API. In the following section
This is a short introduction to the PROJ C API. In the following section
we create two simple programs that illustrate how to transform points
between two different coordinate systems, and how to convert between
projected and geodetic (geographic) coordinates for a single
160 changes: 160 additions & 0 deletions docs/source/development/quickstart_cpp.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,160 @@
.. _dev_quickstart_cpp:

================================================================================
Quick start for C++ API usage
================================================================================

This is a short introduction to the PROJ :ref:`C++ API <cpp>`. In the following section
we create a simple program that illustrate how to transform points
between two different coordinate systems.

Note: the essential osgeo::proj::operation::CoordinateOperation::coordinateTransformer()
method was only added in PROJ 9.3. For earlier versions, coordinate operations
must be exported as a PROJ string with the osgeo::proj::operation::CoordinateOperation::exportToPROJ() method, and
passed to :c:func:`proj_create` to instance a PJ* object to use with :c:func:`proj_trans`
(cf the :ref:`dev_quickstart`). The use of the C API :c:func:`proj_create_crs_to_crs`
might still be easier to let PROJ automatically select the "best" coordinate
operation when several ones are possible, as the corresponding automation is not
currently available in the C++ API.

Before the PROJ API can be used it is necessary to include the various header
files:

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 9-18

For convenience, we also declare using a few namespaces:

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 20-23

We start by creating a database context (:cpp:class:`osgeo::proj::io::DatabaseContext`)
with the default settings to find the PROJ database.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 26
:dedent: 4

We then instantiate a generic authority factory (:cpp:class:`osgeo::proj::io::AuthorityFactory`),
that is not tied to a particular authority, to be able to get transformations
registered by different authorities. This can only be used for a
:cpp:class:`osgeo::proj::operation::CoordinateOperationContext`, and not to
instantiate objects of the database which are all tied to a non-generic authority.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 31

We create a coordinate operation context, that can be customized to ammend
the way coordinate operations are computed. Here we ask for default settings,
as we have a coordinate operation that just involves a "simple" map projection
in the same datum.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 36-37
:dedent: 4

We instantiate a authority factory for EPSG related objects.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 40
:dedent: 4

We instantiate the source CRS from its code: 4326, for WGS 84 latitude/longitude.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 43
:dedent: 4

We instantiate the source CRS from its PROJ.4 string (it would be possible
to instantiate it from its 32631 code, similarly to above), and cast the
generic :cpp:class:`osgeo::proj::util::BaseObject` to the
:cpp:class:`osgeo::proj::crs::CRS` class required later.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 49-51
:dedent: 4

.. warning::

The use of PROJ strings to describe a CRS is not recommended. One of the
main weaknesses of PROJ strings is their inability to describe a geodetic
datum, other than the few ones hardcoded in the ``+datum`` parameter.

We ask for the list of operations available to transform from the source
to the target CRS with the :cpp:func:`osgeo::proj::operation::CoordinateOperationFactory::createOperations`
method.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 55-56
:dedent: 4

We check that we got a non-empty list of operations. The list is sorted from
the most relevant to the less relevant one. Cf :ref:`operations_computation_filtering`
for more details on the sorting of those operations. For a transformation
between a projected CRS and its base CRS, like we do here, there will be only
one operation.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 65
:dedent: 4

We create an execution context (must only be used by one thread at a time)
with the :c:func:`proj_context_create` function.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 68
:dedent: 4

We create a coordinate transformer (:cpp:class:`osgeo::proj::operation::CoordinateTransformer`)
from the first operation of the list:

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 71
:dedent: 4

We can now transform a point with the :cpp:func:`osgeo::proj::operation::CoordinateTransformer::transform`
method. Note that the the expected input values should be passed in the order
and the unit of the successive axis of the input CRS. Similarly the values
returned in the v[] array of the output PJ_COORD are in the order and the unit
of the successive axis of the output CRS. For coordinate operations involving a
time-dependent operation, coord.v[3] is the decimal year of the coordinate epoch
of the input (or HUGE_VAL to indicate none).

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 74-80
:dedent: 4

and output the result:

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 83-85
:dedent: 4

We need to clean up the PJ_CONTEXT handle before exiting with the
:c:func:`proj_context_destroy` function.

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:lines: 88
:dedent: 4

A complete compilable version of the example code can be seen below:

.. literalinclude:: ../../../examples/EPSG_4326_to_32631.cpp
:language: c++
:linenos:
:lines: 9-
2 changes: 2 additions & 0 deletions docs/source/operations/operations_computation.rst
Original file line number Diff line number Diff line change
@@ -109,6 +109,8 @@ does not cover the area of use of the 2 CRSs, a
'Ballpark geographic offset from NAD27 to NAD83' operation is synthetized by PROJ
(see :term:`Ballpark transformation`)

.. _operations_computation_filtering:

Filtering and sorting of coordinate operations
----------------------------------------------

3 changes: 3 additions & 0 deletions examples/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -3,3 +3,6 @@ target_link_libraries(pj_obs_api_mini_demo PRIVATE ${PROJ_LIBRARIES})

add_executable(crs_to_geodetic crs_to_geodetic.c)
target_link_libraries(crs_to_geodetic PRIVATE ${PROJ_LIBRARIES})

add_executable(EPSG_4326_to_32631 EPSG_4326_to_32631.cpp)
target_link_libraries(EPSG_4326_to_32631 PRIVATE ${PROJ_LIBRARIES})
91 changes: 91 additions & 0 deletions examples/EPSG_4326_to_32631.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
/*******************************************************************************
Example code that illustrates how to convert between two CRS

Note: This file is in-lined in the documentation. Any changes must be
reflected in docs/source/development/quickstart.rst

*******************************************************************************/

#include <cassert>
#include <cmath> // for HUGE_VAL
#include <iomanip> // for std::setprecision()

#include <iostream>

#include "proj/coordinateoperation.hpp"
#include "proj/crs.hpp"
#include "proj/io.hpp"
#include "proj/util.hpp" // for nn_dynamic_pointer_cast

using namespace NS_PROJ::crs;
using namespace NS_PROJ::io;
using namespace NS_PROJ::operation;
using namespace NS_PROJ::util;

int main(void) {
auto dbContext = DatabaseContext::create();

// Instantiate a generic authority factory, that is not tied to a particular
// authority, to be able to get transformations registered by different
// authorities. This can only be used for CoordinateOperationContext.
auto authFactory = AuthorityFactory::create(dbContext, std::string());

// Create a coordinate operation context, that can be customized to ammend
// the way coordinate operations are computed. Here we ask for default
// settings.
auto coord_op_ctxt =
CoordinateOperationContext::create(authFactory, nullptr, 0.0);

// Instantiate a authority factory for EPSG related objects.
auto authFactoryEPSG = AuthorityFactory::create(dbContext, "EPSG");

// Instantiate source CRS from EPSG code
auto sourceCRS = authFactoryEPSG->createCoordinateReferenceSystem("4326");

// Instantiate target CRS from PROJ.4 string (commented out, the equivalent
// from the EPSG code)
// auto targetCRS =
// authFactoryEPSG->createCoordinateReferenceSystem("32631");
auto targetCRS =
NN_CHECK_THROW(nn_dynamic_pointer_cast<CRS>(createFromUserInput(
"+proj=utm +zone=31 +datum=WGS84 +type=crs", dbContext)));

// List operations available to transform from EPSG:4326
// (WGS 84 latitude/longitude) to EPSG:32631 (WGS 84 / UTM zone 31N).
auto list = CoordinateOperationFactory::create()->createOperations(
sourceCRS, targetCRS, coord_op_ctxt);

// Check that we got a non-empty list of operations
// The list is sorted from the most relevant to the less relevant one.
// Cf
// https://proj.org/operations/operations_computation.html#filtering-and-sorting-of-coordinate-operations
// for more details on the sorting of those operations.
// For a transformation between a projected CRS and its base CRS, like
// we do here, there will be only one operation.
assert(!list.empty());

// Create an execution context (must only be used by one thread at a time)
PJ_CONTEXT *ctx = proj_context_create();

// Create a coordinate transformer from the first operation of the list
auto transformer = list[0]->coordinateTransformer(ctx);

// Perform the coordinate transformation.
PJ_COORD c = {{
49.0, // latitude in degree
2.0, // longitude in degree
0.0, // z ordinate. unused
HUGE_VAL // time ordinate. unused
}};
c = transformer->transform(c);

// Display result
std::cout << std::fixed << std::setprecision(3);
std::cout << "Easting: " << c.v[0] << std::endl; // should be 426857.988
std::cout << "Northing: " << c.v[1] << std::endl; // should be 5427937.523

// Destroy execution context
proj_context_destroy(ctx);

return 0;
}