aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--CMakeLists.txt23
-rw-r--r--cmake/ProjUtilities.cmake289
-rw-r--r--cmake/ProjVersion.cmake6
-rw-r--r--docs/source/community/rfc/index.rst1
-rw-r--r--docs/source/community/rfc/rfc-2.rst933
-rw-r--r--docs/source/development/migration.rst4
-rw-r--r--docs/source/operations/conversions/unitconvert.rst12
-rw-r--r--docs/source/operations/projections/eqc.rst2
-rw-r--r--docs/source/operations/projections/tmerc.rst2
-rw-r--r--docs/source/operations/transformations/helmert.rst6
-rw-r--r--docs/source/operations/transformations/index.rst1
-rw-r--r--docs/source/operations/transformations/molobadekas.rst147
-rw-r--r--docs/source/resource_files.rst7
-rw-r--r--src/PJ_helmert.c240
-rw-r--r--src/pj_list.h1
-rw-r--r--test/gie/more_builtins.gie24
16 files changed, 1311 insertions, 387 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 69e93deb..e6265dc0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,7 +9,7 @@
#################################################################################
# General settings
#################################################################################
-cmake_minimum_required(VERSION 2.6.0 FATAL_ERROR)
+cmake_minimum_required(VERSION 2.8.3 FATAL_ERROR)
# For historic reasons, the CMake PROJECT-NAME is PROJ4
project(PROJ4 LANGUAGES C CXX)
@@ -82,11 +82,6 @@ set(CMAKE_THREAD_PREFER_PTHREAD TRUE)
find_package (Threads)
include(CheckIncludeFiles)
-include(CheckSymbolExists)
-CHECK_SYMBOL_EXISTS(PTHREAD_MUTEX_RECURSIVE pthread.h HAVE_PTHREAD_MUTEX_RECURSIVE_DEFN)
-if (HAVE_PTHREAD_MUTEX_RECURSIVE_DEFN)
- add_definitions(-DHAVE_PTHREAD_MUTEX_RECURSIVE=1)
-endif()
include (CheckCSourceCompiles)
if (MSVC)
@@ -112,6 +107,22 @@ else ()
add_definitions (-DHAVE_C99_MATH=0)
endif ()
+if (Threads_FOUND AND CMAKE_USE_PTHREADS_INIT)
+ set (CMAKE_REQUIRED_LIBRARIES "${CMAKE_REQUIRED_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT}")
+ check_c_source_compiles("
+#include <pthread.h>
+
+int main(int argc, char* argv[]) {
+ (void)PTHREAD_MUTEX_RECURSIVE;
+ (void)argv;
+ return argc;
+}
+ " HAVE_PTHREAD_MUTEX_RECURSIVE_DEFN)
+ if (HAVE_PTHREAD_MUTEX_RECURSIVE_DEFN)
+ add_definitions(-DHAVE_PTHREAD_MUTEX_RECURSIVE=1)
+ endif()
+endif ()
+
boost_report_value(PROJ_PLATFORM_NAME)
boost_report_value(PROJ_COMPILER_NAME)
diff --git a/cmake/ProjUtilities.cmake b/cmake/ProjUtilities.cmake
index 70d3d518..6969a664 100644
--- a/cmake/ProjUtilities.cmake
+++ b/cmake/ProjUtilities.cmake
@@ -13,228 +13,10 @@
################################################################################
# Macros in this module:
#
-# list_contains: Determine whether a string value is in a list.
-#
-# car: Return the first element in a list
-#
-# cdr: Return all but the first element in a list
-#
-# parse_arguments: Parse keyword arguments for use in other macros.
-#
-# proj_report_directory_property
-#
# proj_target_output_name:
#
################################################################################
-# This utility macro determines whether a particular string value
-# occurs within a list of strings:
-#
-# list_contains(result string_to_find arg1 arg2 arg3 ... argn)
-#
-# This macro sets the variable named by result equal to TRUE if
-# string_to_find is found anywhere in the following arguments.
-macro(list_contains var value)
- set(${var})
- foreach (value2 ${ARGN})
- if (${value} STREQUAL ${value2})
- set(${var} TRUE)
- endif (${value} STREQUAL ${value2})
- endforeach (value2)
-endmacro(list_contains)
-
-# This utility macro extracts the first argument from the list of
-# arguments given, and places it into the variable named var.
-#
-# car(var arg1 arg2 ...)
-macro(car var)
- set(${var} ${ARGV1})
-endmacro(car)
-
-# This utility macro extracts all of the arguments given except the
-# first, and places them into the variable named var.
-#
-# car(var arg1 arg2 ...)
-macro(cdr var junk)
- set(${var} ${ARGN})
-endmacro(cdr)
-
-# The parse_arguments macro will take the arguments of another macro and
-# define several variables. The first argument to parse_arguments is a
-# prefix to put on all variables it creates. The second argument is a
-# list of names, and the third argument is a list of options. Both of
-# these lists should be quoted. The rest of parse_arguments are
-# arguments from another macro to be parsed.
-#
-# parse_arguments(prefix arg_names options arg1 arg2...)
-#
-# For each item in options, parse_arguments will create a variable with
-# that name, prefixed with prefix_. So, for example, if prefix is
-# MY_MACRO and options is OPTION1;OPTION2, then parse_arguments will
-# create the variables MY_MACRO_OPTION1 and MY_MACRO_OPTION2. These
-# variables will be set to true if the option exists in the command line
-# or false otherwise.
-#
-# For each item in arg_names, parse_arguments will create a variable
-# with that name, prefixed with prefix_. Each variable will be filled
-# with the arguments that occur after the given arg_name is encountered
-# up to the next arg_name or the end of the arguments. All options are
-# removed from these lists. parse_arguments also creates a
-# prefix_DEFAULT_ARGS variable containing the list of all arguments up
-# to the first arg_name encountered.
-macro(parse_arguments prefix arg_names option_names)
- set(DEFAULT_ARGS)
- foreach(arg_name ${arg_names})
- set(${prefix}_${arg_name})
- endforeach(arg_name)
- foreach(option ${option_names})
- set(${prefix}_${option} FALSE)
- endforeach(option)
-
- set(current_arg_name DEFAULT_ARGS)
- set(current_arg_list)
- foreach(arg ${ARGN})
- list_contains(is_arg_name ${arg} ${arg_names})
- if (is_arg_name)
- set(${prefix}_${current_arg_name} ${current_arg_list})
- set(current_arg_name ${arg})
- set(current_arg_list)
- else (is_arg_name)
- list_contains(is_option ${arg} ${option_names})
- if (is_option)
- set(${prefix}_${arg} TRUE)
- else (is_option)
- set(current_arg_list ${current_arg_list} ${arg})
- endif (is_option)
- endif (is_arg_name)
- endforeach(arg)
- set(${prefix}_${current_arg_name} ${current_arg_list})
-endmacro(parse_arguments)
-
-# Perform a reverse topological sort on the given LIST.
-#
-# topological_sort(my_list "MY_" "_EDGES")
-#
-# LIST is the name of a variable containing a list of elements to be
-# sorted in reverse topological order. Each element in the list has a
-# set of outgoing edges (for example, those other list elements that
-# it depends on). In the resulting reverse topological ordering
-# (written back into the variable named LIST), an element will come
-# later in the list than any of the elements that can be reached by
-# following its outgoing edges and the outgoing edges of any vertices
-# they target, recursively. Thus, if the edges represent dependencies
-# on build targets, for example, the reverse topological ordering is
-# the order in which one would build those targets.
-#
-# For each element E in this list, the edges for E are contained in
-# the variable named ${PREFIX}${E}${SUFFIX}, where E is the
-# upper-cased version of the element in the list. If no such variable
-# exists, then it is assumed that there are no edges. For example, if
-# my_list contains a, b, and c, one could provide a dependency graph
-# using the following variables:
-#
-# MY_A_EDGES b
-# MY_B_EDGES
-# MY_C_EDGES a b
-#
-# With the involcation of topological_sort shown above and these
-# variables, the resulting reverse topological ordering will be b, a,
-# c.
-function(topological_sort LIST PREFIX SUFFIX)
- # Clear the stack and output variable
- set(VERTICES "${${LIST}}")
- set(STACK)
- set(${LIST})
-
- # Loop over all of the vertices, starting the topological sort from
- # each one.
- foreach(VERTEX ${VERTICES})
- string(TOUPPER ${VERTEX} UPPER_VERTEX)
-
- # If we haven't already processed this vertex, start a depth-first
- # search from where.
- if (NOT FOUND_${UPPER_VERTEX})
- # Push this vertex onto the stack with all of its outgoing edges
- string(REPLACE ";" " " NEW_ELEMENT
- "${VERTEX};${${PREFIX}${UPPER_VERTEX}${SUFFIX}}")
- list(APPEND STACK ${NEW_ELEMENT})
-
- # We've now seen this vertex
- set(FOUND_${UPPER_VERTEX} TRUE)
-
- # While the depth-first search stack is not empty
- list(LENGTH STACK STACK_LENGTH)
- while(STACK_LENGTH GREATER 0)
- # Remove the vertex and its remaining out-edges from the top
- # of the stack
- list(GET STACK -1 OUT_EDGES)
- list(REMOVE_AT STACK -1)
-
- # Get the source vertex and the list of out-edges
- separate_arguments(OUT_EDGES)
- list(GET OUT_EDGES 0 SOURCE)
- list(REMOVE_AT OUT_EDGES 0)
-
- # While there are still out-edges remaining
- list(LENGTH OUT_EDGES OUT_DEGREE)
- while (OUT_DEGREE GREATER 0)
- # Pull off the first outgoing edge
- list(GET OUT_EDGES 0 TARGET)
- list(REMOVE_AT OUT_EDGES 0)
-
- string(TOUPPER ${TARGET} UPPER_TARGET)
- if (NOT FOUND_${UPPER_TARGET})
- # We have not seen the target before, so we will traverse
- # its outgoing edges before coming back to our
- # source. This is the key to the depth-first traversal.
-
- # We've now seen this vertex
- set(FOUND_${UPPER_TARGET} TRUE)
-
- # Push the remaining edges for the current vertex onto the
- # stack
- string(REPLACE ";" " " NEW_ELEMENT
- "${SOURCE};${OUT_EDGES}")
- list(APPEND STACK ${NEW_ELEMENT})
-
- # Setup the new source and outgoing edges
- set(SOURCE ${TARGET})
- string(TOUPPER ${SOURCE} UPPER_SOURCE)
- set(OUT_EDGES
- ${${PREFIX}${UPPER_SOURCE}${SUFFIX}})
- endif(NOT FOUND_${UPPER_TARGET})
-
- list(LENGTH OUT_EDGES OUT_DEGREE)
- endwhile (OUT_DEGREE GREATER 0)
-
- # We have finished all of the outgoing edges for
- # SOURCE; add it to the resulting list.
- list(APPEND ${LIST} ${SOURCE})
-
- # Check the length of the stack
- list(LENGTH STACK STACK_LENGTH)
- endwhile(STACK_LENGTH GREATER 0)
- endif (NOT FOUND_${UPPER_VERTEX})
- endforeach(VERTEX)
-
- set(${LIST} ${${LIST}} PARENT_SCOPE)
-endfunction(topological_sort)
-
-# Small little hack that tweaks a component name (as used for CPack)
-# to make sure to avoid certain names that cause problems. Sets the
-# variable named varname to the "sanitized" name.
-#
-# FIXME: This is a complete hack. We probably need to fix the CPack
-# generators (NSIS in particular) to get rid of the need for this.
-macro(fix_cpack_component_name varname name)
- if (${name} STREQUAL "foreach")
- set(${varname} "boost_foreach")
- else()
- set(${varname} ${name})
- endif()
-endmacro()
-
-
#
# A big shout out to the cmake gurus @ compiz
#
@@ -305,77 +87,6 @@ function(boost_report_value NAME)
colormsg("${NAME}${varpadding} = ${${NAME}}")
endfunction()
-function(trace NAME)
- if(BOOST_CMAKE_TRACE)
- string(LENGTH "${NAME}" varlen)
- math(EXPR padding_len 40-${varlen})
- string(SUBSTRING "........................................"
- 0 ${padding_len} varpadding)
- message("${NAME} ${varpadding} ${${NAME}}")
- endif()
-endfunction()
-
-#
-# pretty-prints the value of a variable so that the
-# equals signs align
-#
-function(boost_report_pretty PRETTYNAME VARNAME)
- string(LENGTH "${PRETTYNAME}" varlen)
- math(EXPR padding_len 30-${varlen})
- string(SUBSTRING " "
- 0 ${padding_len} varpadding)
- message(STATUS "${PRETTYNAME}${varpadding} = ${${VARNAME}}")
-endfunction()
-
-#
-# assert that ARG is actually a library target
-#
-
-macro(dependency_check ARG)
- trace(ARG)
- if (NOT "${ARG}" STREQUAL "")
- get_target_property(deptype ${ARG} TYPE)
- if(NOT deptype MATCHES ".*_LIBRARY$")
- set(DEPENDENCY_OKAY FALSE)
- list(APPEND DEPENDENCY_FAILURES ${ARG})
- endif()
- endif()
-endmacro()
-
-
-
-#
-# Pretty-print of given property of current directory.
-#
-macro(proj_report_directory_property PROPNAME)
- get_directory_property(${PROPNAME} ${PROPNAME})
- boost_report_value(${PROPNAME})
-endmacro()
-
-#
-# Scans the current directory and returns a list of subdirectories.
-# Author: Robert Fleming
-# Source: https://www.cmake.org/pipermail/cmake/2008-February/020114.html
-#
-# Third parameter is 1 if you want relative paths returned.
-# Usage: list_subdirectories(the_list_is_returned_here /path/to/project TRUE)
-#
-
-macro(list_subdirectories retval curdir return_relative)
- file(GLOB sub-dir RELATIVE ${curdir} *)
- set(list_of_dirs "")
- foreach(dir ${sub-dir})
- if(IS_DIRECTORY ${curdir}/${dir})
- if (${return_relative})
- set(list_of_dirs ${list_of_dirs} ${dir})
- else()
- set(list_of_dirs ${list_of_dirs} ${curdir}/${dir})
- endif()
- endif()
- endforeach()
- set(${retval} ${list_of_dirs})
-endmacro()
-
#
# Generates output name for given target depending on platform and version.
# For instance, on Windows, libraries get ABI version suffix proj_X_Y.{dll|lib}.
diff --git a/cmake/ProjVersion.cmake b/cmake/ProjVersion.cmake
index 1121cd36..f761dc47 100644
--- a/cmake/ProjVersion.cmake
+++ b/cmake/ProjVersion.cmake
@@ -17,8 +17,12 @@
# MAJOR.MINOR version is used to set SOVERSION
#
+include(CMakeParseArguments)
+
macro(proj_version)
- parse_arguments(THIS_VERSION "MAJOR;MINOR;PATCH;"
+ cmake_parse_arguments(THIS_VERSION
+ ""
+ "MAJOR;MINOR;PATCH"
""
${ARGN})
diff --git a/docs/source/community/rfc/index.rst b/docs/source/community/rfc/index.rst
index dc0332df..f0f0294c 100644
--- a/docs/source/community/rfc/index.rst
+++ b/docs/source/community/rfc/index.rst
@@ -12,3 +12,4 @@ the project.
:maxdepth: 1
rfc-1
+ rfc-2
diff --git a/docs/source/community/rfc/rfc-2.rst b/docs/source/community/rfc/rfc-2.rst
new file mode 100644
index 00000000..20526f70
--- /dev/null
+++ b/docs/source/community/rfc/rfc-2.rst
@@ -0,0 +1,933 @@
+.. _rfc2:
+
+====================================================================
+PROJ RFC 2: Initial integration of "GDAL SRS barn" work
+====================================================================
+
+:Author: Even Rouault
+:Contact: even.rouault at spatialys.com
+:Status: Adopted (not yet merged into master)
+:Initial version: 2018-10-09
+:Last Updated: 2018-10-31
+
+Summary
+-------
+
+This RFC is the result of a first phase of the `GDAL Coordinate System Barn Raising`_
+efforts. In its current state, this work mostly consists of:
+
+ - a C++ implementation of the ISO-19111:2018 / OGC Topic 2 "Referencing by
+ coordinates" classes to represent Datums, Coordinate systems, CRSs
+ (Coordinate Reference Systems) and Coordinate Operations.
+ - methods to convert between this C++ modeling and WKT1, WKT2 and PROJ string representations of those objects
+ - management and query of a SQLite3 database of CRS and Coordinate Operation definition
+ - a C API binding part of those capabilities
+
+.. _`GDAL Coordinate System Barn Raising`: https://gdalbarn.com/
+
+
+Related standards
+-----------------
+
+Consult `Applicable standards`_
+
+.. _`Applicable standards`: http://even.rouault.free.fr/proj_cpp_api/html/general_doc.html#standards
+
+(They will be linked from the PROJ documentation)
+
+
+Details
+-------
+
+Structure in packages / namespaces
+**********************************
+
+The C++ implementation of the (upcoming) ISO-19111:2018 / OGC Topic 2 "Referencing by
+coordinates" classes follows this abstract modeling as much as possible, using
+package names as C++ namespaces, abstract classes and method names. A new
+BoundCRS class has been added to cover the modeling of the WKT2 BoundCRS
+construct, that is a generalization of the WKT1 TOWGS84 concept. It is
+strongly recommended to have the ISO-19111 standard open to have an introduction
+for the concepts when looking at the code. A few classes have also been
+inspired by the GeoAPI
+
+The classes are organized into several namespaces:
+
+ - osgeo::proj::util
+ A set of base types from ISO 19103, GeoAPI and other PROJ "technical"
+ specific classes
+
+ Template optional<T>, classes BaseObject, IComparable, BoxedValue,
+ ArrayOfBaseObject, PropertyMap, LocalName, NameSpace, GenericName,
+ NameFactory, CodeList, Exception, InvalidValueTypeException,
+ UnsupportedOperationException
+
+ - osgeo::proj::metadata:
+ Common classes from ISO 19115 (Metadata) standard
+
+ Classes Citation, GeographicExtent, GeographicBoundingBox,
+ TemporalExtent, VerticalExtent, Extent, Identifier, PositionalAccuracy,
+
+ - osgeo::proj::common:
+ Common classes: UnitOfMeasure, Measure, Scale, Angle, Length, DateTime,
+ DateEpoch, IdentifiedObject, ObjectDomain, ObjectUsage
+
+ - osgeo::proj::cs:
+ Coordinate systems and their axis
+
+ Classes AxisDirection, Meridian, CoordinateSystemAxis, CoordinateSystem,
+ SphericalCS, EllipsoidalCS, VerticalCS, CartesianCS, OrdinalCS,
+ ParametricCS, TemporalCS, DateTimeTemporalCS, TemporalCountCS,
+ TemporalMeasureCS
+
+ - osgeo::proj::datum:
+ Datum (the relationship of a coordinate system to the body)
+
+ Classes Ellipsoid, PrimeMeridian, Datum, DatumEnsemble,
+ GeodeticReferenceFrame, DynamicGeodeticReferenceFrame,
+ VerticalReferenceFrame, DynamicVerticalReferenceFrame, TemporalDatum,
+ EngineeringDatum, ParametricDatum
+
+ - osgeo::proj::crs:
+ CRS = coordinate reference system = coordinate system with a datum
+
+ Classes CRS, GeodeticCRS, GeographicCRS, DerivedCRS, ProjectedCRS,
+ VerticalCRS, CompoundCRS, BoundCRS, TemporalCRS, EngineeringCRS,
+ ParametricCRS, DerivedGeodeticCRS, DerivedGeographicCRS,
+ DerivedProjectedCRS, DerivedVerticalCRS
+
+ - osgeo::proj::operation:
+ Coordinate operations (relationship between any two coordinate
+ reference systems)
+
+ Classes CoordinateOperation, GeneralOperationParameter,
+ OperationParameter, GeneralParameterValue, ParameterValue,
+ OperationParameterValue, OperationMethod, InvalidOperation,
+ SingleOperation, Conversion, Transformation, PointMotionOperation,
+ ConcatenatedOperation
+
+ - osgeo::proj::io:
+ I/O classes: WKTFormatter, PROJStringFormatter, FormattingException,
+ ParsingException, IWKTExportable, IPROJStringExportable, WKTNode,
+ WKTParser, PROJStringParser, DatabaseContext, AuthorityFactory,
+ FactoryException, NoSuchAuthorityCodeException
+
+What does what?
+***************
+
+The code to parse WKT and PROJ strings and build ISO-19111 objects is
+contained in `io.cpp`_
+
+The code to format WKT and PROJ strings from ISO-19111 objects is mostly
+contained in the related exportToWKT() and exportToPROJString() methods
+overridden in the applicable classes. `io.cpp`_ contains the general mechanics
+to build such strings.
+
+Regarding WKT strings, three variants are handled in import and export:
+
+ - WKT2_2018: variant corresponding to the upcoming ISO-19162:2018 standard
+
+ - WKT2_2015: variant corresponding to the current ISO-19162:2015 standard
+
+ - WKT1_GDAL: variant corresponding to the way GDAL understands the OGC
+ 01-099 and OGC 99-049 standards
+
+Regarding PROJ strings, two variants are handled in import and export:
+
+ - PROJ5: variant used by PROJ >= 5, possibly using pipeline constructs,
+ and avoiding +towgs84 / +nadgrids legacy constructs. This variant honours
+ axis order and input/output units. That is the pipeline for the conversion
+ of EPSG:4326 to EPSG:32631 will assume that the input coordinates are in
+ latitude, longitude order, with degrees.
+
+ - PROJ4: variant used by PROJ 4.x
+
+The raw query of the proj.db database and the upper level construction of
+ISO-19111 objects from the database contents is done in `factory.cpp`_
+
+A few design principles
+***********************
+
+Methods generally take and return xxxNNPtr objects, that is non-null shared
+pointers (pointers with internal reference counting). The advantage of this
+approach is that the user has not to care about the life-cycle of the
+instances (and this makes the code leak-free by design). The only point of
+attention is to make sure no reference cycles are made. This is the case for
+all classes, except the CoordinateOperation class that point to CRS for
+sourceCRS and targetCRS members, whereas DerivedCRS point to a Conversion
+instance (which derives from CoordinateOperation). This issue was detected in
+the ISO-19111 standard. The solution adopted here is to use std::weak_ptr
+in the CoordinateOperation class to avoid the cycle. This design artifact is
+transparent to users.
+
+Another important design point is that all ISO19111 objects are immutable after
+creation, that is they only have getters that do not modify their states.
+Consequently they could possibly use in a thread-safe way. There are however
+classes like PROJStringFormatter, WKTFormatter, DatabaseContext, AuthorityFactory
+and CoordinateOperationContext whose instances are mutable and thus can not be
+used by multiple threads at once.
+
+Example how to build the EPSG:4326 / WGS84 Geographic2D definition from scratch:
+
+::
+
+ auto greenwich = PrimeMeridian::create(
+ util::PropertyMap()
+ .set(metadata::Identifier::CODESPACE_KEY,
+ metadata::Identifier::EPSG)
+ .set(metadata::Identifier::CODE_KEY, 8901)
+ .set(common::IdentifiedObject::NAME_KEY, "Greenwich"),
+ common::Angle(0));
+ // actually predefined as PrimeMeridian::GREENWICH constant
+
+ auto ellipsoid = Ellipsoid::createFlattenedSphere(
+ util::PropertyMap()
+ .set(metadata::Identifier::CODESPACE_KEY, metadata::Identifier::EPSG)
+ .set(metadata::Identifier::CODE_KEY, 7030)
+ .set(common::IdentifiedObject::NAME_KEY, "WGS 84"),
+ common::Length(6378137),
+ common::Scale(298.257223563));
+ // actually predefined as Ellipsoid::WGS84 constant
+
+ auto datum = GeodeticReferenceFrame::create(
+ util::PropertyMap()
+ .set(metadata::Identifier::CODESPACE_KEY, metadata::Identifier::EPSG)
+ .set(metadata::Identifier::CODE_KEY, 6326)
+ .set(common::IdentifiedObject::NAME_KEY, "World Geodetic System 1984");
+ ellipsoid
+ util::optional<std::string>(), // anchor
+ greenwich);
+ // actually predefined as GeodeticReferenceFrame::EPSG_6326 constant
+
+ auto geogCRS = GeographicCRS::create(
+ util::PropertyMap()
+ .set(metadata::Identifier::CODESPACE_KEY, metadata::Identifier::EPSG)
+ .set(metadata::Identifier::CODE_KEY, 4326)
+ .set(common::IdentifiedObject::NAME_KEY, "WGS 84"),
+ datum,
+ cs::EllipsoidalCS::createLatitudeLongitude(scommon::UnitOfMeasure::DEGREE));
+ // actually predefined as GeographicCRS::EPSG_4326 constant
+
+Algorithmic focus
+*****************
+
+On the algorithmic side, a somewhat involved logic is the
+CoordinateOperationFactory::createOperations() in `coordinateoperation.cpp`_
+that takes a pair of source and target CRS and returns a set of possible
+`coordinate operations`_ (either single operations like a Conversion or a
+Transformation, or concatenated operations). It uses the intrinsinc structure
+of those objects to create the coordinate operation pipeline. That is, if
+going from a ProjectedCRS to another one, by doing first the inverse conversion
+from the source ProjectedCRS to its base GeographicCRS, then finding the
+appropriate transformation(s) from this base GeographicCRS to the base
+GeographicCRS of the target CRS, and then appyling the conversion from this
+base GeographicCRS to the target ProjectedCRS. At each step, it queries the
+database to find if one or several transformations are available. The
+resulting coordinate operations are filtered, and sorted, with user provided hints:
+
+ - desired accuracy
+ - area of use, defined as a bounding box in longitude, latitude space (its
+ actual CRS does not matter for the intended use)
+ - if no area of use is defined, if and how the area of use of the source
+ and target CRS should be used. By default, the smallest area of use is
+ used. The rationale is for example when transforming between a national
+ ProjectedCRS and a world-scope GeographicCRS to use the are of use of
+ this ProjectedCRS to select the appropriate datum shifts.
+ - how the area of use of the candidate transformations and the desired area of
+ use (either explicitly or implicitly defined, as explained above) are
+ compared. By default, only transformations whose area of use is fully
+ contained in the desired area of use are selected. It is also possible
+ to relax this test by specifying that only an intersection test must be used.
+ - whether `PROJ transformation grid`_ names should be susbstituted to the
+ official names, when a match is found in the `grid_alternatives` table
+ of the database. Defaults to true
+ - whether the availability of those grids should be used to filter and sort
+ the results. By default, the transformations using grids available in the
+ system will be presented first.
+
+The results are sorted, with the most relevant ones appearing first in the
+result vector. The criteria used are in that order
+
+ - grid actual availability: operations referencing grids not available will be
+ listed after ones with available grids
+ - grid potential availability: operation referencing grids not known at
+ all in the proj.db will be listed after operations with grids known, but
+ not available.
+ - known accuracy: operations with unknown accuracies will be listed
+ after operations with known accuracy
+ - area of use: operations with smaller area of use (the intersection of the
+ operation area of used with the desired area of use) will be listed after
+ the ones with larger area of use
+ - accuracy: operations with lower accuracy will be listed after operations
+ with higher accuracy (caution: lower accuracy actually means a higher numeric
+ value of the accuracy property, since it is a precision in metre)
+
+All those settings can be specified in the CoordinateOperationContext instance
+passed to createOperations().
+
+An interesting example to understand how those parameters play together is
+to use `projinfo -s EPSG:4267 -t EPSG:4326` (NAD27 to WGS84 conversions),
+and see how specifying desired area of use, spatial criterion, grid availability,
+etc. affects the results.
+
+The following command currently returns 78 results:
+
+::
+
+ projinfo -s EPSG:4267 -t EPSG:4326 --summary --spatial-test intersects
+
+The createOperations() algorithm also does a kind of "CRS routing".
+A typical example is if wanting to transform between CRS A and
+CRS B, but no direct transformation is referenced in proj.db between those.
+But if there are transformations between A <--> C and B <--> C, then it
+is possible to build a concatenated operation A --> C --> B. The typical
+example is when C is WGS84, but the implementation is generic and just finds
+a common pivot from the database. An example of finding a non-WGS84 pivot is
+when searching a transformation between EPSG:4326 and EPSG:6668 (JGD2011 -
+Japanese Geodetic Datum 2011), which has no direct transformation registered
+in the EPSG database . However there are transformations between those two
+CRS and JGD2000 (and also Tokyo datum, but that one involves less accurate
+transformations)
+
+::
+
+ projinfo -s EPSG:4326 -t EPSG:6668 --grid-check none --bbox 135.42,34.84,142.14,41.58 --summary
+
+ Candidate operations found: 7
+ unknown id, Inverse of JGD2000 to WGS 84 (1) + JGD2000 to JGD2011 (1), 1.2 m, Japan - northern Honshu
+ unknown id, Inverse of JGD2000 to WGS 84 (1) + JGD2000 to JGD2011 (2), 2 m, Japan excluding northern main province
+ unknown id, Inverse of Tokyo to WGS 84 (108) + Tokyo to JGD2011 (2), 9.2 m, Japan onshore excluding northern main province
+ unknown id, Inverse of Tokyo to WGS 84 (108) + Tokyo to JGD2000 (2) + JGD2000 to JGD2011 (1), 9.4 m, Japan - northern Honshu
+ unknown id, Inverse of Tokyo to WGS 84 (2) + Tokyo to JGD2011 (2), 13.2 m, Japan - onshore mainland and adjacent islands
+ unknown id, Inverse of Tokyo to WGS 84 (2) + Tokyo to JGD2000 (2) + JGD2000 to JGD2011 (1), 13.4 m, Japan - northern Honshu
+ unknown id, Inverse of Tokyo to WGS 84 (1) + Tokyo to JGD2011 (2), 29.2 m, Asia - Japan and South Korea
+
+
+.. _`coordinate operations`: https://proj4.org/operations/index.html
+
+.. _`PROJ transformation grid`: https://proj4.org/resource_files.html#transformation-grids
+
+
+Code repository
+---------------
+
+The current state of the work can be found in the
+`iso19111 branch of rouault/proj.4 repository`_ , and is also available as
+a GitHub pull request at https://github.com/OSGeo/proj.4/pull/1040
+
+Here is a not-so-usable `comparison with a fixed snapshot of master branch`_
+
+.. _`iso19111 branch of rouault/proj.4 repository`: https://github.com/rouault/proj.4/tree/iso19111
+
+.. _`comparison with a fixed snapshot of master branch`: https://github.com/OSGeo/proj.4/compare/iso19111_dev...rouault:iso19111
+
+
+Database
+--------
+
+Content
+*******
+
+The database contains CRS and coordinate operation definitions from the `EPSG`_
+database (IOGP’s EPSG Geodetic Parameter Dataset) v9.5.3,
+`IGNF registry`_ (French National Geographic Institute), ESRI database, as well
+as a few customizations.
+
+.. _`EPSG`: http://www.epsg.org/
+.. _`IGNF registry`: https://geodesie.ign.fr/index.php?page=documentation#titre3
+
+Building (for PROJ developers creating the database)
+****************************************************
+
+The building of the database is a several stage process:
+
+Construct SQL scripts for EPSG
+++++++++++++++++++++++++++++++
+
+The first stage consists in constructing .sql scripts mostly with
+CREATE TABLE and INSERT statements to create the database structure and
+populate it. There is one .sql file for each database table, populated
+with the content of the EPSG database, automatically
+generated with the `build_db.py`_ script, which processes the PostgreSQL
+dumps issued by IOGP. A number of other scripts are dedicated to manual
+editing, for example `grid_alternatives.sql`_ file that binds official
+grid names to PROJ grid names
+
+Concert UTF8 SQL to sqlite3 db
+++++++++++++++++++++++++++++++
+
+The second stage is done automatically by the make process. It pipes the
+.sql script, in the right order, to the sqlite3 binary to generate a
+first version of the proj.db SQLite3 database.
+
+Add extra registries
+++++++++++++++++++++
+
+The third stage consists in creating additional .sql files from the
+content of other registries. For that process, we need to bind some
+definitions of those registries to those of the EPSG database, to be
+able to link to existing objects and detect some boring duplicates.
+The `ignf.sql`_ file has been generated using
+the `build_db_create_ignf.py`_ script from the current data/IGNF file
+that contains CRS definitions (and implicit transformations to WGS84)
+as PROJ.4 strings.
+The `esri.sql`_ file has been generated using the `build_db_from_esri.py`_
+script, from the .csv files in
+https://github.com/Esri/projection-engine-db-doc/tree/master/csv
+
+Finalize proj.db
+++++++++++++++++
+
+The last stage runs make again to incorporate the new .sql files generated
+in the previous stage (so the process of building the database involves
+a kind of bootstrapping...)
+
+Building (for PROJ users)
+*************************
+
+The make process just runs the second stage mentionned above from the .sql
+files. The resulting proj.db is currently 5.3 MB large.
+
+Structure
+*********
+
+The database is structured into the following tables and views. They generally
+match a ISO-19111 concept, and is generally close to the general structure of
+the EPSG database. Regarding identification of objects, where the EPSG database
+only contains a 'code' numeric column, the PROJ database identifies objects
+with a (auth_name, code) tuple of string values, allowing several registries to
+be combined together.
+
+- Technical:
+ - `authority_list`: view enumerating the authorities present in the database. Currently: EPSG, IGNF, PROJ
+ - `metadata`: a few key/value pairs, for example to indicate the version of the registries imported in the database
+ - `object_view`: synthetic view listing objects (ellipsoids, datums, CRS, coordinate operations...) code and name, and the table name where they are further described
+ - `alias_names`: list possible alias for the `name` field of object table
+ - `link_from_deprecated_to_non_deprecated`: to handle the link between old ESRI to new ESRI/EPSG codes
+
+- Commmon:
+ - `unit_of_measure`: table with UnitOfMeasure definitions.
+ - `area`: table with area-of-use (bounding boxes) applicable to CRS and coordinate operations.
+
+- Coordinate systems:
+ - `axis`: table with CoordinateSystemAxis definitions.
+ - `coordinate_system`: table with CoordinateSystem definitions.
+
+- Ellipsoid and datums:
+ - `ellipsoid`: table with ellipsoid definitions.
+ - `prime_meridian`: table with PrimeMeridian definitions.
+ - `geodetic_datum`: table with GeodeticReferenceFrame definitions.
+ - `vertical_datum`: table with VerticalReferenceFrame definitions.
+
+- CRS:
+ - `geodetic_crs`: table with GeodeticCRS and GeographicCRS definitions.
+ - `projected_crs`: table with ProjectedCRS definitions.
+ - `vertical_crs`: table with VerticalCRS definitions.
+ - `compound_crs`: table with CompoundCRS definitions.
+
+- Coordinate operations:
+ - `coordinate_operation_view`: view giving a number of common attributes shared by the concrete tables implementing CoordinateOperation
+ - `conversion`: table with definitions of Conversion (mostly parameter and values of Projection)
+ - `concatenated_operation`: table with definitions of ConcatenatedOperation.
+ - `grid_transformation`: table with all grid-based transformations.
+ - `grid_packages`: table listing packages in which grids can be found. ie "proj-datumgrid", "proj-datumgrid-europe", ...
+ - `grid_alternatives`: table binding official grid names to PROJ grid names. e.g "Und_min2.5x2.5_egm2008_isw=82_WGS84_TideFree.gz" --> "egm08_25.gtx"
+ - `helmert_transformation`: table with all Helmert-based transformations.
+ - `other_transformation`: table with other type of transformations.
+
+The main departure with the structure of the EPSG database is the split of
+the various coordinate operations over several tables. This was done mostly
+for human-readability as the EPSG organization of coordoperation,
+coordoperationmethod, coordoperationparam, coordoperationparamusage,
+coordoperationparamvalue tables makes it hard to grasp at once all the
+parameters and values for a given operation.
+
+
+Utilities
+---------
+
+A new `projinfo` utility has been added. It enables the user to enter a CRS or
+coordinate operation by a AUTHORITY:CODE, PROJ string or WKT string, and see
+it translated in the different flavors of PROJ and WKT strings.
+It also enables to build coordinate operations between two CRSs.
+
+Usage
+*****
+
+::
+
+ usage: projinfo [-o formats] [-k crs|operation] [--summary] [-q]
+ [--bbox min_long,min_lat,max_long,max_lat]
+ [--spatial-test contains|intersects]
+ [--crs-extent-use none|both|intersection|smallest]
+ [--grid-check none|discard_missing|sort]
+ [--boundcrs-to-wgs84]
+ {object_definition} | (-s {srs_def} -t {srs_def})
+
+ -o: formats is a comma separated combination of: all,default,PROJ4,PROJ,WKT_ALL,WKT2_2015,WKT2_2018,WKT1_GDAL
+ Except 'all' and 'default', other format can be preceded by '-' to disable them
+
+Examples
+********
+
+Specify CRS by AUTHORITY:CODE
++++++++++++++++++++++++++++++
+
+::
+
+ $ projinfo EPSG:4326
+
+ PROJ string:
+ +proj=pipeline +step +proj=longlat +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
+
+ WKT2_2015 string:
+ GEODCRS["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]],
+ AREA["World"],
+ BBOX[-90,-180,90,180],
+ ID["EPSG",4326]]
+
+
+Specify CRS by PROJ string and specify output formats
+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+::
+
+ $ projinfo -o PROJ4,PROJ,WKT1_GDAL,WKT2_2018 "+title=IGN 1972 Nuku Hiva - UTM fuseau 7 Sud +proj=tmerc +towgs84=165.7320,216.7200,180.5050,-0.6434,-0.4512,-0.0791,7.420400 +a=6378388.0000 +rf=297.0000000000000 +lat_0=0.000000000 +lon_0=-141.000000000 +k_0=0.99960000 +x_0=500000.000 +y_0=10000000.000 +units=m +no_defs"
+
+ PROJ string:
+ Error when exporting to PROJ string: BoundCRS cannot be exported as a PROJ.5 string, but its baseCRS might
+
+ PROJ.4 string:
+ +proj=utm +zone=7 +south +ellps=intl +towgs84=165.732,216.72,180.505,-0.6434,-0.4512,-0.0791,7.4204
+
+ WKT2_2018 string:
+ BOUNDCRS[
+ SOURCECRS[
+ PROJCRS["IGN 1972 Nuku Hiva - UTM fuseau 7 Sud",
+ BASEGEOGCRS["unknown",
+ DATUM["unknown",
+ ELLIPSOID["International 1909 (Hayford)",6378388,297,
+ LENGTHUNIT["metre",1,
+ ID["EPSG",9001]]]],
+ PRIMEM["Greenwich",0,
+ ANGLEUNIT["degree",0.0174532925199433],
+ ID["EPSG",8901]]],
+ CONVERSION["unknown",
+ METHOD["Transverse Mercator",
+ ID["EPSG",9807]],
+ PARAMETER["Latitude of natural origin",0,
+ ANGLEUNIT["degree",0.0174532925199433],
+ ID["EPSG",8801]],
+ PARAMETER["Longitude of natural origin",-141,
+ ANGLEUNIT["degree",0.0174532925199433],
+ ID["EPSG",8802]],
+ PARAMETER["Scale factor at natural origin",0.9996,
+ SCALEUNIT["unity",1],
+ ID["EPSG",8805]],
+ PARAMETER["False easting",500000,
+ LENGTHUNIT["metre",1],
+ ID["EPSG",8806]],
+ PARAMETER["False northing",10000000,
+ LENGTHUNIT["metre",1],
+ ID["EPSG",8807]]],
+ CS[Cartesian,2],
+ AXIS["(E)",east,
+ ORDER[1],
+ LENGTHUNIT["metre",1,
+ ID["EPSG",9001]]],
+ AXIS["(N)",north,
+ ORDER[2],
+ LENGTHUNIT["metre",1,
+ ID["EPSG",9001]]]]],
+ TARGETCRS[
+ 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["latitude",north,
+ ORDER[1],
+ ANGLEUNIT["degree",0.0174532925199433]],
+ AXIS["longitude",east,
+ ORDER[2],
+ ANGLEUNIT["degree",0.0174532925199433]],
+ ID["EPSG",4326]]],
+ ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
+ METHOD["Position Vector transformation (geog2D domain)",
+ ID["EPSG",9606]],
+ PARAMETER["X-axis translation",165.732,
+ ID["EPSG",8605]],
+ PARAMETER["Y-axis translation",216.72,
+ ID["EPSG",8606]],
+ PARAMETER["Z-axis translation",180.505,
+ ID["EPSG",8607]],
+ PARAMETER["X-axis rotation",-0.6434,
+ ID["EPSG",8608]],
+ PARAMETER["Y-axis rotation",-0.4512,
+ ID["EPSG",8609]],
+ PARAMETER["Z-axis rotation",-0.0791,
+ ID["EPSG",8610]],
+ PARAMETER["Scale difference",1.0000074204,
+ ID["EPSG",8611]]]]
+
+ WKT1_GDAL:
+ PROJCS["IGN 1972 Nuku Hiva - UTM fuseau 7 Sud",
+ GEOGCS["unknown",
+ DATUM["unknown",
+ SPHEROID["International 1909 (Hayford)",6378388,297],
+ TOWGS84[165.732,216.72,180.505,-0.6434,-0.4512,-0.0791,7.4204]],
+ PRIMEM["Greenwich",0,
+ AUTHORITY["EPSG","8901"]],
+ UNIT["degree",0.0174532925199433,
+ AUTHORITY["EPSG","9122"]],
+ AXIS["Longitude",EAST],
+ AXIS["Latitude",NORTH]],
+ PROJECTION["Transverse_Mercator"],
+ PARAMETER["latitude_of_origin",0],
+ PARAMETER["central_meridian",-141],
+ PARAMETER["scale_factor",0.9996],
+ PARAMETER["false_easting",500000],
+ PARAMETER["false_northing",10000000],
+ UNIT["metre",1,
+ AUTHORITY["EPSG","9001"]],
+ AXIS["Easting",EAST],
+ AXIS["Northing",NORTH]]
+
+
+Find transformations between 2 CRS
+++++++++++++++++++++++++++++++++++
+
+Between "Poland zone I" (based on Pulkovo 42 datum) and "UTM WGS84 zone 34N"
+
+Summary view:
+
+::
+
+ $ projinfo -s EPSG:2171 -t EPSG:32634 --summary
+
+ Candidate operations found: 1
+ unknown id, Inverse of Poland zone I + Pulkovo 1942(58) to WGS 84 (1) + UTM zone 34N, 1 m, Poland - onshore
+
+Display of pipelines:
+
+::
+
+ $ PROJ_LIB=data src/projinfo -s EPSG:2171 -t EPSG:32634 -o PROJ
+
+ PROJ string:
+ +proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=sterea +lat_0=50.625 +lon_0=21.0833333333333 +k=0.9998 +x_0=4637000 +y_0=5647000 +ellps=krass +step +proj=cart +ellps=krass +step +proj=helmert +x=33.4 +y=-146.6 +z=-76.3 +rx=-0.359 +ry=-0.053 +rz=0.844 +s=-0.84 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=utm +zone=34 +ellps=WGS84
+
+
+Impacted files
+--------------
+
+New files (excluding makefile.am, CMakeLists.txt and other build infrastructure
+artifacts):
+
+ * include/proj/: Public installed C++ headers
+ - `common.hpp`_: declarations of osgeo::proj::common namespace.
+ - `coordinateoperation.hpp`_: declarations of osgeo::proj::operation namespace.
+ - `coordinatesystem.hpp`_: declarations of osgeo::proj::cs namespace.
+ - `crs.hpp`_: declarations of osgeo::proj::crs namespace.
+ - `datum.hpp`_: declarations of osgeo::proj::datum namespace.
+ - `io.hpp`_: declarations of osgeo::proj::io namespace.
+ - `metadata.hpp`_: declarations of osgeo::proj::metadata namespace.
+ - `util.hpp`_: declarations of osgeo::proj::util namespace.
+ - `nn.hpp`_: Code from https://github.com/dropbox/nn to manage Non-nullable pointers for C++
+
+ .. _`common.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/common.hpp
+ .. _`coordinateoperation.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/coordinateoperation.hpp
+ .. _`coordinatesystem.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/coordinatesystem.hpp
+ .. _`crs.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/crs.hpp
+ .. _`datum.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/datum.hpp
+ .. _`io.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/io.hpp
+ .. _`metadata.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/metadata.hpp
+ .. _`nn.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/nn.hpp
+ .. _`util.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/util.hpp
+
+ * include/proj/internal: Private non-installed C++ headers
+ - `coordinateoperation_internal.hpp`_: classes InverseCoordinateOperation, InverseConversion, InverseTransformation, PROJBasedOperation, and functions to get conversion mappings between WKT and PROJ syntax
+ - `coordinateoperation_constants.hpp`_: Select subset of conversion/transformation EPSG names and codes for the purpose of translating them to PROJ strings
+ - `coordinatesystem_internal.hpp`_: classes AxisDirectionWKT1, AxisName and AxisAbbreviation
+ - `internal.hpp`_: a few helper functions, mostly to do string-based operations
+ - `io_internal.hpp`_: class WKTConstants
+ - `helmert_constants.hpp`_: Helmert-based transformation & parameters names and codes.
+ - `lru_cache.hpp`_: code from https://github.com/mohaps/lrucache11 to have a generic Least-Recently-Used cache of objects
+
+ .. _`coordinateoperation_internal.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/coordinateoperation_internal.hpp
+ .. _`coordinatesystem_internal.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/coordinatesystem_internal.hpp
+ .. _`internal.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/internal.hpp
+ .. _`io_internal.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/io_internal.hpp
+ .. _`coordinateoperation_constants.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/coordinateoperation_constants.hpp
+ .. _`helmert_constants.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/helmert_constants.hpp
+ .. _`lru_cache.hpp`: https://github.com/rouault/proj.4/blob/iso19111/include/proj/internal/lru_cache.hpp
+
+ * src/:
+ - `c_api.cpp`_: C++ API mapped to C functions
+ - `common.cpp`_: implementation of `common.hpp`_
+ - `coordinateoperation.cpp`_: implementation of `coordinateoperation.hpp`_
+ - `coordinatesystem.cpp`_: implementation of `coordinatesystem.hpp`_
+ - `crs.cpp`_: implementation of `crs.hpp`_
+ - `datum.cpp`_: implementation of `datum.hpp`_
+ - `factory.cpp`_: implementation of AuthorityFactory class (from `io.hpp`_)
+ - `internal.cpp`_: implementation of `internal.hpp`_
+ - `io.cpp`_: implementation of `io.hpp`_
+ - `metadata.cpp`_: implementation of `metadata.hpp`_
+ - `static.cpp`_: a number of static constants (like pre-defined well-known ellipsoid, datum and CRS), put in the right order for correct static initializations
+ - `util.cpp`_: implementation of `util.hpp`_
+ - `projinfo.cpp`_: new 'projinfo' binary
+ - `general.dox`_: generic introduction documentation.
+
+ .. _`c_api.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/c_api.cpp
+ .. _`common.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/common.cpp
+ .. _`coordinateoperation.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/coordinateoperation.cpp
+ .. _`coordinatesystem.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/coordinatesystem.cpp
+ .. _`crs.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/crs.cpp
+ .. _`datum.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/datum.cpp
+ .. _`factory.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/factory.cpp
+ .. _`internal.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/internal.cpp
+ .. _`io.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/io.cpp
+ .. _`metadata.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/metadata.cpp
+ .. _`projinfo.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/projinfo.cpp
+ .. _`static.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/static.cpp
+ .. _`util.cpp`: https://github.com/rouault/proj.4/blob/iso19111/src/util.cpp
+ .. _`general.dox`: https://github.com/rouault/proj.4/blob/iso19111/src/general.dox
+
+ * data/sql/:
+ - `area.sql`_: generated by `build_db.py`_
+ - `axis.sql`_: generated by `build_db.py`_
+ - `begin.sql`_: hand generated (trivial)
+ - `commit.sql`_: hand generated (trivial)
+ - `compound_crs.sql`_: generated by `build_db.py`_
+ - `concatenated_operation.sql`_: generated by `build_db.py`_
+ - `conversion.sql`_: generated by `build_db.py`_
+ - `coordinate_operation.sql`_: generated by `build_db.py`_
+ - `coordinate_system.sql`_: generated by `build_db.py`_
+ - `crs.sql`_: generated by `build_db.py`_
+ - `customizations.sql`_: hand generated (empty)
+ - `ellipsoid.sql`_: generated by `build_db.py`_
+ - `geodetic_crs.sql`_: generated by `build_db.py`_
+ - `geodetic_datum.sql`_: generated by `build_db.py`_
+ - `grid_alternatives.sql`_: hand-generated. Contains links between official registry grid names and PROJ ones
+ - `grid_transformation.sql`_: generated by `build_db.py`_
+ - `grid_transformation_custom.sql`_: hand-generated
+ - `helmert_transformation.sql`_: generated by `build_db.py`_
+ - `ignf.sql`_: generated by `build_db_create_ignf.py`_
+ - `esri.sql`_: generated by `build_db_from_esri.py`_
+ - `metadata.sql`_: hand-generated
+ - `other_transformation.sql`_: generated by `build_db.py`_
+ - `prime_meridian.sql`_: generated by `build_db.py`_
+ - `proj_db_table_defs.sql`_: hand-generated. Database structure: CREATE TABLE / CREATE VIEW / CREATE TRIGGER
+ - `projected_crs.sql`_: generated by `build_db.py`_
+ - `unit_of_measure.sql`_: generated by `build_db.py`_
+ - `vertical_crs.sql`_: generated by `build_db.py`_
+ - `vertical_datum.sql`_: generated by `build_db.py`_
+
+ .. _`area.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/area.sql
+ .. _`axis.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/axis.sql
+ .. _`begin.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/begin.sql
+ .. _`commit.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/commit.sql
+ .. _`compound_crs.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/compound_crs.sql
+ .. _`concatenated_operation.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/concatenated_operation.sql
+ .. _`conversion.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/conversion.sql
+ .. _`coordinate_operation.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/coordinate_operation.sql
+ .. _`coordinate_system.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/coordinate_system.sql
+ .. _`crs.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/crs.sql
+ .. _`customizations.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/customizations.sql
+ .. _`ellipsoid.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/ellipsoid.sql
+ .. _`geodetic_crs.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/geodetic_crs.sql
+ .. _`geodetic_datum.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/geodetic_datum.sql
+ .. _`grid_alternatives.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/grid_alternatives.sql
+ .. _`grid_transformation_custom.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/grid_transformation_custom.sql
+ .. _`grid_transformation.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/grid_transformation.sql
+ .. _`helmert_transformation.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/helmert_transformation.sql
+ .. _`ignf.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/ignf.sql
+ .. _`esri.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/esri.sql
+ .. _`metadata.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/metadata.sql
+ .. _`other_transformation.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/other_transformation.sql
+ .. _`prime_meridian.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/prime_meridian.sql
+ .. _`proj_db_table_defs.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/proj_db_table_defs.sql
+ .. _`projected_crs.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/projected_crs.sql
+ .. _`unit_of_measure.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/unit_of_measure.sql
+ .. _`vertical_crs.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/vertical_crs.sql
+ .. _`vertical_datum.sql`: https://github.com/rouault/proj.4/blob/iso19111/data/sql/vertical_datum.sql
+
+ * scripts/:
+ - `build_db.py`_ : generate .sql files from EPSG database dumps
+ - `build_db_create_ignf.py`_: generates data/sql/`ignf.sql`_
+ - `build_db_from_esri.py`_: generates data/sql/`esri.sql`_
+ - `doxygen.sh`_: generates Doxygen documentation
+ - `gen_html_coverage.sh`_: generates HTML report of the coverage for --coverage build
+ - `filter_lcov_info.py`_: utility used by gen_html_coverage.sh
+ - `reformat.sh`_: used by reformat_cpp.sh
+ - `reformat_cpp.sh`_: reformat all .cpp/.hpp files according to LLVM-style formatting rules
+
+ .. _`build_db.py`: https://github.com/rouault/proj.4/blob/iso19111/scripts/build_db.py
+ .. _`build_db_create_ignf.py`: https://github.com/rouault/proj.4/blob/iso19111/scripts/build_db_create_ignf.py
+ .. _`build_db_from_esri.py`: https://github.com/rouault/proj.4/blob/iso19111/scripts/build_db_from_esri.py
+ .. _`doxygen.sh`: https://github.com/rouault/proj.4/blob/iso19111/scripts/doxygen.sh
+ .. _`gen_html_coverage.sh`: https://github.com/rouault/proj.4/blob/iso19111/scripts/gen_html_coverage.sh
+ .. _`filter_lcov_info.py`: https://github.com/rouault/proj.4/blob/iso19111/scripts/filter_lcov_info.py
+ .. _`reformat.sh`: https://github.com/rouault/proj.4/blob/iso19111/scripts/reformat.sh
+ .. _`reformat_cpp.sh`: https://github.com/rouault/proj.4/blob/iso19111/scripts/reformat_cpp.sh
+
+ * tests/unit/
+ - `test_c_api.cpp`_: test of src/c_api.cpp
+ - `test_common.cpp`_: test of src/common.cpp
+ - `test_util.cpp`_: test of src/util.cpp
+ - `test_crs.cpp`_: test of src/crs.cpp
+ - `test_datum.cpp`_: test of src/datum.cpp
+ - `test_factory.cpp`_: test of src/factory.cpp
+ - `test_io.cpp`_: test of src/io.cpp
+ - `test_metadata.cpp`_: test of src/metadata.cpp
+ - `test_operation.cpp`_: test of src/operation.cpp
+
+ .. _`test_c_api.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_c_api.cpp
+ .. _`test_common.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_common.cpp
+ .. _`test_util.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_util.cpp
+ .. _`test_crs.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_crs.cpp
+ .. _`test_datum.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_datum.cpp
+ .. _`test_factory.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_factory.cpp
+ .. _`test_io.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_io.cpp
+ .. _`test_metadata.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_metadata.cpp
+ .. _`test_operation.cpp`: https://github.com/rouault/proj.4/blob/iso19111/test/unit/test_operation.cpp
+
+C API
+-----
+
+`proj.h`_ has been extended to bind a number of C++ classes/methods to a C API.
+
+The main structure is an opaque PJ_OBJ* roughly encapsulating a osgeo::proj::BaseObject,
+that can represent a CRS or a CoordinateOperation object. A number of the
+C functions will work only if the right type of underlying C++ object is used
+with them. Misuse will be properly handled at runtime. If a user passes
+a PJ_OBJ* representing a coordinate operation to a pj_obj_crs_xxxx() function,
+it will properly error out. This design has been chosen over creating a
+dedicate PJ_xxx object for each C++ class, because such an approach would
+require adding many conversion and free functions for little benefit.
+
+This C API is incomplete. In particular, it does not allow to
+build ISO19111 objects at hand. However it currently permits a number of
+actions:
+
+ - building CRS and coordinate operations from WKT and PROJ strings, or
+ from the proj.db database
+ - exporting CRS and coordinate operations as WKT and PROJ strings
+ - querying main attributes of those objects
+ - finding coordinate operations between two CRS.
+
+`test_c_api.cpp`_ should demonstrates simple usage of the API (note:
+for the conveniency of writing the tests in C++, test_c_api.cpp wraps the C PJ_OBJ*
+instances in C++ 'keeper' objects that automatically call the pj_obj_unref()
+function at function end. In a pure C use, the caller must use pj_obj_unref()
+to prevent leaks.)
+
+.. _`proj.h`: http://even.rouault.free.fr/proj_cpp_api/html/proj_8h.html
+
+
+Documentation
+-------------
+
+All public C++ classes and methods and C functions are documented with
+Doxygen.
+
+`Current snapshot of Class list`_
+
+`Spaghetti inheritance diagram`_
+
+.. _`Current snapshot of Class list`: http://even.rouault.free.fr/proj_cpp_api/html/annotated.html
+.. _`Spaghetti inheritance diagram`: http://even.rouault.free.fr/proj_cpp_api/html/inherits.html
+
+A basic integration of the Doxygen XML output into the general PROJ
+documentation (using reStructuredText format) has been done with the the
+Sphinx extension `Breathe`_, producing:
+
+ * `One section with the C++ API`_
+ * `One section with the C API`_
+
+.. _`Breathe`: https://breathe.readthedocs.io/en/latest/
+.. _`One section with the C++ API`: http://even.rouault.free.fr/proj_cpp_api/rst_generated/html/development/reference/cpp/index.html
+.. _`One section with the C API`: http://even.rouault.free.fr/proj_cpp_api/rst_generated/html/development/reference/functions.html#c-api-for-iso-19111-functionality
+
+Testing
+-------
+
+Nearly all exported methods are tested by a unit test. Global line coverage
+of the new files is 92%. Those tests represent 16k lines of codes.
+
+
+Build requirements
+------------------
+
+The new code leverages on a number of C++11 features (auto keyword, constexpr,
+initializer list, std::shared_ptr, lambda functions, etc.), which means that
+a C++11-compliant compiler must be used to generate PROJ:
+
+ * gcc >= 4.8
+ * clang >= 3.3
+ * Visual Studio >= 2015.
+
+Compilers tested by the Travis-CI and AppVeyor continuous integration
+environments:
+
+ * GCC 4.8
+ * mingw-w64-x86-64 4.8
+ * clang 5.0
+ * Apple LLVM version 9.1.0 (clang-902.0.39.2)
+ * MSVC 2015 32 and 64 bit
+ * MSVC 2017 32 and 64 bit
+
+The libsqlite3 >= 3.7 development package must also be available. And the sqlite3
+binary must be available to build the proj.db files from the .sql files.
+
+Runtime requirements
+--------------------
+
+* libc++/libstdc++/MSVC runtime consistent with the compiler used
+* libsqlite3 >= 3.7
+
+
+Backward compatibility
+----------------------
+
+At this stage, no backward compatibility issue is foreseen, as no
+existing functional C code has been modified to use the new capabilities
+
+Future work
+-----------
+
+The work described in this RFC will be pursued in a number of directions.
+Non-exhaustively:
+
+ - Support for ESRI WKT1 dialect (PROJ currently ingest the ProjectedCRS in
+ `esri.sql`_ in that dialect, but there is no mapping between it and EPSG
+ operation and parameter names, so conversion to PROJ strings does not
+ always work.
+
+ - closer integration with the existing code base. In particular, the +init=dict:code
+ syntax should now go first to the database (then the `epsg` and `IGNF`
+ files can be removed). Similarly proj_create_crs_to_crs() could use the
+ new capabilities to find an appropriate coordinate transformation.
+
+ - and whatever else changes are needed to address GDAL and libgeotiff needs
+
+
+Adoption status
+---------------
+
+The RFC has been adopted with support from PSC members Kurt Schwehr, Kristian
+Evers, Howard Butler and Even Rouault.
diff --git a/docs/source/development/migration.rst b/docs/source/development/migration.rst
index 34e166bf..99143f20 100644
--- a/docs/source/development/migration.rst
+++ b/docs/source/development/migration.rst
@@ -153,8 +153,8 @@ Function mapping from old to new API
###############################################################################
+---------------------------------------+---------------------------------------+
-| **Old API functions** | **New API functions** |
-+---------------------------------------+---------------------------------------+
+| Old API functions | New API functions |
++=======================================+=======================================+
| pj_fwd | proj_trans |
+---------------------------------------+---------------------------------------+
| pj_inv | proj_trans |
diff --git a/docs/source/operations/conversions/unitconvert.rst b/docs/source/operations/conversions/unitconvert.rst
index e99ca837..f1488d9a 100644
--- a/docs/source/operations/conversions/unitconvert.rst
+++ b/docs/source/operations/conversions/unitconvert.rst
@@ -87,8 +87,8 @@ The same list can also be produced on the command line with :program:`proj` or
:program:`cs2cs`, by adding the `-lu` flag when calling the utility.
+----------+---------------------------------+
-| **Label**| **Name** |
-+----------+---------------------------------+
+| Label | Name |
++==========+=================================+
| km | Kilometer |
+----------+---------------------------------+
| m | Meter |
@@ -142,8 +142,8 @@ Angular units
In the table below all angular units supported by PROJ `unitconvert` are listed.
+----------+---------------------------------+
-| **Label**| **Name** |
-+----------+---------------------------------+
+| Label | Name |
++==========+=================================+
| deg | Degree |
+----------+---------------------------------+
| grad | Grad |
@@ -159,8 +159,8 @@ Time units
In the table below all time units supported by PROJ are listed.
+--------------+-----------------------------+
-| **label** | **Name** |
-+--------------+-----------------------------+
+| Label | Name |
++==============+=============================+
| mjd | Modified Julian date |
+--------------+-----------------------------+
| decimalyear | Decimal year |
diff --git a/docs/source/operations/projections/eqc.rst b/docs/source/operations/projections/eqc.rst
index e5c57f47..d36c1854 100644
--- a/docs/source/operations/projections/eqc.rst
+++ b/docs/source/operations/projections/eqc.rst
@@ -40,7 +40,7 @@ The following table gives special cases of the cylindrical equidistant projectio
+---------------------------------------------------------+--------------------------+
| Projection Name | (lat ts=) :math:`\phi_0` |
-+---------------------------------------------------------+--------------------------+
++=========================================================+==========================+
| Plain/Plane Chart | 0° |
+---------------------------------------------------------+--------------------------+
| Simple Cylindrical | 0° |
diff --git a/docs/source/operations/projections/tmerc.rst b/docs/source/operations/projections/tmerc.rst
index 4c9cdab0..6d1c145d 100644
--- a/docs/source/operations/projections/tmerc.rst
+++ b/docs/source/operations/projections/tmerc.rst
@@ -44,7 +44,7 @@ The following table gives special cases of the Transverse Mercator projection.
+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
| Projection Name | Areas | Central meridian | Zone width | Scale Factor |
-+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
++=====================================+=====================================================+================================+==========================================+==============+
| Transverse Mercator | World wide | Various | less than 6° | Various |
+-------------------------------------+-----------------------------------------------------+--------------------------------+------------------------------------------+--------------+
| Transverse Mercator south oriented | Southern Africa | 2° intervals E of 11°E | 2° | 1.000 |
diff --git a/docs/source/operations/transformations/helmert.rst b/docs/source/operations/transformations/helmert.rst
index 85e140af..73449b40 100644
--- a/docs/source/operations/transformations/helmert.rst
+++ b/docs/source/operations/transformations/helmert.rst
@@ -7,7 +7,7 @@ Helmert transform
.. versionadded:: 5.0.0
The Helmert transformation changes coordinates from one reference frame to
-anoether by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
+another by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
14-parameter kinematic counterparts.
@@ -25,7 +25,7 @@ anoether by means of 3-, 4-and 7-parameter shifts, or one of their 6-, 8- and
| **Output type** | Cartesian coordinates |
+-----------------+-------------------------------------------------------------------+
-The Helmert transform, in all it's various incarnations, is used to perform reference
+The Helmert transform, in all its various incarnations, is used to perform reference
frame shifts. The transformation operates in cartesian space. It can be used to transform
planar coordinates from one datum to another, transform 3D cartesian
coordinates from one static reference frame to another or it can be used to do fully
@@ -117,7 +117,7 @@ Parameters
.. option:: +y=<value>
- Translation of the x-axis given in meters.
+ Translation of the y-axis given in meters.
.. option:: +z=<value>
diff --git a/docs/source/operations/transformations/index.rst b/docs/source/operations/transformations/index.rst
index 45c4a80b..d4bbf4e0 100644
--- a/docs/source/operations/transformations/index.rst
+++ b/docs/source/operations/transformations/index.rst
@@ -16,5 +16,6 @@ systems are based on different datums.
helmert
horner
molodensky
+ molobadekas
hgridshift
vgridshift
diff --git a/docs/source/operations/transformations/molobadekas.rst b/docs/source/operations/transformations/molobadekas.rst
new file mode 100644
index 00000000..b7d638bf
--- /dev/null
+++ b/docs/source/operations/transformations/molobadekas.rst
@@ -0,0 +1,147 @@
+.. _molobadekas:
+
+================================================================================
+Molodensky-Badekas transform
+================================================================================
+
+.. versionadded:: 6.0.0
+
+The Molodensky-Badekas transformation changes coordinates from one reference frame to
+another by means of a 10-parameter shift.
+
+.. note::
+
+ It should not be confused with the :ref:`Molodensky` transform which
+ operates directly in the geodetic coordinates. Molodensky-Badekas can rather
+ be seen as a variation of :ref:`Helmert`
+
++-----------------+-------------------------------------------------------------------+
+| **Alias** | molobadekas |
++-----------------+-------------------------------------------------------------------+
+| **Domain** | 3D |
++-----------------+-------------------------------------------------------------------+
+| **Input type** | Cartesian coordinates |
++-----------------+-------------------------------------------------------------------+
+| **Output type** | Cartesian coordinates |
++-----------------+-------------------------------------------------------------------+
+
+The Molodensky-Badekas transformation is a variation of the :ref:`Helmert` where
+the rotational terms are not directly applied to the ECEF coordinates, but on
+cartesian coordinates relative to a reference point (usually close to Earth surface,
+and to the area of use of the transformation). When ``px`` = ``py`` = ``pz`` = 0,
+this is equivalent to a 7-parameter Helmert transformation.
+
+Example
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+Transforming coordinates from La Canoa to REGVEN:
+
+::
+
+ proj=molobadekas convention=coordinate_frame
+ x=-270.933 y=115.599 z=-360.226 rx=-5.266 ry=-1.238 rz=2.381
+ s=-5.109 px=2464351.59 py=-5783466.61 pz=974809.81
+
+
+Parameters
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+.. note::
+
+ All parameters (except convention) are optional but at least one should be
+ used, otherwise the operation will return the coordinates unchanged.
+
+.. option:: +convention=coordinate_frame/position_vector
+
+ Indicates the convention to express the rotational terms when a 3D-Helmert /
+ 7-parameter more transform is involved.
+
+ The two conventions are equally popular and a frequent source of confusion.
+ The coordinate frame convention is also described as an clockwise
+ rotation of the coordinate frame. It corresponds to EPSG method code
+ 1034 (in the geocentric domain) or 9636 (in the geographic domain)
+ The position vector convention is also described as an anticlockwise
+ (counter-clockwise) rotation of the coordinate frame.
+ It corresponds to as EPSG method code 1061 (in the geocentric domain) or
+ 1063 (in the geographic domain).
+
+ The result obtained with parameters specified in a given convention
+ can be obtained in the other convention by negating the rotational parameters
+ (``rx``, ``ry``, ``rz``)
+
+.. option:: +x=<value>
+
+ Translation of the x-axis given in meters.
+
+.. option:: +y=<value>
+
+ Translation of the y-axis given in meters.
+
+.. option:: +z=<value>
+
+ Translation of the z-axis given in meters.
+
+.. option:: +s=<value>
+
+ Scale factor given in ppm.
+
+.. option:: +rx=<value>
+
+ X-axis rotation given arc seconds.
+
+.. option:: +ry=<value>
+
+ Y-axis rotation given in arc seconds.
+
+.. option:: +rz=<value>
+
+ Z-axis rotation given in arc seconds.
+
+.. option:: +px=<value>
+
+ Coordinate along the x-axis of the reference point given in meters.
+
+.. option:: +py=<value>
+
+ Coordinate along the y-axis of the reference point given in meters.
+
+.. option:: +pz=<value>
+
+ Coordinate along the z-axis of the reference point given in meters.
+
+Mathematical description
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+In the *Position Vector* convention, we define :math:`R_x = radians \left( rx \right)`,
+:math:`R_z = radians \left( ry \right)` and :math:`R_z = radians \left( rz \right)`
+
+In the *Coordinate Frame* convention, :math:`R_x = - radians \left( rx \right)`,
+:math:`R_z = - radians \left( ry \right)` and :math:`R_z = - radians \left( rz \right)`
+
+.. math::
+ :label: 10param
+
+ \begin{align}
+ \begin{bmatrix}
+ X \\
+ Y \\
+ Z \\
+ \end{bmatrix}^{output} =
+ \begin{bmatrix}
+ T_x + P_x\\
+ T_y + P_y\\
+ T_z + P_z\\
+ \end{bmatrix} +
+ \left(1 + s \times 10^{-6}\right)
+ \begin{bmatrix}
+ 1 & -R_z & R_y \\
+ Rz & 1 & -R_x \\
+ -Ry & R_x & 1 \\
+ \end{bmatrix}
+ \begin{bmatrix}
+ X^{input} - P_x\\
+ Y^{input} - P_y\\
+ Z^{input} - P_z\\
+ \end{bmatrix}
+ \end{align}
diff --git a/docs/source/resource_files.rst b/docs/source/resource_files.rst
index 4f7721c0..bd21486c 100644
--- a/docs/source/resource_files.rst
+++ b/docs/source/resource_files.rst
@@ -103,13 +103,6 @@ Brazil
`Brazilian grids <http://www.ibge.gov.br/home/geociencias/geodesia/param_transf/default_param_transf.shtm>`__ for datums Corrego Alegre 1961, Corrego Alegre 1970-72, SAD69 and SAD69(96)
-Great Britain
-................................................................................
-
-`Great Britain's OSTN15_NTv2: OSGB 1936 => ETRS89 <https://www.ordnancesurvey.co.uk/docs/gps/OSTN15_NTv2.zip>`__
-
-`Great Britain's OSTN02_NTv2: OSGB 1936 => ETRS89 <http://www.ordnancesurvey.co.uk/business-and-government/help-and-support/navigation-technology/os-net/ostn02-ntv2-format.html>`__
-
Netherlands
................................................................................
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index 7acbbb36..757cf950 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -10,6 +10,10 @@
Implements 3(6)-, 4(8) and 7(14)-parameter Helmert transformations for
3D data.
+ Also incorporates Molodensky-Badekas variant of 7-parameter Helmert
+ transformation, where the rotation is not applied regarding the centre
+ of the spheroid, but given a reference point.
+
Primarily useful for implementation of datum shifts in transformation
pipelines.
@@ -17,7 +21,8 @@
Thomas Knudsen, thokn@sdfe.dk, 2016-05-24/06-05
Kristian Evers, kreve@sdfe.dk, 2017-05-01
-Last update: 2017-05-15
+Even Rouault, even.roault@spatialys.com
+Last update: 2018-10-26
************************************************************************
* Copyright (c) 2016, Thomas Knudsen / SDFE
@@ -52,6 +57,7 @@ Last update: 2017-05-15
#include "geocent.h"
PROJ_HEAD(helmert, "3(6)-, 4(8)- and 7(14)-parameter Helmert shift");
+PROJ_HEAD(molobadekas, "Molodensky-Badekas transformation");
static XYZ helmert_forward_3d (LPZ lpz, PJ *P);
static LPZ helmert_reverse_3d (XYZ xyz, PJ *P);
@@ -66,6 +72,7 @@ struct pj_opaque_helmert {
XYZ xyz;
XYZ xyz_0;
XYZ dxyz;
+ XYZ refp;
PJ_OPK opk;
PJ_OPK opk_0;
PJ_OPK dopk;
@@ -375,16 +382,16 @@ static XYZ helmert_forward_3d (LPZ lpz, PJ *P) {
scale = 1 + Q->scale * 1e-6;
- X = lpz.lam;
- Y = lpz.phi;
- Z = lpz.z;
+ X = lpz.lam - Q->refp.x;
+ Y = lpz.phi - Q->refp.y;
+ Z = lpz.z - Q->refp.z;
point.xyz.x = scale * ( R00 * X + R01 * Y + R02 * Z);
point.xyz.y = scale * ( R10 * X + R11 * Y + R12 * Z);
point.xyz.z = scale * ( R20 * X + R21 * Y + R22 * Z);
- point.xyz.x += Q->xyz.x;
+ point.xyz.x += Q->xyz.x; /* for Molodensky-Badekas, Q->xyz already incorporates the Q->refp offset */
point.xyz.y += Q->xyz.y;
point.xyz.z += Q->xyz.z;
@@ -421,9 +428,9 @@ static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) {
Z = (xyz.z - Q->xyz.z) / scale;
/* Inverse rotation through transpose multiplication */
- point.xyz.x = ( R00 * X + R10 * Y + R20 * Z);
- point.xyz.y = ( R01 * X + R11 * Y + R21 * Z);
- point.xyz.z = ( R02 * X + R12 * Y + R22 * Z);
+ point.xyz.x = ( R00 * X + R10 * Y + R20 * Z) + Q->refp.x;
+ point.xyz.y = ( R01 * X + R11 * Y + R21 * Z) + Q->refp.y;
+ point.xyz.z = ( R02 * X + R12 * Y + R22 * Z) + Q->refp.z;
return point.lpz;
}
@@ -467,30 +474,108 @@ static PJ_COORD helmert_reverse_4d (PJ_COORD point, PJ *P) {
/* Arcsecond to radians */
#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0)
-/***********************************************************************/
-PJ *TRANSFORMATION(helmert, 0) {
-/***********************************************************************/
+
+static PJ* init_helmert_six_parameters(PJ* P) {
struct pj_opaque_helmert *Q = pj_calloc (1, sizeof (struct pj_opaque_helmert));
if (0==Q)
return pj_default_destructor (P, ENOMEM);
P->opaque = (void *) Q;
- P->fwd4d = helmert_forward_4d;
- P->inv4d = helmert_reverse_4d;
- P->fwd3d = helmert_forward_3d;
- P->inv3d = helmert_reverse_3d;
- P->fwd = helmert_forward;
- P->inv = helmert_reverse;
-
/* In most cases, we work on 3D cartesian coordinates */
P->left = PJ_IO_UNITS_CARTESIAN;
P->right = PJ_IO_UNITS_CARTESIAN;
- /* But in the 2D case, the coordinates are projected */
+
+ /* Translations */
+ if (pj_param (P->ctx, P->params, "tx").i)
+ Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
+
+ if (pj_param (P->ctx, P->params, "ty").i)
+ Q->xyz_0.y = pj_param (P->ctx, P->params, "dy").f;
+
+ if (pj_param (P->ctx, P->params, "tz").i)
+ Q->xyz_0.z = pj_param (P->ctx, P->params, "dz").f;
+
+ /* Rotations */
+ if (pj_param (P->ctx, P->params, "trx").i)
+ Q->opk_0.o = pj_param (P->ctx, P->params, "drx").f * ARCSEC_TO_RAD;
+
+ if (pj_param (P->ctx, P->params, "try").i)
+ Q->opk_0.p = pj_param (P->ctx, P->params, "dry").f * ARCSEC_TO_RAD;
+
+ if (pj_param (P->ctx, P->params, "trz").i)
+ Q->opk_0.k = pj_param (P->ctx, P->params, "drz").f * ARCSEC_TO_RAD;
+
+ /* Use small angle approximations? */
+ if (pj_param (P->ctx, P->params, "bexact").i)
+ Q->exact = 1;
+
+ return P;
+}
+
+
+static PJ* read_convention(PJ* P) {
+
+ struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *)P->opaque;
+
+ /* In case there are rotational terms, we require an explicit convention
+ * to be provided. */
+ if (!Q->no_rotation) {
+ const char* convention = pj_param (P->ctx, P->params, "sconvention").s;
+ if( !convention ) {
+ proj_log_error (P, "helmert: missing 'convention' argument");
+ return pj_default_destructor (P, PJD_ERR_MISSING_ARGS);
+ }
+ if( strcmp(convention, "position_vector") == 0 ) {
+ Q->is_position_vector = 1;
+ }
+ else if( strcmp(convention, "coordinate_frame") == 0 ) {
+ Q->is_position_vector = 0;
+ }
+ else {
+ proj_log_error (P, "helmert: invalid value for 'convention' argument");
+ return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
+ }
+
+ /* historically towgs84 in PROJ has always been using position_vector
+ * convention. Accepting coordinate_frame would be confusing. */
+ if (pj_param_exists (P->params, "towgs84")) {
+ if( !Q->is_position_vector ) {
+ proj_log_error (P, "helmert: towgs84 should only be used with "
+ "convention=position_vector");
+ return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
+ }
+ }
+ }
+
+ return P;
+}
+
+
+/***********************************************************************/
+PJ *TRANSFORMATION(helmert, 0) {
+/***********************************************************************/
+
+ struct pj_opaque_helmert *Q;
+
+ if( !init_helmert_six_parameters(P) ) {
+ return 0;
+ }
+
+ /* In the 2D case, the coordinates are projected */
if (pj_param_exists (P->params, "theta")) {
P->left = PJ_IO_UNITS_PROJECTED;
P->right = PJ_IO_UNITS_PROJECTED;
}
+ P->fwd4d = helmert_forward_4d;
+ P->inv4d = helmert_reverse_4d;
+ P->fwd3d = helmert_forward_3d;
+ P->inv3d = helmert_reverse_3d;
+ P->fwd = helmert_forward;
+ P->inv = helmert_reverse;
+
+ Q = (struct pj_opaque_helmert *)P->opaque;
+
/* Detect obsolete transpose flag and error out if found */
if (pj_param (P->ctx, P->params, "ttranspose").i) {
proj_log_error (P, "helmert: 'transpose' argument is no longer valid. "
@@ -517,26 +602,6 @@ PJ *TRANSFORMATION(helmert, 0) {
Q->scale_0 = (P->datum_params[6] - 1) * 1e6;
}
- /* Translations */
- if (pj_param (P->ctx, P->params, "tx").i)
- Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
-
- if (pj_param (P->ctx, P->params, "ty").i)
- Q->xyz_0.y = pj_param (P->ctx, P->params, "dy").f;
-
- if (pj_param (P->ctx, P->params, "tz").i)
- Q->xyz_0.z = pj_param (P->ctx, P->params, "dz").f;
-
- /* Rotations */
- if (pj_param (P->ctx, P->params, "trx").i)
- Q->opk_0.o = pj_param (P->ctx, P->params, "drx").f * ARCSEC_TO_RAD;
-
- if (pj_param (P->ctx, P->params, "try").i)
- Q->opk_0.p = pj_param (P->ctx, P->params, "dry").f * ARCSEC_TO_RAD;
-
- if (pj_param (P->ctx, P->params, "trz").i)
- Q->opk_0.k = pj_param (P->ctx, P->params, "drz").f * ARCSEC_TO_RAD;
-
if (pj_param (P->ctx, P->params, "ttheta").i) {
Q->theta_0 = pj_param (P->ctx, P->params, "dtheta").f * ARCSEC_TO_RAD;
Q->fourparam = 1;
@@ -585,10 +650,6 @@ PJ *TRANSFORMATION(helmert, 0) {
if (pj_param(P->ctx, P->params, "tt_obs").i)
Q->t_obs = pj_param (P->ctx, P->params, "dt_obs").f;
- /* Use small angle approximations? */
- if (pj_param (P->ctx, P->params, "bexact").i)
- Q->exact = 1;
-
Q->xyz = Q->xyz_0;
Q->opk = Q->opk_0;
Q->scale = Q->scale_0;
@@ -599,34 +660,8 @@ PJ *TRANSFORMATION(helmert, 0) {
Q->no_rotation = 1;
}
- /* In case there are rotational terms, we require an explicit convention
- * to be provided. */
- if (!Q->no_rotation) {
- const char* convention = pj_param (P->ctx, P->params, "sconvention").s;
- if( !convention ) {
- proj_log_error (P, "helmert: missing 'convention' argument");
- return pj_default_destructor (P, PJD_ERR_MISSING_ARGS);
- }
- if( strcmp(convention, "position_vector") == 0 ) {
- Q->is_position_vector = 1;
- }
- else if( strcmp(convention, "coordinate_frame") == 0 ) {
- Q->is_position_vector = 0;
- }
- else {
- proj_log_error (P, "helmert: invalid value for 'convention' argument");
- return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
- }
-
- /* historically towgs84 in PROJ has always been using position_vector
- * convention. Accepting coordinate_frame would be confusing. */
- if (pj_param_exists (P->params, "towgs84")) {
- if( !Q->is_position_vector ) {
- proj_log_error (P, "helmert: towgs84 should only be used with "
- "convention=position_vector");
- return pj_default_destructor (P, PJD_ERR_INVALID_ARG);
- }
- }
+ if( !read_convention(P) ) {
+ return 0;
}
/* Let's help with debugging */
@@ -653,3 +688,66 @@ PJ *TRANSFORMATION(helmert, 0) {
return P;
}
+
+
+/***********************************************************************/
+PJ *TRANSFORMATION(molobadekas, 0) {
+/***********************************************************************/
+
+ struct pj_opaque_helmert *Q;
+
+ if( !init_helmert_six_parameters(P) ) {
+ return 0;
+ }
+
+ P->fwd3d = helmert_forward_3d;
+ P->inv3d = helmert_reverse_3d;
+
+ Q = (struct pj_opaque_helmert *)P->opaque;
+
+ /* Scale */
+ if (pj_param (P->ctx, P->params, "ts").i) {
+ Q->scale_0 = pj_param (P->ctx, P->params, "ds").f;
+ }
+
+ Q->opk = Q->opk_0;
+ Q->scale = Q->scale_0;
+
+ if( !read_convention(P) ) {
+ return 0;
+ }
+
+ /* Reference point */
+ if (pj_param (P->ctx, P->params, "tpx").i)
+ Q->refp.x = pj_param (P->ctx, P->params, "dpx").f;
+
+ if (pj_param (P->ctx, P->params, "tpy").i)
+ Q->refp.y = pj_param (P->ctx, P->params, "dpy").f;
+
+ if (pj_param (P->ctx, P->params, "tpz").i)
+ Q->refp.z = pj_param (P->ctx, P->params, "dpz").f;
+
+
+ /* Let's help with debugging */
+ if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
+ proj_log_debug(P, "Molodensky-Badekas parameters:");
+ proj_log_debug(P, "x= %8.5f y= %8.5f z= %8.5f", Q->xyz_0.x, Q->xyz_0.y, Q->xyz_0.z);
+ proj_log_debug(P, "rx= %8.5f ry= %8.5f rz= %8.5f",
+ Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD);
+ proj_log_debug(P, "s= %8.5f exact=%d%s", Q->scale, Q->exact,
+ Q->is_position_vector ? " convention=position_vector" :
+ " convention=coordinate_frame");
+ proj_log_debug(P, "px= %8.5f py= %8.5f pz= %8.5f", Q->refp.x, Q->refp.y, Q->refp.z);
+ }
+
+ /* as an optimization, we incorporate the refp in the translation terms */
+ Q->xyz_0.x += Q->refp.x;
+ Q->xyz_0.y += Q->refp.y;
+ Q->xyz_0.z += Q->refp.z;
+
+ Q->xyz = Q->xyz_0;
+
+ build_rot_matrix(P);
+
+ return P;
+}
diff --git a/src/pj_list.h b/src/pj_list.h
index 8c72dd1f..3592dcc4 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -94,6 +94,7 @@ PROJ_HEAD(mil_os, "Miller Oblated Stereographic")
PROJ_HEAD(mill, "Miller Cylindrical")
PROJ_HEAD(misrsom, "Space oblique for MISR")
PROJ_HEAD(moll, "Mollweide")
+PROJ_HEAD(molobadekas, "Molodensky-Badekas transform") /* implemented in PJ_helmert.c */
PROJ_HEAD(molodensky, "Molodensky transform")
PROJ_HEAD(murd1, "Murdoch I")
PROJ_HEAD(murd2, "Murdoch II")
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 6ac0c70c..13b77b0a 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -411,6 +411,30 @@ expect failure errno invalid_arg
-------------------------------------------------------------------------------
+Molodensky-Badekas from IOGP Guidance 7.2, Transformation from La Canoa to REGVEN
+between geographic 2D coordinate reference systems (EPSG Dataset transformation code 1771).
+Here just taking the Cartesian step of the transformation.
+-------------------------------------------------------------------------------
+operation proj=molobadekas convention=coordinate_frame
+ x=-270.933 y=115.599 z=-360.226 rx=-5.266 ry=-1.238 rz=2.381
+ s=-5.109 px=2464351.59 py=-5783466.61 pz=974809.81
+-------------------------------------------------------------------------------
+tolerance 1 cm
+roundtrip 1
+accept 2550408.96 -5749912.26 1054891.11
+expect 2550138.45 -5749799.87 1054530.82
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+Test error cases of molobadekas
+-------------------------------------------------------------------------------
+
+# Missing convention
+operation proj=molobadekas
+expect failure errno missing_arg
+
+
+-------------------------------------------------------------------------------
geocentric latitude
-------------------------------------------------------------------------------
operation proj=geoc ellps=GRS80