diff options
34 files changed, 842 insertions, 122 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 12a24b7b..02f49711 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,8 +19,6 @@ if(NOT CMAKE_VERSION VERSION_LESS 3.1) cmake_policy(SET CMP0054 NEW) endif() -add_definitions(-DPROJ_COMPILATION=1) - # Set C++ version # Make CMAKE_CXX_STANDARD available as cache option overridable by user set(CMAKE_CXX_STANDARD 11 @@ -34,19 +34,18 @@ RUN apt-get update; \ libsqlite3-0 \ curl unzip -COPY --from=builder /build/usr/bin/ /usr/bin/ -COPY --from=builder /build/usr/lib/ /usr/lib/ -COPY --from=builder /build/usr/include/ /usr/include/ -COPY --from=builder /build/usr/share/proj/ /usr/share/proj/ - +# Put this first as this is rarely changing RUN \ + mkdir -p /usr/share/proj; \ curl -LOs http://download.osgeo.org/proj/proj-datumgrid-1.8.zip && unzip -j -u -o proj-datumgrid-1.8.zip -d /usr/share/proj; \ curl -LOs http://download.osgeo.org/proj/proj-datumgrid-europe-1.2.zip && unzip -j -u -o proj-datumgrid-europe-1.2.zip -d /usr/share/proj; \ curl -LOs http://download.osgeo.org/proj/proj-datumgrid-oceania-1.0.zip && unzip -j -u -o proj-datumgrid-oceania-1.0.zip -d /usr/share/proj; \ curl -LOs http://download.osgeo.org/proj/proj-datumgrid-world-1.0.zip && unzip -j -u -o proj-datumgrid-world-1.0.zip -d /usr/share/proj; \ - curl -LOs http://download.osgeo.org/proj/proj-datumgrid-north-america-1.2.zip && unzip -j -u -o proj-datumgrid-north-america-1.2.zip -d /usr/share/proj; - - - + curl -LOs http://download.osgeo.org/proj/proj-datumgrid-north-america-1.2.zip && unzip -j -u -o proj-datumgrid-north-america-1.2.zip -d /usr/share/proj; \ + rm *.zip +COPY --from=builder /build/usr/share/proj/ /usr/share/proj/ +COPY --from=builder /build/usr/include/ /usr/include/ +COPY --from=builder /build/usr/bin/ /usr/bin/ +COPY --from=builder /build/usr/lib/ /usr/lib/ @@ -4,7 +4,7 @@ [](https://ci.appveyor.com/project/OSGeo/proj-4) [](https://coveralls.io/github/OSGeo/proj.4?branch=master) [](https://gitter.im/OSGeo/proj.4) -[](http://lists.maptools.org/mailman/listinfo/proj) +[](http://lists.osgeo.org/mailman/listinfo/proj) PROJ is a generic coordinate transformation software, that transforms coordinates from one coordinate reference system (CRS) to another. diff --git a/docs/source/apps/projinfo.rst b/docs/source/apps/projinfo.rst index c8219838..955b9a18 100644 --- a/docs/source/apps/projinfo.rst +++ b/docs/source/apps/projinfo.rst @@ -29,12 +29,17 @@ Synopsis | {object_definition} | (-s {srs_def} -t {srs_def}) | - where {object_definition} or {object_definition} is a PROJ string, a - WKT string, an object name, AUTHORITY:CODE - (where AUTHORITY is the name of a CRS authority and CODE the code of a CRS - found in the proj.db database) or a OGC URN (such as "urn:ogc:def:crs:EPSG::4326", - "urn:ogc:def:coordinateOperation:EPSG::1671", "urn:ogc:def:ellipsoid:EPSG::7001" - or "urn:ogc:def:datum:EPSG::6326") + where {object_definition} or {srs_def} is + + - a proj-string, + - a WKT string, + - an object code (like "EPSG:4326", "urn:ogc:def:crs:EPSG::4326", + "urn:ogc:def:coordinateOperation:EPSG::1671"), + - a OGC URN combining references for compound coordinate reference systems + (e.g "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717" or custom abbreviated + syntax "EPSG:2393+5717"), + - a OGC URN combining references for concatenated operations + (e.g. "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618") Description *********** diff --git a/docs/source/development/migration.rst b/docs/source/development/migration.rst index 99143f20..df622ecb 100644 --- a/docs/source/development/migration.rst +++ b/docs/source/development/migration.rst @@ -1,6 +1,151 @@ .. _API_migration: ================================================================================ +Version 4 to 6 API Migration +================================================================================ + +This is a transition guide for developers wanting to migrate their code to use +PROJ version 6. + +Code example +############################################################################### + +The difference between the old and new API is shown here with a few examples. Below +we implement the same program with the two different API's. The program reads +input longitude and latitude from the command line and convert them to +projected coordinates with the Mercator projection. + +We start by writing the program for PROJ 4: + +.. code-block:: C + + #include <proj_api.h> + + main(int argc, char **argv) { + projPJ pj_merc, pj_latlong; + double x, y; + + if (!(pj_longlat = pj_init_plus("+proj=longlat +ellps=clrk66")) ) + return 1; + if (!(pj_merc = pj_init_plus("+proj=merc +datum=clrk66 +lat_ts=33")) ) + return 1; + + while (scanf("%lf %lf", &x, &y) == 2) { + x *= DEG_TO_RAD; /* longitude */ + y *= DEG_TO_RAD; /* latitude */ + p = pj_transform(pj_longlat, pj_merc, 1, 1, &x, &y, NULL ); + printf("%.2f\t%.2f\n", x, y); + } + + pj_free(pj_longlat); + pj_free(pj_merc); + + return 0; + } + +The same program implemented using PROJ 6: + +.. code-block:: C + + #include <proj.h> + + main(int argc, char **argv) { + PJ *P; + PJ_COORD c; + + /* NOTE: the use of PROJ strings to describe CRS is strongly discouraged */ + /* in PROJ 6, as PROJ strings are a poor way of describing a CRS, and */ + /* more precise its geodetic datum. */ + /* Use of codes provided by authorities (such as "EPSG:4326", etc...) */ + /* or WKT strings will bring the full power of the "transformation */ + /* engine" used by PROJ to determine the best transformation(s) between */ + /* two CRS. */ + P = proj_create_crs_to_crs(PJ_DEFAULT_CTX, + "+proj=longlat +ellps=clrs66", + "+proj=merc +ellps=clrk66 +lat_ts=33", + NULL); + if (P==0) + return 1; + + { + /* For that particular use case, this is not needed. */ + /* proj_normalize_for_visualization() ensures that the coordinate */ + /* order expected and returned by proj_trans() will be longitude, */ + /* latitude for geographic CRS, and easting, northing for projected */ + /* CRS. If instead of using PROJ strings as above, "EPSG:XXXX" codes */ + /* had been used, this might had been necessary. */ + PJ* P_for_GIS = proj_normalize_for_visualization(C, P); + if( 0 == P_for_GIS ) { + proj_destroy(P); + return 1; + } + proj_destroy(P); + P = P_for_GIS; + } + + while (scanf("%lf %lf", &c.lp.lam, &c.lp.phi) == 2) { + /* No need to convert to radian */ + c = proj_trans(P, PJ_FWD, c); + printf("%.2f\t%.2f\n", c.xy.x, c.xy.y); + } + + proj_destroy(P); + } + + +Function mapping from old to new API +############################################################################### + ++---------------------------------------+-------------------------------------------------+ +| Old API functions | New API functions | ++=======================================+=================================================+ +| pj_fwd | :c:func:`proj_trans` | ++---------------------------------------+-------------------------------------------------+ +| pj_inv | :c:func:`proj_trans` | ++---------------------------------------+-------------------------------------------------+ +| pj_fwd3 | :c:func:`proj_trans` | ++---------------------------------------+-------------------------------------------------+ +| pj_inv3 | :c:func:`proj_trans` | ++---------------------------------------+-------------------------------------------------+ +| pj_transform | :c:func:`proj_create_crs_to_crs` + | +| | (:c:func:`proj_normalize_for_visualization` +) | +| | :c:func:`proj_trans`, | +| | :c:func:`proj_trans_array` or | +| | :c:func:`proj_trans_generic` | ++---------------------------------------+-------------------------------------------------+ +| pj_init | :c:func:`proj_create` / | +| | :c:func:`proj_create_crs_to_crs` | ++---------------------------------------+-------------------------------------------------+ +| pj_init | :c:func:`proj_create` / | +| | :c:func:`proj_create_crs_to_crs` | ++---------------------------------------+-------------------------------------------------+ +| pj_free | :c:func:`proj_destroy` | ++---------------------------------------+-------------------------------------------------+ +| pj_is_latlong | :c:func:`proj_get_type` | ++---------------------------------------+-------------------------------------------------+ +| pj_is_geocent | :c:func:`proj_get_type` | ++---------------------------------------+-------------------------------------------------+ +| pj_get_def | :c:func:`proj_pj_info` | ++---------------------------------------+-------------------------------------------------+ +| pj_latlong_from_proj | *No direct equivalent*, but can be accomplished | +| | by chaining :c:func:`proj_create`, | +| | :c:func:`proj_crs_get_horizontal_datum` and | +| | :c:func:`proj_create_geographic_crs_from_datum` | ++---------------------------------------+-------------------------------------------------+ +| pj_set_finder | :c:func:`proj_context_set_file_finder` | ++---------------------------------------+-------------------------------------------------+ +| pj_set_searchpath | :c:func:`proj_context_set_search_paths` | ++---------------------------------------+-------------------------------------------------+ +| pj_deallocate_grids | *No equivalent* | ++---------------------------------------+-------------------------------------------------+ +| pj_strerrno | *No equivalent* | ++---------------------------------------+-------------------------------------------------+ +| pj_get_errno_ref | :c:func:`proj_errno` | ++---------------------------------------+-------------------------------------------------+ +| pj_get_release | :c:func:`proj_info` | ++---------------------------------------+-------------------------------------------------+ + +================================================================================ Version 4 to 5 API Migration ================================================================================ @@ -66,7 +211,7 @@ Code example The difference between the old and new API is shown here with a few examples. Below we implement the same program with the two different API's. The program reads -input latitude and longitude from the command line and convert them to +input longitude and latitude from the command line and convert them to projected coordinates with the Mercator projection. We start by writing the program for PROJ v. 4: @@ -76,21 +221,24 @@ We start by writing the program for PROJ v. 4: #include <proj_api.h> main(int argc, char **argv) { - projPJ pj_merc, pj_latlong; + projPJ pj_merc, pj_longlat; double x, y; - if (!(pj_merc = pj_init_plus("+proj=merc +ellps=clrk66 +lat_ts=33")) ) + if (!(pj_longlat = pj_init_plus("+proj=longlat +ellps=clrk66")) ) return 1; - if (!(pj_latlong = pj_init_plus("+proj=latlong +ellps=clrk66")) ) + if (!(pj_merc = pj_init_plus("+proj=merc +ellps=clrk66 +lat_ts=33")) ) return 1; while (scanf("%lf %lf", &x, &y) == 2) { - x *= DEG_TO_RAD; - y *= DEG_TO_RAD; - p = pj_transform(pj_latlong, pj_merc, 1, 1, &x, &y, NULL ); + x *= DEG_TO_RAD; /* longitude */ + y *= DEG_TO_RAD; /* latitude */ + p = pj_transform(pj_longlat, pj_merc, 1, 1, &x, &y, NULL ); printf("%.2f\t%.2f\n", x, y); } + pj_free(pj_longlat); + pj_free(pj_merc); + return 0; } @@ -115,6 +263,7 @@ The same program implemented using PROJ v. 5: printf("%.2f\t%.2f\n", c.xy.x, c.xy.y); } + proj_destroy(P); } Looking at the two different programs, there's a few immediate @@ -155,27 +304,27 @@ Function mapping from old to new API +---------------------------------------+---------------------------------------+ | Old API functions | New API functions | +=======================================+=======================================+ -| pj_fwd | proj_trans | +| pj_fwd | :c:func:`proj_trans` | +---------------------------------------+---------------------------------------+ -| pj_inv | proj_trans | +| pj_inv | :c:func:`proj_trans` | +---------------------------------------+---------------------------------------+ -| pj_fwd3 | proj_trans | +| pj_fwd3 | :c:func:`proj_trans` | +---------------------------------------+---------------------------------------+ -| pj_inv3 | proj_trans | +| pj_inv3 | :c:func:`proj_trans` | +---------------------------------------+---------------------------------------+ | pj_transform | proj_trans_array or proj_trans_generic| +---------------------------------------+---------------------------------------+ -| pj_init | proj_create | +| pj_init | :c:func:`proj_create` | +---------------------------------------+---------------------------------------+ -| pj_init_plus | proj_create | +| pj_init_plus | :c:func:`proj_create` | +---------------------------------------+---------------------------------------+ -| pj_free | proj_destroy | +| pj_free | :c:func:`proj_destroy` | +---------------------------------------+---------------------------------------+ -| pj_is_latlong | proj_angular_output | +| pj_is_latlong | :c:func:`proj_angular_output` | +---------------------------------------+---------------------------------------+ -| pj_is_geocent | proj_angular_outout | +| pj_is_geocent | :c:func:`proj_angular_output` | +---------------------------------------+---------------------------------------+ -| pj_get_def | proj_pj_info | +| pj_get_def | :c:func:`proj_pj_info` | +---------------------------------------+---------------------------------------+ | pj_latlong_from_proj | *No equivalent* | +---------------------------------------+---------------------------------------+ @@ -187,7 +336,7 @@ Function mapping from old to new API +---------------------------------------+---------------------------------------+ | pj_strerrno | *No equivalent* | +---------------------------------------+---------------------------------------+ -| pj_get_errno_ref | proj_errno | +| pj_get_errno_ref | :c:func:`proj_errno` | +---------------------------------------+---------------------------------------+ -| pj_get_release | proj_info | +| pj_get_release | :c:func:`proj_info` | +---------------------------------------+---------------------------------------+ diff --git a/docs/source/development/quickstart.rst b/docs/source/development/quickstart.rst index b6c5d17c..d53d98fd 100644 --- a/docs/source/development/quickstart.rst +++ b/docs/source/development/quickstart.rst @@ -25,7 +25,7 @@ See the :doc:`reference for more info on data types <reference/datatypes>`. .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 43-45 + :lines: 43-46 :dedent: 4 For use in multi-threaded programs the :c:type:`PJ_CONTEXT` threading-context is used. @@ -35,13 +35,43 @@ this in detail. .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 49 + :lines: 50 :dedent: 4 Next we create the :c:type:`PJ` transformation object ``P`` with the function -:c:func:`proj_create`. :c:func:`proj_create` takes the threading context ``C`` -created above, and a proj-string that defines the desired transformation. -Here we transform from geodetic coordinate to UTM zone 32N. +:c:func:`proj_create_crs_to_crs`. :c:func:`proj_create_crs_to_crs` takes the threading context ``C`` +created above, a string that describes the source coordinate reference system (CRS), +a string that describes the target CRS and an optional description of the area of +use. +The strings for the source or target CRS may be PROJ strings (``+proj=longlat +datum=WGS84``), +CRS identified by their code (``EPSG:4326`` or ``urn:ogc:def:crs:EPSG::4326``) or +by a well-known text (WKT) string ( +:: + + GEOGCRS["WGS 84", + DATUM["World Geodetic System 1984", + ELLIPSOID["WGS 84",6378137,298.257223563, + LENGTHUNIT["metre",1]]], + PRIMEM["Greenwich",0, + ANGLEUNIT["degree",0.0174532925199433]], + CS[ellipsoidal,2], + AXIS["geodetic latitude (Lat)",north, + ORDER[1], + ANGLEUNIT["degree",0.0174532925199433]], + AXIS["geodetic longitude (Lon)",east, + ORDER[2], + ANGLEUNIT["degree",0.0174532925199433]], + USAGE[ + SCOPE["unknown"], + AREA["World"], + BBOX[-90,-180,90,180]], + ID["EPSG",4326]] + +). +The use of PROJ strings to describe a CRS is considered as legacy (one of the +main weakness of PROJ strings is their inability to describe a geodetic datum, +other than the few ones hardcoded in the ``+datum`` parameter). +Here we transform from geographic coordinates to UTM zone 32N. It is recommended to create one threading-context per thread used by the program. This ensures that all :c:type:`PJ` objects created in the same context will be sharing resources such as error-numbers and loaded grids. @@ -51,43 +81,72 @@ details. .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 51-53 + :lines: 52-60 :dedent: 4 -PROJ uses it's own data structures for handling coordinates. Here we use a +:c:func:`proj_create_crs_to_crs` creates a transformation object, which accepts +coordinates expressed in the units and axis order of the definition of the +source CRS, and return transformed coordinates in the units and axis order of +the definition of the target CRS. +For almost most geographic CRS, the units will be in most cases degrees (in +rare cases, such as EPSG:4807 / NTF (Paris), this can be grads). For geographic +CRS defined by the EPSG authority, the order of coordinates is latitude first, +longitude second. When using a PROJ string, on contrary the order will be +longitude first, latitude second. +For projected CRS, the units may vary (metre, us-foot, etc..). For projected +CRS defined by the EPSG authority, and with EAST / NORTH directions, the order +might be easting first, northing second, or the reverse. When using a PROJ string, +the order will be easting first, northing second, except if the ``+axis`` +parameter modifies it. + +If for the needs of your software, you want +a uniform axis order (and thus do not care about axis order mandated by the +authority defining the CRS), the :c:func:`proj_normalize_for_visualization` +function can be used to modify the PJ* object returned by +:c:func:`proj_create_crs_to_crs` so that it accepts as input and returns as +output coordinates using the traditional GIS order, that is longitude, latitude +(followed by elevation, time) for geographic CRS and easting, northing for most +projected CRS. + +.. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c + :language: c + :lines: 65-71 + :dedent: 4 + +PROJ uses its own data structures for handling coordinates. Here we use a :c:type:`PJ_COORD` which is easily assigned with the function :c:func:`proj_coord`. -Note that the input values are converted to radians with :c:func:`proj_torad`. -This is necessary since PROJ is using radians internally. See :doc:`transformations` -for further details. +When using +proj=longlat, the order of coordinates is longitude, latitude, +and values are expressed in degrees. If you used instead a EPSG geographic CRS, +like EPSG:4326 (WGS84), it would be latitude, longitude. .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 57 + :lines: 76 :dedent: 4 The coordinate defined above is transformed with :c:func:`proj_trans`. For this a :c:type:`PJ` object, a transformation direction (either forward or inverse) and the coordinate is needed. The transformed coordinate is returned in ``b``. -Here the forward (:c:type:`PJ_FWD`) transformation from geodetic to UTM is made. +Here the forward (:c:type:`PJ_FWD`) transformation from geographic to UTM is made. .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 60-61 + :lines: 79-80 :dedent: 4 -The inverse transformation (UTM to geodetic) is done similar to above, +The inverse transformation (UTM to geographic) is done similar to above, this time using :c:type:`PJ_INV` as the direction. .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 62-63 + :lines: 81-82 :dedent: 4 Before ending the program the allocated memory needs to be released again: .. literalinclude:: ../../../examples/pj_obs_api_mini_demo.c :language: c - :lines: 66-67 + :lines: 85-86 :dedent: 4 diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst index 00653ad9..64a4e8ca 100644 --- a/docs/source/development/reference/functions.rst +++ b/docs/source/development/reference/functions.rst @@ -29,9 +29,17 @@ paragraph for more details. .. c:function:: PJ* proj_create(PJ_CONTEXT *ctx, const char *definition) - Create a transformation object, or a CRS object, from a proj-string, - a WKT string, or object code (like "EPSG:4326", "urn:ogc:def:crs:EPSG::4326", - "urn:ogc:def:coordinateOperation:EPSG::1671"). + Create a transformation object, or a CRS object, from: + + - a proj-string, + - a WKT string, + - an object code (like "EPSG:4326", "urn:ogc:def:crs:EPSG::4326", + "urn:ogc:def:coordinateOperation:EPSG::1671"), + - a OGC URN combining references for compound coordinate reference systems + (e.g "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717" or custom abbreviated + syntax "EPSG:2393+5717"), + - a OGC URN combining references for concatenated operations + (e.g. "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618") Example call: @@ -110,7 +118,8 @@ paragraph for more details. - the name of a CRS as found in the PROJ database, e.g "WGS84", "NAD27", etc. - - more generally any string accepted by :c:func:`proj_create` + - more generally any string accepted by :c:func:`proj_create` representing + a CRS An "area of use" can be specified in area. When it is supplied, the more accurate transformation between two given systems can be chosen. @@ -144,6 +153,25 @@ paragraph for more details. :type `area`: PJ_AREA :returns: :c:type:`PJ*` +.. c:function:: PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ* obj) + + .. versionadded:: 6.1.0 + + Returns a PJ* object whose axis order is the one expected for + visualization purposes. + + The input object must be a coordinate operation, that has been created with + proj_create_crs_to_crs(). + If the axis order of its source or target CRS is northing,easting, then an + axis swap operation will be inserted. + + The returned :c:type:`PJ`-pointer should be deallocated with :c:func:`proj_destroy`. + + :param PJ_CONTEXT* ctx: Threading context. + :param `obj`: Object of type CoordinateOperation + :returns: :c:type:`PJ*` + + .. c:function:: PJ* proj_destroy(PJ *P) Deallocate a :c:type:`PJ` transformation object. diff --git a/docs/source/resource_files.rst b/docs/source/resource_files.rst index 53843571..d6f64175 100644 --- a/docs/source/resource_files.rst +++ b/docs/source/resource_files.rst @@ -106,7 +106,7 @@ Brazil Netherlands ................................................................................ -`Dutch grid <https://www.kadaster.nl/web/Themas/Registraties/Rijksdriehoeksmeting/Transformatie-van-coordinaten.htm>`__ (Registration required before download) +`Dutch grid <https://zakelijk.kadaster.nl/transformatie-van-coordinaten>`__ (Registration required before download) Portugal ................................................................................ diff --git a/examples/pj_obs_api_mini_demo.c b/examples/pj_obs_api_mini_demo.c index 94520490..3df94e2d 100644 --- a/examples/pj_obs_api_mini_demo.c +++ b/examples/pj_obs_api_mini_demo.c @@ -42,23 +42,42 @@ int main (void) { PJ_CONTEXT *C; PJ *P; + PJ* P_for_GIS; PJ_COORD a, b; /* or you may set C=PJ_DEFAULT_CTX if you are sure you will */ /* use PJ objects from only one thread */ C = proj_context_create(); - P = proj_create (C, "+proj=utm +zone=32 +ellps=GRS80"); - if (0==P) - return puts ("Oops"), 0; + P = proj_create_crs_to_crs (C, + "EPSG:4326", + "+proj=utm +zone=32 +datum=WGS84", /* or EPSG:32632 */ + NULL); + + if (0==P) { + fprintf(stderr, "Oops\n"); + return 1; + } + + /* This will ensure that the order of coordinates for the input CRS */ + /* will be longitude, latitude, whereas EPSG:4326 mandates latitude, */ + /* longitude */ + P_for_GIS = proj_normalize_for_visualization(C, P); + if( 0 == P_for_GIS ) { + fprintf(stderr, "Oops\n"); + return 1; + } + proj_destroy(P); + P = P_for_GIS; /* a coordinate union representing Copenhagen: 55d N, 12d E */ - /* note: PROJ.4 works in radians, hence the proj_torad() calls */ - a = proj_coord (proj_torad(12), proj_torad(55), 0, 0); + /* Given that we have used proj_normalize_for_visualization(), the order of + /* coordinates is longitude, latitude, and values are expressed in degrees. */ + a = proj_coord (12, 55, 0, 0); /* transform to UTM zone 32, then back to geographical */ b = proj_trans (P, PJ_FWD, a); - printf ("easting: %g, northing: %g\n", b.enu.e, b.enu.n); + printf ("easting: %.3f, northing: %.3f\n", b.enu.e, b.enu.n); b = proj_trans (P, PJ_INV, b); printf ("longitude: %g, latitude: %g\n", b.lp.lam, b.lp.phi); diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 5386ca6c..a603533e 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -993,16 +993,16 @@ class PROJ_GCC_DLL AuthorityFactory { PROJ_INTERNAL std::list<crs::CompoundCRSNNPtr> createCompoundCRSFromExisting(const crs::CompoundCRSNNPtr &crs) const; + + PROJ_INTERNAL crs::CRSNNPtr + createCoordinateReferenceSystem(const std::string &code, + bool allowCompound) const; //! @endcond protected: PROJ_INTERNAL AuthorityFactory(const DatabaseContextNNPtr &context, const std::string &authorityName); - PROJ_INTERNAL crs::CRSNNPtr - createCoordinateReferenceSystem(const std::string &code, - bool allowCompound) const; - PROJ_INTERNAL crs::GeodeticCRSNNPtr createGeodeticCRS(const std::string &code, bool geographicOnly) const; diff --git a/src/Makefile.am b/src/Makefile.am index 9858d78f..aa1787c1 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -8,7 +8,7 @@ check_PROGRAMS = geodtest AM_CPPFLAGS = -DPROJ_LIB=\"$(pkgdatadir)\" \ -DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ -I$(top_srcdir)/include @SQLITE3_CFLAGS@ -AM_CXXFLAGS = @CXX_WFLAGS@ @FLTO_FLAG@ -DPROJ_COMPILATION +AM_CXXFLAGS = @CXX_WFLAGS@ @FLTO_FLAG@ include_HEADERS = proj.h proj_experimental.h proj_constants.h proj_api.h geodesic.h \ org_proj4_PJ.h proj_symbol_rename.h @@ -16,11 +16,11 @@ include_HEADERS = proj.h proj_experimental.h proj_constants.h proj_api.h geodesi EXTRA_DIST = bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \ bin_geod.cmake bin_proj.cmake bin_projinfo.cmake \ lib_proj.cmake CMakeLists.txt bin_geodtest.cmake tests/geodtest.cpp \ - wkt1_grammar.y wkt2_grammar.y apps/emess.h + wkt1_grammar.y wkt2_grammar.y apps/emess.h apps/utils.h -proj_SOURCES = apps/proj.cpp apps/emess.cpp +proj_SOURCES = apps/proj.cpp apps/emess.cpp apps/utils.cpp projinfo_SOURCES = apps/projinfo.cpp -cs2cs_SOURCES = apps/cs2cs.cpp apps/emess.cpp +cs2cs_SOURCES = apps/cs2cs.cpp apps/emess.cpp apps/utils.cpp cct_SOURCES = apps/cct.cpp apps/proj_strtod.cpp apps/proj_strtod.h apps/optargpm.h geod_SOURCES = apps/geod.cpp apps/geod_set.cpp apps/geod_interface.cpp apps/geod_interface.h apps/emess.cpp diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp index 877a68ff..40b0d584 100644 --- a/src/apps/cs2cs.cpp +++ b/src/apps/cs2cs.cpp @@ -45,6 +45,7 @@ #include "proj.h" #include "proj_internal.h" #include "emess.h" +#include "utils.h" // clang-format on #define MAX_LINE 1000 @@ -522,6 +523,13 @@ int main(int argc, char **argv) { if (eargc == 0) /* if no specific files force sysin */ eargv[eargc++] = const_cast<char *>("-"); + if( oform ) { + if( !validate_form_string_for_numbers(oform) ) { + emess(3, "invalid format string"); + exit(0); + } + } + /* * If the user has requested inverse, then just reverse the * coordinate systems. diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp index 2af49c34..888d723f 100644 --- a/src/apps/proj.cpp +++ b/src/apps/proj.cpp @@ -7,6 +7,7 @@ #include <string.h> #include <math.h> #include "emess.h" +#include "utils.h" #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__WIN32__) # include <fcntl.h> @@ -461,6 +462,13 @@ int main(int argc, char **argv) { if (eargc == 0) /* if no specific files force sysin */ eargv[eargc++] = const_cast<char*>("-"); + if( oform ) { + if( !validate_form_string_for_numbers(oform) ) { + emess(3, "invalid format string"); + exit(0); + } + } + /* done with parameter and control input */ if (inverse && postscale) { prescale = 1; @@ -487,7 +495,6 @@ int main(int argc, char **argv) { proj.inv = pj_inv; } else proj.fwd = pj_fwd; - /* set input formatting control */ if (mon) { pj_pr_list(Proj); diff --git a/src/apps/utils.cpp b/src/apps/utils.cpp new file mode 100644 index 00000000..7dc809c9 --- /dev/null +++ b/src/apps/utils.cpp @@ -0,0 +1,58 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Utilities for command line arguments + * Author: Even Rouault <even dot rouault at spatialys dot com> + * + ****************************************************************************** + * Copyright (c) 2019, Even Rouault <even dot rouault at spatialys dot com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +#include "utils.h" + +#include <string.h> + +bool validate_form_string_for_numbers(const char* formatString) { + /* Only accepts '%[+]?[number]?[.]?[number]?[e|E|f|F|g|G]' */ + bool valid = true; + if( formatString[0] != '%' ) + valid = false; + else { + auto oformLen = strlen(formatString); + for( int i = 1; i < static_cast<int>(oformLen) - 1; i++ ) { + if( !(formatString[i] == '.' || + formatString[i] == '+' || + (formatString[i] >= '0' && formatString[i] <= '9')) ) { + valid = false; + break; + } + } + if( valid ) { + valid = formatString[oformLen-1] == 'e' || + formatString[oformLen-1] == 'E' || + formatString[oformLen-1] == 'f' || + formatString[oformLen-1] == 'F' || + formatString[oformLen-1] == 'g' || + formatString[oformLen-1] == 'G'; + } + } + return valid; +} diff --git a/src/apps/utils.h b/src/apps/utils.h new file mode 100644 index 00000000..99c14091 --- /dev/null +++ b/src/apps/utils.h @@ -0,0 +1,29 @@ +/****************************************************************************** + * + * Project: PROJ + * Purpose: Utilities for command line arguments + * Author: Even Rouault <even dot rouault at spatialys dot com> + * + ****************************************************************************** + * Copyright (c) 2019, Even Rouault <even dot rouault at spatialys dot com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************/ + +bool validate_form_string_for_numbers(const char* formatString); diff --git a/src/bin_cs2cs.cmake b/src/bin_cs2cs.cmake index 7ee26673..d2bb3b97 100644 --- a/src/bin_cs2cs.cmake +++ b/src/bin_cs2cs.cmake @@ -1,6 +1,7 @@ set(CS2CS_SRC apps/cs2cs.cpp apps/emess.cpp + apps/utils.cpp ) source_group("Source Files\\Bin" FILES ${CS2CS_SRC}) diff --git a/src/bin_proj.cmake b/src/bin_proj.cmake index b9ae03e5..ce282fc6 100644 --- a/src/bin_proj.cmake +++ b/src/bin_proj.cmake @@ -1,6 +1,7 @@ set(PROJ_SRC apps/proj.cpp apps/emess.cpp + apps/utils.cpp ) source_group("Source Files\\Bin" FILES ${PROJ_SRC}) diff --git a/src/conversions/cart.cpp b/src/conversions/cart.cpp index e6942d65..c1f6f09d 100644 --- a/src/conversions/cart.cpp +++ b/src/conversions/cart.cpp @@ -162,6 +162,12 @@ static PJ_LPZ geodetic (PJ_XYZ cart, PJ *P) { c = cos(theta); s = sin(theta); lpz.phi = atan2 (cart.z + P->e2s*P->b*s*s*s, p - P->es*P->a*c*c*c); + if( fabs(lpz.phi) > M_HALFPI ) { + // this happen on non-sphere ellipsoid when x,y,z is very close to 0 + // there is no single solution to the cart->geodetic conversion in + // that case, so arbitrarily pickup phi = 0. + lpz.phi = 0; + } lpz.lam = atan2 (cart.y, cart.x); N = normal_radius_of_curvature (P->a, P->es, lpz.phi); diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 8fea51e9..6a2a1ae0 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -5964,7 +5964,8 @@ int proj_coordoperation_get_param( * @param out_values Pointer to an array of value_count double values. * @param value_count Size of out_values array. The suggested size is 7 to get * translation, rotation and scale difference parameters. Rotation and scale - * difference terms might be zero if the transformation only includes translation + * difference terms might be zero if the transformation only includes + * translation * parameters. In that case, value_count could be set to 3. * @param emit_error_if_incompatible Boolean to inicate if an error must be * logged if coordoperation is not compatible with a WKT1 TOWGS84 diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 7eab849c..9c393aba 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -7484,6 +7484,8 @@ Transformation::Private::registerInv(util::BaseObjectNNPtr thisIn, TransformationNNPtr invTransform) { invTransform->d->forwardOperation_ = util::nn_dynamic_pointer_cast<Transformation>(thisIn); + invTransform->setHasBallparkTransformation( + invTransform->d->forwardOperation_->hasBallparkTransformation()); return invTransform; } //! @endcond @@ -9490,6 +9492,7 @@ CoordinateOperationNNPtr ConcatenatedOperation::inverse() const { auto op = create(properties, inversedOperations, coordinateOperationAccuracies()); op->d->computedName_ = d->computedName_; + op->setHasBallparkTransformation(hasBallparkTransformation()); return op; } @@ -12584,11 +12587,11 @@ getResolvedCRS(const crs::CRSNNPtr &crs, const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), *ids.front()->codeSpace()); try { - crs::CRSNNPtr resolvedCrs( + auto resolvedCrs( tmpAuthFactory->createProjectedCRS(ids.front()->code())); - if (resolvedCrs->_isEquivalentTo( + if (resolvedCrs->isEquivalentTo( crs.get(), util::IComparable::Criterion::EQUIVALENT)) { - return resolvedCrs; + return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs); } } catch (const std::exception &) { } @@ -12616,12 +12619,13 @@ getResolvedCRS(const crs::CRSNNPtr &crs, const auto tmpAuthFactory = io::AuthorityFactory::create( authFactory->databaseContext(), *ids.front()->codeSpace()); try { - crs::CRSNNPtr resolvedCrs( + auto resolvedCrs( tmpAuthFactory->createCompoundCRS(ids.front()->code())); - if (resolvedCrs->_isEquivalentTo( + if (resolvedCrs->isEquivalentTo( crs.get(), util::IComparable::Criterion::EQUIVALENT)) { - return resolvedCrs; + return util::nn_static_pointer_cast<crs::CRS>( + resolvedCrs); } } catch (const std::exception &) { } @@ -12665,6 +12669,17 @@ CoordinateOperationFactory::createOperations( auto l_resolvedTargetCRS = getResolvedCRS(l_targetCRS, context); Private::Context contextPrivate(l_resolvedSourceCRS, l_resolvedTargetCRS, context); + + if (context->getSourceAndTargetCRSExtentUse() == + CoordinateOperationContext::SourceTargetCRSExtentUse::INTERSECTION) { + auto sourceCRSExtent(getExtent(l_resolvedSourceCRS)); + auto targetCRSExtent(getExtent(l_resolvedTargetCRS)); + if (sourceCRSExtent && targetCRSExtent && + !sourceCRSExtent->intersects(NN_NO_CHECK(targetCRSExtent))) { + return std::vector<CoordinateOperationNNPtr>(); + } + } + return filterAndSort(Private::createOperations(l_resolvedSourceCRS, l_resolvedTargetCRS, contextPrivate), @@ -12703,6 +12718,8 @@ void InverseCoordinateOperation::setPropertiesFromForward() { if (forwardOperation_->sourceCRS() && forwardOperation_->targetCRS()) { setCRSs(forwardOperation_.get(), true); } + setHasBallparkTransformation( + forwardOperation_->hasBallparkTransformation()); } // --------------------------------------------------------------------------- @@ -12785,7 +12802,7 @@ PROJBasedOperationNNPtr PROJBasedOperation::create( auto method = OperationMethod::create( util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, - "PROJ-based operation method (approximate) : " + + "PROJ-based operation method (approximate): " + projString), std::vector<GeneralOperationParameterNNPtr>{}); auto op = PROJBasedOperation::nn_make_shared<PROJBasedOperation>(method); diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index f24b3457..9f7467a2 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -2484,6 +2484,8 @@ crs::CRSNNPtr AuthorityFactory::createCoordinateReferenceSystem( return createCoordinateReferenceSystem(code, true); } +//! @cond Doxygen_Suppress + crs::CRSNNPtr AuthorityFactory::createCoordinateReferenceSystem(const std::string &code, bool allowCompound) const { @@ -2513,6 +2515,9 @@ AuthorityFactory::createCoordinateReferenceSystem(const std::string &code, } throw FactoryException("unhandled CRS type: " + type); } + +//! @endcond + // --------------------------------------------------------------------------- //! @cond Doxygen_Suppress diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 360d55a2..0a32bb7c 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -1342,7 +1342,7 @@ struct WKTParser::Private { CRSPtr buildCRS(const WKTNodeNNPtr &node); - CoordinateOperationNNPtr buildCoordinateOperation(const WKTNodeNNPtr &node); + TransformationNNPtr buildCoordinateOperation(const WKTNodeNNPtr &node); ConcatenatedOperationNNPtr buildConcatenatedOperation(const WKTNodeNNPtr &node); @@ -2940,7 +2940,7 @@ WKTParser::Private::buildConversion(const WKTNodeNNPtr &node, // --------------------------------------------------------------------------- -CoordinateOperationNNPtr +TransformationNNPtr WKTParser::Private::buildCoordinateOperation(const WKTNodeNNPtr &node) { const auto *nodeP = node->GP(); auto &methodNode = nodeP->lookForChild(WKTConstants::METHOD); @@ -2991,11 +2991,10 @@ WKTParser::Private::buildCoordinateOperation(const WKTNodeNNPtr &node) { stripQuotes(accuracyNode->GP()->children()[0]))); } - return util::nn_static_pointer_cast<CoordinateOperation>( - Transformation::create(buildProperties(node), NN_NO_CHECK(sourceCRS), - NN_NO_CHECK(targetCRS), interpolationCRS, - buildProperties(methodNode), parameters, values, - accuracies)); + return Transformation::create(buildProperties(node), NN_NO_CHECK(sourceCRS), + NN_NO_CHECK(targetCRS), interpolationCRS, + buildProperties(methodNode), parameters, + values, accuracies); } // --------------------------------------------------------------------------- @@ -4307,16 +4306,31 @@ BaseObjectNNPtr WKTParser::Private::build(const WKTNodeNNPtr &node) { } if (ci_equal(name, WKTConstants::COORDINATEOPERATION)) { - return util::nn_static_pointer_cast<BaseObject>( - buildCoordinateOperation(node)); + auto transf = buildCoordinateOperation(node); + + const char *prefixes[] = { + "PROJ-based operation method: ", + "PROJ-based operation method (approximate): "}; + for (const char *prefix : prefixes) { + if (starts_with(transf->method()->nameStr(), prefix)) { + auto projString = + transf->method()->nameStr().substr(strlen(prefix)); + return util::nn_static_pointer_cast<BaseObject>( + PROJBasedOperation::create( + PropertyMap(), projString, transf->sourceCRS(), + transf->targetCRS(), + transf->coordinateOperationAccuracies())); + } + } + + return util::nn_static_pointer_cast<BaseObject>(transf); } if (ci_equal(name, WKTConstants::CONVERSION)) { auto conv = buildConversion(node, UnitOfMeasure::METRE, UnitOfMeasure::DEGREE); - if (conv->nameStr() == "PROJ-based coordinate operation" && - starts_with(conv->method()->nameStr(), + if (starts_with(conv->method()->nameStr(), "PROJ-based operation method: ")) { auto projString = conv->method()->nameStr().substr( strlen("PROJ-based operation method: ")); @@ -4397,15 +4411,97 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text, const auto authorities = dbContextNNPtr->getAuthorities(); for (const auto &authCandidate : authorities) { if (ci_equal(authCandidate, authName)) { - return AuthorityFactory::create(dbContextNNPtr, - authCandidate) - ->createCoordinateReferenceSystem(code); + factory = + AuthorityFactory::create(dbContextNNPtr, authCandidate); + try { + return factory->createCoordinateReferenceSystem(code); + } catch (...) { + // EPSG:4326+3855 + auto tokensCode = split(code, '+'); + if (tokensCode.size() == 2) { + auto crs1(factory->createCoordinateReferenceSystem( + tokensCode[0], false)); + auto crs2(factory->createCoordinateReferenceSystem( + tokensCode[1], false)); + return CompoundCRS::create( + util::PropertyMap().set( + IdentifiedObject::NAME_KEY, + crs1->nameStr() + " + " + crs2->nameStr()), + {crs1, crs2}); + } + throw; + } } } throw; } } + // OGC 07-092r2: para 7.5.2 + // URN combined references for compound coordinate reference systems + if (starts_with(text, "urn:ogc:def:crs,")) { + if (!dbContext) { + throw ParsingException("no database context specified"); + } + auto tokensComma = split(text, ','); + std::vector<CRSNNPtr> components; + std::string name; + for (size_t i = 1; i < tokensComma.size(); i++) { + tokens = split(tokensComma[i], ':'); + if (tokens.size() != 4) { + throw ParsingException( + concat("invalid crs component: ", tokensComma[i])); + } + const auto &type = tokens[0]; + auto factory = + AuthorityFactory::create(NN_NO_CHECK(dbContext), tokens[1]); + const auto &code = tokens[3]; + if (type == "crs") { + auto crs(factory->createCoordinateReferenceSystem(code, false)); + components.emplace_back(crs); + if (!name.empty()) { + name += " + "; + } + name += crs->nameStr(); + } else { + throw ParsingException( + concat("unexpected object type: ", type)); + } + } + return CompoundCRS::create( + util::PropertyMap().set(IdentifiedObject::NAME_KEY, name), + components); + } + + // OGC 07-092r2: para 7.5.3 + // 7.5.3 URN combined references for concatenated operations + if (starts_with(text, "urn:ogc:def:coordinateOperation,")) { + if (!dbContext) { + throw ParsingException("no database context specified"); + } + auto tokensComma = split(text, ','); + std::vector<CoordinateOperationNNPtr> components; + for (size_t i = 1; i < tokensComma.size(); i++) { + tokens = split(tokensComma[i], ':'); + if (tokens.size() != 4) { + throw ParsingException(concat( + "invalid coordinateOperation component: ", tokensComma[i])); + } + const auto &type = tokens[0]; + auto factory = + AuthorityFactory::create(NN_NO_CHECK(dbContext), tokens[1]); + const auto &code = tokens[3]; + if (type == "coordinateOperation") { + auto op(factory->createCoordinateOperation(code, false)); + components.emplace_back(op); + } else { + throw ParsingException( + concat("unexpected object type: ", type)); + } + } + return ConcatenatedOperation::createComputeMetadata(components, true); + } + // urn:ogc:def:crs:EPSG::4326 if (tokens.size() == 7) { if (!dbContext) { @@ -4510,6 +4606,13 @@ static BaseObjectNNPtr createFromUserInput(const std::string &text, * "urn:ogc:def:coordinateOperation:EPSG::1671", * "urn:ogc:def:ellipsoid:EPSG::7001" * or "urn:ogc:def:datum:EPSG::6326"</li> + * <li> OGC URN combining references for compound coordinate reference systems + * e.g. "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717" + * We also accept a custom abbreviated syntax EPSG:2393+5717 + * </li> + * <li> OGC URN combining references for concatenated operations + * e.g. + * "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618"</li> * <li>an Object name. e.g "WGS 84", "WGS 84 / UTM zone 31N". In that case as * uniqueness is not guaranteed, the function may apply heuristics to * determine the appropriate best match.</li> @@ -4546,6 +4649,13 @@ BaseObjectNNPtr createFromUserInput(const std::string &text, * "urn:ogc:def:coordinateOperation:EPSG::1671", * "urn:ogc:def:ellipsoid:EPSG::7001" * or "urn:ogc:def:datum:EPSG::6326"</li> + * <li> OGC URN combining references for compound coordinate reference systems + * e.g. "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717" + * We also accept a custom abbreviated syntax EPSG:2393+5717 + * </li> + * <li> OGC URN combining references for concatenated operations + * e.g. + * "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618"</li> * <li>an Object name. e.g "WGS 84", "WGS 84 / UTM zone 31N". In that case as * uniqueness is not guaranteed, the function may apply heuristics to * determine the appropriate best match.</li> diff --git a/src/projections/airy.cpp b/src/projections/airy.cpp index f7068061..ba6a40ff 100644 --- a/src/projections/airy.cpp +++ b/src/projections/airy.cpp @@ -79,6 +79,10 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ } if (fabs(s = 1. - cosz) > EPS) { t = 0.5 * (1. + cosz); + if(t == 0) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } Krho = -log(t)/s - Q->Cb / t; } else Krho = 0.5 - Q->Cb; diff --git a/src/projections/aitoff.cpp b/src/projections/aitoff.cpp index 127841ff..23554605 100644 --- a/src/projections/aitoff.cpp +++ b/src/projections/aitoff.cpp @@ -116,15 +116,20 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ sl = sin(lp.lam * 0.5); cl = cos(lp.lam * 0.5); sp = sin(lp.phi); cp = cos(lp.phi); D = cp * cl; - C = 1. - D * D; - D = acos(D) / pow(C, 1.5); - f1 = 2. * D * C * cp * sl; - f2 = D * C * sp; - f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl); - f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp; - f2p = sp * sp * cl / C + D * sl * sl * cp; - f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl); - if (Q->mode == WINKEL_TRIPEL) { + C = 1. - D * D; + const double denom = pow(C, 1.5); + if( denom == 0 ) { + proj_errno_set(P, PJD_ERR_NON_CONVERGENT); + return lp; + } + D = acos(D) / denom; + f1 = 2. * D * C * cp * sl; + f2 = D * C * sp; + f1p = 2.* (sl * cl * sp * cp / C - D * sp * sl); + f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp; + f2p = sp * sp * cl / C + D * sl * sl * cp; + f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl); + if (Q->mode == WINKEL_TRIPEL) { f1 = 0.5 * (f1 + lp.lam * Q->cosphi1); f2 = 0.5 * (f2 + lp.phi); f1p *= 0.5; @@ -156,10 +161,10 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ } while (((fabs(xy.x-x) > EPSILON) || (fabs(xy.y-y) > EPSILON)) && (round++ < MAXROUND)); if (iter == MAXITER && round == MAXROUND) - { - pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT ); - /* fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl); */ - } + { + pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT ); + /* fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl); */ + } return lp; } diff --git a/src/projections/isea.cpp b/src/projections/isea.cpp index 28510cb0..e8720b27 100644 --- a/src/projections/isea.cpp +++ b/src/projections/isea.cpp @@ -898,6 +898,10 @@ static int isea_hex(struct isea_dgg *g, int tri, quad = isea_ptdi(g, tri, pt, &v); + if( v.x < (INT_MIN >> 4) || v.x > (INT_MAX >> 4) ) + { + throw "Invalid shift"; + } hex->x = ((int)v.x << 4) + quad; hex->y = v.y; diff --git a/src/projections/mod_ster.cpp b/src/projections/mod_ster.cpp index b26ea289..50f66839 100644 --- a/src/projections/mod_ster.cpp +++ b/src/projections/mod_ster.cpp @@ -35,7 +35,12 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ pow((1. - esphi) / (1. + esphi), P->e * .5)) - M_HALFPI; schi = sin(chi); cchi = cos(chi); - s = 2. / (1. + Q->schio * schi + Q->cchio * cchi * coslon); + const double denom = 1. + Q->schio * schi + Q->cchio * cchi * coslon; + if( denom == 0 ) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return xy; + } + s = 2. / denom; p.r = s * cchi * sinlon; p.i = s * (Q->cchio * schi - Q->schio * cchi * coslon); p = pj_zpoly1(p, Q->zcoeff, Q->n); diff --git a/src/projections/nsper.cpp b/src/projections/nsper.cpp index a0bb5686..37938924 100644 --- a/src/projections/nsper.cpp +++ b/src/projections/nsper.cpp @@ -96,7 +96,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ PJ_LP lp = {0.0,0.0}; struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); - double rh, cosz, sinz; + double rh; if (Q->tilt) { double bm, bq, yt; @@ -108,16 +108,18 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ xy.y = bq * Q->cg - bm * Q->sg; } rh = hypot(xy.x, xy.y); - if ((sinz = 1. - rh * rh * Q->pfact) < 0.) { - proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); - return lp; - } - sinz = (Q->p - sqrt(sinz)) / (Q->pn1 / rh + rh / Q->pn1); - cosz = sqrt(1. - sinz * sinz); if (fabs(rh) <= EPS10) { lp.lam = 0.; lp.phi = P->phi0; } else { + double cosz, sinz; + sinz = 1. - rh * rh * Q->pfact; + if (sinz < 0.) { + proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION); + return lp; + } + sinz = (Q->p - sqrt(sinz)) / (Q->pn1 / rh + rh / Q->pn1); + cosz = sqrt(1. - sinz * sinz); switch (Q->mode) { case OBLIQ: lp.phi = asin(cosz * Q->sinph0 + xy.y * sinz * Q->cosph0 / rh); diff --git a/src/projections/omerc.cpp b/src/projections/omerc.cpp index e07c209e..0de3aa7d 100644 --- a/src/projections/omerc.cpp +++ b/src/projections/omerc.cpp @@ -154,6 +154,8 @@ PJ *PROJECTION(omerc) { phi1 = pj_param(P->ctx, P->params, "rlat_1").f; lam2 = pj_param(P->ctx, P->params, "rlon_2").f; phi2 = pj_param(P->ctx, P->params, "rlat_2").f; + if (fabs(phi1) > M_HALFPI || fabs(phi2) > M_HALFPI) + return pj_default_destructor(P, PJD_ERR_LAT_LARGER_THAN_90); if (fabs(phi1 - phi2) <= TOL || (con = fabs(phi1)) <= TOL || fabs(con - M_HALFPI) <= TOL || diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 631767a3..08b8c9b2 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -324,7 +324,7 @@ Airy =============================================================================== ------------------------------------------------------------------------------- -operation +proj=airy +a=6400000 +lat_1=0 +lat_2=2 +operation +proj=airy +a=6400000 ------------------------------------------------------------------------------- tolerance 0.1 mm accept 2 1 @@ -404,6 +404,11 @@ expect 0 0 accept 25 25 expect 0.3821 0.4216 +------------------------------------------------------------------------------- +operation +proj=airy +R=1 +no_cut +------------------------------------------------------------------------------- +accept -180 0 +expect failure errno tolerance_condition =============================================================================== Aitoff @@ -1848,6 +1853,10 @@ expect 4030931.833981509 1323687.864777399 accept -80.000000000 36.000000000 expect 3450764.261536101 -175619.041820732 +# For some reason, does not fail on MacOSX +#accept 60 -45 +#expect failure errno tolerance_condition + direction inverse accept -1800000.000000000 2600000.000000000 expect -157.989285000 64.851559610 @@ -2167,6 +2176,9 @@ expect -1575486.353641554 3442168.342028188 accept -2 -1 expect -1575486.353880283 3234352.695594706 +operation +proj=isea +mode=hex +resolution=31 +accept 0 0 +expect failure =============================================================================== Kavraisky V @@ -4010,6 +4022,16 @@ operation +proj=omerc +lat_1=0.8 +a=6400000 +b=.4 ------------------------------------------------------------------------------- expect failure errno invalid_eccentricity +------------------------------------------------------------------------------- +operation +proj=omerc +lat_1=91 +------------------------------------------------------------------------------- +expect failure errno lat_larger_than_90 + +------------------------------------------------------------------------------- +operation +proj=omerc +lat_2=91 +------------------------------------------------------------------------------- +expect failure errno lat_larger_than_90 + =============================================================================== Ortelius Oval @@ -5430,6 +5452,8 @@ accept -200 100 expect -0.001376321 0.001453641 accept -200 -100 expect -0.001988706 -0.000228872 +accept 0 0 +expect 0 0 ------------------------------------------------------------------------------- operation +proj=tpers +a=6400000 +h=1000000 +tilt=20 diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie index ba146b06..75ba55e9 100644 --- a/test/gie/more_builtins.gie +++ b/test/gie/more_builtins.gie @@ -794,5 +794,64 @@ expect 25 25 25 25 operation +proj=aeqd +R=1 +lat_0=91 expect failure errno lat_larger_than_90 +------------------------------------------------------------------------------- +# cart +------------------------------------------------------------------------------- + +operation +proj=cart +ellps=GRS80 +tolerance 0.001mm + +accept 0 0 0 +expect 6378137 0 0 + +accept 0 90 0 +expect 0 0 6356752.314140347 + +accept 0 -90 0 +expect 0 0 -6356752.314140347 + +accept 90 0 0 +expect 0 6378137 0 + +accept -90 0 0 +expect 0 -6378137 0 + +accept 180 0 0 +expect -6378137 0 0 + +accept -180 0 0 +expect -6378137 0 0 + +# Center of Earth ! +accept 0 0 -6378137 +expect 0 0 0 + +accept 0 90 -6356752.314140347 +expect 0 0 0 + +direction inverse + +accept 6378137 0 0 +expect 0 0 0 + +accept 0 0 6356752.314140347 +expect 0 90 0 + +accept 0 0 -6356752.314140347 +expect 0 -90 0 + +accept 0 6378137 0 +expect 90 0 0 + +accept 0 -6378137 0 +expect -90 0 0 + +accept -6378137 0 0 +expect 180 0 0 + +# Center of Earth ! +accept 0 0 0 +expect 0 0 -6378137 + </gie> diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 30e0b427..df2481e6 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -8775,6 +8775,39 @@ TEST(io, createFromUserInput) { EXPECT_THROW( createFromUserInput("urn:ogc:def:unhandled:EPSG::4326", dbContext), ParsingException); + EXPECT_THROW(createFromUserInput("urn:ogc:def:crs:non_existing_auth::4326", + dbContext), + NoSuchAuthorityCodeException); + EXPECT_THROW(createFromUserInput( + "urn:ogc:def:crs,crs:EPSG::2393,unhandled_type:EPSG::5717", + dbContext), + ParsingException); + EXPECT_THROW(createFromUserInput( + "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::unexisting_code", + dbContext), + NoSuchAuthorityCodeException); + EXPECT_THROW( + createFromUserInput( + "urn:ogc:def:crs,crs:EPSG::2393::extra_element,crs:EPSG::EPSG", + dbContext), + ParsingException); + EXPECT_THROW(createFromUserInput("urn:ogc:def:coordinateOperation," + "coordinateOperation:EPSG::3895," + "unhandled_type:EPSG::1618", + dbContext), + ParsingException); + EXPECT_THROW( + createFromUserInput("urn:ogc:def:coordinateOperation," + "coordinateOperation:EPSG::3895," + "coordinateOperation:EPSG::unexisting_code", + dbContext), + NoSuchAuthorityCodeException); + EXPECT_THROW( + createFromUserInput("urn:ogc:def:coordinateOperation," + "coordinateOperation:EPSG::3895::extra_element," + "coordinateOperation:EPSG::1618", + dbContext), + ParsingException); EXPECT_NO_THROW(createFromUserInput("+proj=longlat", nullptr)); EXPECT_NO_THROW(createFromUserInput("EPSG:4326", dbContext)); @@ -8789,6 +8822,31 @@ TEST(io, createFromUserInput) { createFromUserInput("urn:ogc:def:meridian:EPSG::8901", dbContext)); EXPECT_NO_THROW( createFromUserInput("urn:ogc:def:ellipsoid:EPSG::7030", dbContext)); + { + auto obj = createFromUserInput("EPSG:2393+5717", dbContext); + auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), + "KKJ / Finland Uniform Coordinate System + N60 height"); + } + { + auto obj = createFromUserInput( + "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717", dbContext); + auto crs = nn_dynamic_pointer_cast<CompoundCRS>(obj); + ASSERT_TRUE(crs != nullptr); + EXPECT_EQ(crs->nameStr(), + "KKJ / Finland Uniform Coordinate System + N60 height"); + } + { + auto obj = createFromUserInput("urn:ogc:def:coordinateOperation," + "coordinateOperation:EPSG::3895," + "coordinateOperation:EPSG::1618", + dbContext); + auto concat = nn_dynamic_pointer_cast<ConcatenatedOperation>(obj); + ASSERT_TRUE(concat != nullptr); + EXPECT_EQ(concat->nameStr(), + "MGI (Ferro) to MGI (1) + MGI to WGS 84 (3)"); + } EXPECT_NO_THROW(createFromUserInput( "GEOGCRS[\"WGS 84\",\n" " DATUM[\"World Geodetic System 1984\",\n" diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index b7b87d76..a38e9df2 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -5396,6 +5396,39 @@ TEST(operation, projCRS_to_projCRS_context_incompatible_areas) { // --------------------------------------------------------------------------- +TEST(operation, projCRS_to_projCRS_context_incompatible_areas_ballpark) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + auto list = CoordinateOperationFactory::create()->createOperations( + authFactory->createCoordinateReferenceSystem("26711"), // UTM 11 NAD27 + authFactory->createCoordinateReferenceSystem( + "3034"), // ETRS89 / LCC Europe + ctxt); + ASSERT_GE(list.size(), 1U); + EXPECT_TRUE(list[0]->hasBallparkTransformation()); +} + +// --------------------------------------------------------------------------- + +TEST( + operation, + projCRS_to_projCRS_context_incompatible_areas_crs_extent_use_intersection) { + auto authFactory = + AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto ctxt = CoordinateOperationContext::create(authFactory, nullptr, 0.0); + ctxt->setSourceAndTargetCRSExtentUse( + CoordinateOperationContext::SourceTargetCRSExtentUse::INTERSECTION); + auto list = CoordinateOperationFactory::create()->createOperations( + authFactory->createCoordinateReferenceSystem("26711"), // UTM 11 NAD27 + authFactory->createCoordinateReferenceSystem( + "3034"), // ETRS89 / LCC Europe + ctxt); + ASSERT_GE(list.size(), 0U); +} + +// --------------------------------------------------------------------------- + TEST(operation, projCRS_to_projCRS_north_pole_inverted_axis) { auto authFactory = @@ -6444,6 +6477,18 @@ TEST(operation, compoundCRS_to_compoundCRS_context) { "+step +proj=hgridshift +grids=conus +step " "+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap " "+order=2,1"); + { + // Test that we can round-trip this through WKT and still get the same + // PROJ string. + auto wkt = list[0]->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT2_2018).get()); + auto obj = WKTParser().createFromWKT(wkt); + auto co = nn_dynamic_pointer_cast<CoordinateOperation>(obj); + ASSERT_TRUE(co != nullptr); + EXPECT_EQ( + list[0]->exportToPROJString(PROJStringFormatter::create().get()), + co->exportToPROJString(PROJStringFormatter::create().get())); + } bool foundApprox = false; for (size_t i = 0; i < list.size(); i++) { diff --git a/travis/install.sh b/travis/install.sh index c8da9046..3c1e92da 100755 --- a/travis/install.sh +++ b/travis/install.sh @@ -2,6 +2,17 @@ set -e +UNAME="$(uname)" || UNAME="" +if test "${UNAME}" = "Linux" ; then + NPROC=$(nproc); +elif test "${UNAME}" = "Darwin" ; then + NPROC=$(sysctl -n hw.ncpu); +fi +if test "x${NPROC}" = "x"; then + NPROC=2; +fi +echo "NPROC=${NPROC}" + # Download grid files wget http://download.osgeo.org/proj/proj-datumgrid-1.8.zip @@ -27,7 +38,7 @@ if [ -f /usr/lib/jvm/java-7-openjdk-amd64/include/jni.h ]; then else ../configure --prefix=/tmp/proj_autoconf_install_from_dist_all fi -make -j2 +make -j${NPROC} make check make install find /tmp/proj_autoconf_install_from_dist_all @@ -37,7 +48,7 @@ cd .. mkdir build_cmake cd build_cmake cmake .. -DCMAKE_INSTALL_PREFIX=/tmp/proj_cmake_install -VERBOSE=1 make -j2 +VERBOSE=1 make -j${NPROC} make install # The cmake build is not able to generate the null file, so copy it at hand cp /tmp/proj_autoconf_install_from_dist_all/share/proj/null /tmp/proj_cmake_install/share/proj @@ -55,7 +66,7 @@ cd ../.. mkdir build_autoconf_grids cd build_autoconf_grids ../configure --prefix=/tmp/proj_autoconf_install_grids -make -j2 +make -j${NPROC} make check (cd src && make multistresstest && make test228) PROJ_LIB=../data src/multistresstest @@ -74,7 +85,7 @@ if [ "$BUILD_NAME" != "linux_clang" ]; then else ./configure fi -make -j2 +make -j${NPROC} make check # Rerun tests without grids not included in proj-datumgrid diff --git a/travis/osx/before_install.sh b/travis/osx/before_install.sh index 3e9d8d46..a5d40af5 100755 --- a/travis/osx/before_install.sh +++ b/travis/osx/before_install.sh @@ -10,6 +10,7 @@ brew install sqlite3 brew install doxygen export PATH=$HOME/Library/Python/2.7/bin:$PATH -pip install --user sphinx sphinx-rtd-theme sphinxcontrib-bibtex breathe +# breathe=4.12.0 is the last version to work for us with sphinx 1.8.5 / Python 2 +pip install --user sphinx sphinx-rtd-theme sphinxcontrib-bibtex breathe==4.12.0 which sphinx-build (cd docs; make html) |