diff options
38 files changed, 1535 insertions, 1101 deletions
diff --git a/configure.ac b/configure.ac index f0caf4d8..1135b4a4 100644 --- a/configure.ac +++ b/configure.ac @@ -353,7 +353,6 @@ if ! test "x$with_external_gtest" = "xyes" ; then test/googletest/include/gtest/internal/custom/Makefile test/googletest/src/Makefile]) fi -AC_CONFIG_FILES([data/install], [chmod +x data/install]) AC_CONFIG_FILES([proj.pc]) diff --git a/data/install.in b/data/install.in deleted file mode 100755 index 98b2c846..00000000 --- a/data/install.in +++ /dev/null @@ -1,27 +0,0 @@ -: -# SCCSID @(#)install.in 4.8 95/09/24 GIE REL -# -# Installation script all PROJ.4 system reference files. -# -# This script must be executed after compilation of proj library. -# -# Execute as: -# -# install -# -# ------------------------------------------------------------ -prefix=/usr/local -# -lib=${prefix}/lib/proj.4 -# Copy in "init" files -for x in proj_def.dat nad27 nad83 world GL27 -do - cp $x ${lib}/$x - if [ $? -ne 0 ] ; then - echo "init copying failed for file $x" - else - echo "file $x installed" - fi -done -echo "normal completion" -exit 0 diff --git a/docs/source/operations/projections/etmerc.rst b/docs/source/operations/projections/etmerc.rst deleted file mode 100644 index 2a007130..00000000 --- a/docs/source/operations/projections/etmerc.rst +++ /dev/null @@ -1,32 +0,0 @@ -.. _etmerc: - -******************************************************************************** -Extended Transverse Mercator -******************************************************************************** - -.. figure:: ./images/etmerc.png - :width: 500 px - :align: center - :alt: Extended Transverse Mercator - - proj-string: ``+proj=etmerc +lon_0=-40`` - -Parameters -################################################################################ - -.. note:: All parameters for the projection are optional. - -.. include:: ../options/lon_0.rst - -.. include:: ../options/lat_0.rst - -.. include:: ../options/ellps.rst - -.. include:: ../options/R.rst - -.. include:: ../options/k_0.rst - -.. include:: ../options/x_0.rst - -.. include:: ../options/y_0.rst - diff --git a/docs/source/operations/projections/images/etmerc.png b/docs/source/operations/projections/images/etmerc.png Binary files differdeleted file mode 100644 index 35bce083..00000000 --- a/docs/source/operations/projections/images/etmerc.png +++ /dev/null diff --git a/docs/source/operations/projections/index.rst b/docs/source/operations/projections/index.rst index 572deb1a..e6ec9076 100644 --- a/docs/source/operations/projections/index.rst +++ b/docs/source/operations/projections/index.rst @@ -43,7 +43,6 @@ Projections map the spherical 3D space to a flat 2D space. eqdc eqearth euler - etmerc fahey fouc fouc_s diff --git a/docs/source/operations/projections/tmerc.rst b/docs/source/operations/projections/tmerc.rst index 6d1c145d..fd600001 100644 --- a/docs/source/operations/projections/tmerc.rst +++ b/docs/source/operations/projections/tmerc.rst @@ -68,13 +68,19 @@ Example using Gauss-Kruger on Germany area (aka EPSG:31467) :: Example using Gauss Boaga on Italy area (EPSG:3004) :: $ echo 15 42 | proj +proj=tmerc +lat_0=0 +lon_0=15 +k_0=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +no_defs - 2520000.00 4649858.60 + 2520000.00 4649858.60 Parameters ################################################################################ .. note:: All parameters for the projection are optional. +.. option:: +approx + + .. versionadded:: 6.0.0 + + Use faster, less accurate algorithm for the Transverse Mercator. + .. include:: ../options/lon_0.rst .. include:: ../options/lat_0.rst diff --git a/docs/source/operations/projections/utm.rst b/docs/source/operations/projections/utm.rst index ac65aa1a..8759f4eb 100644 --- a/docs/source/operations/projections/utm.rst +++ b/docs/source/operations/projections/utm.rst @@ -66,6 +66,12 @@ Required Add this flag when using the UTM on the southern hemisphere. +.. option:: +approx + + .. versionadded:: 6.0.0 + + Use faster, less accurate algorithm for the Transverse Mercator. + Optional ------------------------------------------------------------------------------- diff --git a/include/proj/io.hpp b/include/proj/io.hpp index 4120d707..091efcb7 100644 --- a/include/proj/io.hpp +++ b/include/proj/io.hpp @@ -365,7 +365,7 @@ class PROJ_GCC_DLL PROJStringFormatter { PROJ_DLL ~PROJStringFormatter(); //! @endcond - PROJ_DLL void setUseETMercForTMerc(bool flag); + PROJ_DLL void setUseApproxTMerc(bool flag); PROJ_DLL const std::string &toString() const; @@ -378,7 +378,7 @@ class PROJ_GCC_DLL PROJStringFormatter { PROJ_DLL void startInversion(); PROJ_DLL void stopInversion(); PROJ_INTERNAL bool isInverted() const; - PROJ_INTERNAL bool getUseETMercForTMerc(bool &settingSetOut) const; + PROJ_INTERNAL bool getUseApproxTMerc() const; PROJ_INTERNAL void setCoordinateOperationOptimizations(bool enable); PROJ_DLL void @@ -888,6 +888,38 @@ class PROJ_GCC_DLL AuthorityFactory { PROJ_DLL std::string getDescriptionText(const std::string &code) const; + /** CRS information */ + struct CRSInfo { + /** Authority name */ + std::string authName{}; + /** Code */ + std::string code{}; + /** Name */ + std::string name{}; + /** Type */ + ObjectType type{ObjectType::CRS}; + /** Whether the object is deprecated */ + bool deprecated{}; + /** Whereas the west_lon_degree, south_lat_degree, east_lon_degree and + * north_lat_degree fields are valid. */ + bool bbox_valid{}; + /** Western-most longitude of the area of use, in degrees. */ + double west_lon_degree{}; + /** Southern-most latitude of the area of use, in degrees. */ + double south_lat_degree{}; + /** Eastern-most longitude of the area of use, in degrees. */ + double east_lon_degree{}; + /** Northern-most latitude of the area of use, in degrees. */ + double north_lat_degree{}; + /** Name of the area of use. */ + std::string areaName{}; + /** Name of the projection method for a projected CRS. Might be empty + * even for projected CRS in some cases. */ + std::string projectionMethodName{}; + }; + + PROJ_DLL std::list<CRSInfo> getCRSInfoList() const; + // non-standard PROJ_DLL static AuthorityFactoryNNPtr create(const DatabaseContextNNPtr &context, diff --git a/m4/pkg.m4 b/m4/pkg.m4 new file mode 100644 index 00000000..452488c8 --- /dev/null +++ b/m4/pkg.m4 @@ -0,0 +1,157 @@ +# pkg.m4 - Macros to locate and utilize pkg-config. -*- Autoconf -*- +# +# Copyright © 2004 Scott James Remnant <scott@netsplit.com>. +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +# +# As a special exception to the GNU General Public License, if you +# distribute this file as part of a program that contains a +# configuration script generated by Autoconf, you may include it under +# the same distribution terms that you use for the rest of that program. + +# PKG_PROG_PKG_CONFIG([MIN-VERSION]) +# ---------------------------------- +AC_DEFUN([PKG_PROG_PKG_CONFIG], +[m4_pattern_forbid([^_?PKG_[A-Z_]+$]) +m4_pattern_allow([^PKG_CONFIG(_PATH)?$]) +AC_ARG_VAR([PKG_CONFIG], [path to pkg-config utility])dnl +if test "x$ac_cv_env_PKG_CONFIG_set" != "xset"; then + AC_PATH_TOOL([PKG_CONFIG], [pkg-config]) +fi +if test -n "$PKG_CONFIG"; then + _pkg_min_version=m4_default([$1], [0.9.0]) + AC_MSG_CHECKING([pkg-config is at least version $_pkg_min_version]) + if $PKG_CONFIG --atleast-pkgconfig-version $_pkg_min_version; then + AC_MSG_RESULT([yes]) + else + AC_MSG_RESULT([no]) + PKG_CONFIG="" + fi + +fi[]dnl +])# PKG_PROG_PKG_CONFIG + +# PKG_CHECK_EXISTS(MODULES, [ACTION-IF-FOUND], [ACTION-IF-NOT-FOUND]) +# +# Check to see whether a particular set of modules exists. Similar +# to PKG_CHECK_MODULES(), but does not set variables or print errors. +# +# +# Similar to PKG_CHECK_MODULES, make sure that the first instance of +# this or PKG_CHECK_MODULES is called, or make sure to call +# PKG_CHECK_EXISTS manually +# -------------------------------------------------------------- +AC_DEFUN([PKG_CHECK_EXISTS], +[AC_REQUIRE([PKG_PROG_PKG_CONFIG])dnl +if test -n "$PKG_CONFIG" && \ + AC_RUN_LOG([$PKG_CONFIG --exists --print-errors "$1"]); then + m4_ifval([$2], [$2], [:]) +m4_ifvaln([$3], [else + $3])dnl +fi]) + + +# _PKG_CONFIG([VARIABLE], [COMMAND], [MODULES]) +# --------------------------------------------- +m4_define([_PKG_CONFIG], +[if test -n "$PKG_CONFIG"; then + if test -n "$$1"; then + pkg_cv_[]$1="$$1" + else + PKG_CHECK_EXISTS([$3], + [pkg_cv_[]$1=`$PKG_CONFIG --[]$2 "$3" 2>/dev/null`], + [pkg_failed=yes]) + fi +else + pkg_failed=untried +fi[]dnl +])# _PKG_CONFIG + +# _PKG_SHORT_ERRORS_SUPPORTED +# ----------------------------- +AC_DEFUN([_PKG_SHORT_ERRORS_SUPPORTED], +[AC_REQUIRE([PKG_PROG_PKG_CONFIG]) +if $PKG_CONFIG --atleast-pkgconfig-version 0.20; then + _pkg_short_errors_supported=yes +else + _pkg_short_errors_supported=no +fi[]dnl +])# _PKG_SHORT_ERRORS_SUPPORTED + + +# PKG_CHECK_MODULES(VARIABLE-PREFIX, MODULES, [ACTION-IF-FOUND], +# [ACTION-IF-NOT-FOUND]) +# +# +# Note that if there is a possibility the first call to +# PKG_CHECK_MODULES might not happen, you should be sure to include an +# explicit call to PKG_PROG_PKG_CONFIG in your configure.ac +# +# +# -------------------------------------------------------------- +AC_DEFUN([PKG_CHECK_MODULES], +[AC_REQUIRE([PKG_PROG_PKG_CONFIG])dnl +AC_ARG_VAR([$1][_CFLAGS], [C compiler flags for $1, overriding pkg-config])dnl +AC_ARG_VAR([$1][_LIBS], [linker flags for $1, overriding pkg-config])dnl + +pkg_failed=no +AC_MSG_CHECKING([for $1]) + +_PKG_CONFIG([$1][_CFLAGS], [cflags], [$2]) +_PKG_CONFIG([$1][_LIBS], [libs], [$2]) + +m4_define([_PKG_TEXT], [Alternatively, you may set the environment variables $1[]_CFLAGS +and $1[]_LIBS to avoid the need to call pkg-config. +See the pkg-config man page for more details.]) + +if test $pkg_failed = yes; then + _PKG_SHORT_ERRORS_SUPPORTED + if test $_pkg_short_errors_supported = yes; then + $1[]_PKG_ERRORS=`$PKG_CONFIG --short-errors --errors-to-stdout --print-errors "$2"` + else + $1[]_PKG_ERRORS=`$PKG_CONFIG --errors-to-stdout --print-errors "$2"` + fi + # Put the nasty error message in config.log where it belongs + echo "$$1[]_PKG_ERRORS" >&AS_MESSAGE_LOG_FD + + ifelse([$4], , [AC_MSG_ERROR(dnl +[Package requirements ($2) were not met: + +$$1_PKG_ERRORS + +Consider adjusting the PKG_CONFIG_PATH environment variable if you +installed software in a non-standard prefix. + +_PKG_TEXT +])], + [AC_MSG_RESULT([no]) + $4]) +elif test $pkg_failed = untried; then + ifelse([$4], , [AC_MSG_FAILURE(dnl +[The pkg-config script could not be found or is too old. Make sure it +is in your PATH or set the PKG_CONFIG environment variable to the full +path to pkg-config. + +_PKG_TEXT + +To get pkg-config, see <http://pkg-config.freedesktop.org/>.])], + [$4]) +else + $1[]_CFLAGS=$pkg_cv_[]$1[]_CFLAGS + $1[]_LIBS=$pkg_cv_[]$1[]_LIBS + AC_MSG_RESULT([yes]) + ifelse([$3], , :, [$3]) +fi[]dnl +])# PKG_CHECK_MODULES diff --git a/src/Makefile.am b/src/Makefile.am index 43bd91df..698baf97 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -122,7 +122,6 @@ libproj_la_SOURCES = \ projections/wag7.cpp \ projections/lcca.cpp \ projections/geos.cpp \ - projections/etmerc.cpp \ projections/boggs.cpp \ projections/collg.cpp \ projections/comill.cpp \ @@ -184,15 +183,15 @@ libproj_la_SOURCES = \ transformations/molodensky.cpp \ transformations/vgridshift.cpp \ \ - aasincos.cpp adjlon.cpp bch2bps.cpp bchgen.cpp \ - biveval.cpp dmstor.cpp mk_cheby.cpp auth.cpp \ + aasincos.cpp adjlon.cpp \ + dmstor.cpp auth.cpp \ deriv.cpp ell_set.cpp ellps.cpp errno.cpp \ factors.cpp fwd.cpp init.cpp inv.cpp \ list.cpp malloc.cpp mlfn.cpp msfn.cpp proj_mdist.cpp \ open_lib.cpp param.cpp phi2.cpp pr_list.cpp \ qsfn.cpp strerrno.cpp \ tsfn.cpp units.cpp ctx.cpp log.cpp zpoly1.cpp rtodms.cpp \ - vector1.cpp release.cpp gauss.cpp \ + release.cpp gauss.cpp \ fileapi.cpp \ \ gc_reader.cpp gridcatalog.cpp \ diff --git a/src/apps/cct.cpp b/src/apps/cct.cpp index 4deefba6..34bf0777 100644 --- a/src/apps/cct.cpp +++ b/src/apps/cct.cpp @@ -410,6 +410,8 @@ int main(int argc, char **argv) { ); } + proj_destroy(P); + if (stdout != fout) fclose (fout); free (o); diff --git a/src/apps/gie.cpp b/src/apps/gie.cpp index e912a076..c3622d52 100644 --- a/src/apps/gie.cpp +++ b/src/apps/gie.cpp @@ -1058,7 +1058,7 @@ static const struct errno_vs_err_const lookup[] = { {"pjd_err_lat_0_or_alpha_eq_90" , -33}, {"pjd_err_ellipsoid_use_required" , -34}, {"pjd_err_invalid_utm_zone" , -35}, - {"pjd_err_tcheby_val_out_of_range" , -36}, + {"" , -36}, /* no longer used */ {"pjd_err_failed_to_find_proj" , -37}, {"pjd_err_failed_to_load_grid" , -38}, {"pjd_err_invalid_m_or_n" , -39}, @@ -1081,6 +1081,7 @@ static const struct errno_vs_err_const lookup[] = { {"pjd_err_ellipsoidal_unsupported" , -56}, {"pjd_err_too_many_inits" , -57}, {"pjd_err_invalid_arg" , -58}, + {"pjd_err_inconsistent_unit" , -59}, {"pjd_err_dont_skip" , 5555}, {"pjd_err_unknown" , 9999}, {"pjd_err_enomem" , ENOMEM}, @@ -1140,7 +1141,8 @@ static int errno_from_err_const (const char *err_const) { /* First try to find a match excluding the PJD_ERR_ prefix */ for (i = 0; i < n; i++) { - if (0==strncmp (lookup[i].the_err_const + 8, err_const, len)) + if (strlen(lookup[i].the_err_const) > 8 && + 0==strncmp (lookup[i].the_err_const + 8, err_const, len)) return lookup[i].the_errno; } diff --git a/src/bch2bps.cpp b/src/bch2bps.cpp deleted file mode 100644 index 9346457c..00000000 --- a/src/bch2bps.cpp +++ /dev/null @@ -1,154 +0,0 @@ -/* convert bivariate w Chebyshev series to w Power series */ - -#include "proj.h" -#include "proj_internal.h" -/* basic support procedures */ - static void /* clear vector to zero */ -clear(PJ_UV *p, int n) { static const PJ_UV c = {0., 0.}; while (n--) *p++ = c; } - static void /* clear matrix rows to zero */ -bclear(PJ_UV **p, int n, int m) { while (n--) clear(*p++, m); } - static void /* move vector */ -bmove(PJ_UV *a, PJ_UV *b, int n) { while (n--) *a++ = *b++; } - static void /* a <- m * b - c */ -submop(PJ_UV *a, double m, PJ_UV *b, PJ_UV *c, int n) { - while (n--) { - a->u = m * b->u - c->u; - a++->v = m * b++->v - c++->v; - } -} - static void /* a <- b - c */ -subop(PJ_UV *a, PJ_UV *b, PJ_UV *c, int n) { - while (n--) { - a->u = b->u - c->u; - a++->v = b++->v - c++->v; - } -} - static void /* multiply vector a by scalar m */ -dmult(PJ_UV *a, double m, int n) { while(n--) { a->u *= m; a->v *= m; ++a; } } - static void /* row adjust a[] <- a[] - m * b[] */ -dadd(PJ_UV *a, PJ_UV *b, double m, int n) { - while(n--) { - a->u -= m * b->u; - a++->v -= m * b++->v; - } -} - static int /* convert row to power series */ -rows(PJ_UV *c, PJ_UV *d, int n) { - PJ_UV sv, *dd; - int j, k; - - dd = (PJ_UV *)vector1(n-1, sizeof(PJ_UV)); - if (!dd) - return 0; - sv.u = sv.v = 0.; - for (j = 0; j < n; ++j) d[j] = dd[j] = sv; - d[0] = c[n-1]; - for (j = n-2; j >= 1; --j) { - for (k = n-j; k >= 1; --k) { - sv = d[k]; - d[k].u = 2. * d[k-1].u - dd[k].u; - d[k].v = 2. * d[k-1].v - dd[k].v; - dd[k] = sv; - } - sv = d[0]; - d[0].u = -dd[0].u + c[j].u; - d[0].v = -dd[0].v + c[j].v; - dd[0] = sv; - } - for (j = n-1; j >= 1; --j) { - d[j].u = d[j-1].u - dd[j].u; - d[j].v = d[j-1].v - dd[j].v; - } - d[0].u = -dd[0].u + .5 * c[0].u; - d[0].v = -dd[0].v + .5 * c[0].v; - pj_dalloc(dd); - return 1; -} - static int /* convert columns to power series */ -cols(PJ_UV **c, PJ_UV **d, int nu, int nv) { - PJ_UV *sv, **dd; - int j, k; - - dd = (PJ_UV **)vector2(nu, nv, sizeof(PJ_UV)); - if (!dd) - return 0; - sv = (PJ_UV *)vector1(nv, sizeof(PJ_UV)); - if (!sv) { - freev2((void **)dd, nu); - return 0; - } - bclear(d, nu, nv); - bclear(dd, nu, nv); - bmove(d[0], c[nu-1], nv); - for (j = nu-2; j >= 1; --j) { - for (k = nu-j; k >= 1; --k) { - bmove(sv, d[k], nv); - submop(d[k], 2., d[k-1], dd[k], nv); - bmove(dd[k], sv, nv); - } - bmove(sv, d[0], nv); - subop(d[0], c[j], dd[0], nv); - bmove(dd[0], sv, nv); - } - for (j = nu-1; j >= 1; --j) - subop(d[j], d[j-1], dd[j], nv); - submop(d[0], .5, c[0], dd[0], nv); - freev2((void **) dd, nu); - pj_dalloc(sv); - return 1; -} - static void /* row adjust for range -1 to 1 to a to b */ -rowshft(double a, double b, PJ_UV *d, int n) { - int k, j; - double fac, cnst; - - cnst = 2. / (b - a); - fac = cnst; - for (j = 1; j < n; ++j) { - d[j].u *= fac; - d[j].v *= fac; - fac *= cnst; - } - cnst = .5 * (a + b); - for (j = 0; j <= n-2; ++j) - for (k = n - 2; k >= j; --k) { - d[k].u -= cnst * d[k+1].u; - d[k].v -= cnst * d[k+1].v; - } -} - static void /* column adjust for range -1 to 1 to a to b */ -colshft(double a, double b, PJ_UV **d, int n, int m) { - int k, j; - double fac, cnst; - - cnst = 2. / (b - a); - fac = cnst; - for (j = 1; j < n; ++j) { - dmult(d[j], fac, m); - fac *= cnst; - } - cnst = .5 * (a + b); - for (j = 0; j <= n-2; ++j) - for (k = n - 2; k >= j; --k) - dadd(d[k], d[k+1], cnst, m); -} - int /* entry point */ -bch2bps(PJ_UV a, PJ_UV b, PJ_UV **c, int nu, int nv) { - PJ_UV **d; - int i; - - if (nu < 1 || nv < 1 || !(d = (PJ_UV **)vector2(nu, nv, sizeof(PJ_UV)))) - return 0; - /* do rows to power series */ - for (i = 0; i < nu; ++i) { - if (!rows(c[i], d[i], nv)) - return 0; - rowshft(a.v, b.v, d[i], nv); - } - /* do columns to power series */ - if (!cols(d, c, nu, nv)) - return 0; - colshft(a.u, b.u, c, nu, nv); - freev2((void **) d, nu); - return 1; -} diff --git a/src/bchgen.cpp b/src/bchgen.cpp deleted file mode 100644 index 9677b6f2..00000000 --- a/src/bchgen.cpp +++ /dev/null @@ -1,59 +0,0 @@ -/* generate double bivariate Chebychev polynomial */ -#include "proj.h" -#include "proj_internal.h" - int -bchgen(PJ_UV a, PJ_UV b, int nu, int nv, PJ_UV **f, PJ_UV(*func)(PJ_UV)) { - int i, j, k; - PJ_UV arg, *t, bma, bpa, *c; - double d, fac; - - bma.u = 0.5 * (b.u - a.u); bma.v = 0.5 * (b.v - a.v); - bpa.u = 0.5 * (b.u + a.u); bpa.v = 0.5 * (b.v + a.v); - for ( i = 0; i < nu; ++i) { - arg.u = cos(M_PI * (i + 0.5) / nu) * bma.u + bpa.u; - for ( j = 0; j < nv; ++j) { - arg.v = cos(M_PI * (j + 0.5) / nv) * bma.v + bpa.v; - f[i][j] = (*func)(arg); - if ((f[i][j]).u == HUGE_VAL) - return(1); - } - } - if (!(c = (PJ_UV *) vector1(nu, sizeof(PJ_UV)))) return 1; - fac = 2. / nu; - for ( j = 0; j < nv ; ++j) { - for ( i = 0; i < nu; ++i) { - arg.u = arg.v = 0.; - for (k = 0; k < nu; ++k) { - d = cos(M_PI * i * (k + .5) / nu); - arg.u += f[k][j].u * d; - arg.v += f[k][j].v * d; - } - arg.u *= fac; - arg.v *= fac; - c[i] = arg; - } - for (i = 0; i < nu; ++i) - f[i][j] = c[i]; - } - pj_dalloc(c); - if (!(c = (PJ_UV*) vector1(nv, sizeof(PJ_UV)))) return 1; - fac = 2. / nv; - for ( i = 0; i < nu; ++i) { - t = f[i]; - for (j = 0; j < nv; ++j) { - arg.u = arg.v = 0.; - for (k = 0; k < nv; ++k) { - d = cos(M_PI * j * (k + .5) / nv); - arg.u += t[k].u * d; - arg.v += t[k].v * d; - } - arg.u *= fac; - arg.v *= fac; - c[j] = arg; - } - f[i] = c; - c = t; - } - pj_dalloc(c); - return(0); -} diff --git a/src/biveval.cpp b/src/biveval.cpp deleted file mode 100644 index 9ead2fb7..00000000 --- a/src/biveval.cpp +++ /dev/null @@ -1,88 +0,0 @@ -/* procedures for evaluating Tseries */ -#include "proj.h" -#include "proj_internal.h" -# define NEAR_ONE 1.00001 -static double ceval(struct PW_COEF *C, int n, PJ_UV w, PJ_UV w2) { - double d=0, dd=0, vd, vdd, tmp, *c; - int j; - - for (C += n ; n-- ; --C ) { - if ( (j = C->m) != 0) { - vd = vdd = 0.; - for (c = C->c + --j; j ; --j ) { - vd = w2.v * (tmp = vd) - vdd + *c--; - vdd = tmp; - } - d = w2.u * (tmp = d) - dd + w.v * vd - vdd + 0.5 * *c; - } else - d = w2.u * (tmp = d) - dd; - dd = tmp; - } - if ( (j = C->m) != 0 ) { - vd = vdd = 0.; - for (c = C->c + --j; j ; --j ) { - vd = w2.v * (tmp = vd) - vdd + *c--; - vdd = tmp; - } - return (w.u * d - dd + 0.5 * ( w.v * vd - vdd + 0.5 * *c )); - } else - return (w.u * d - dd); -} - -PJ_UV /* bivariate Chebyshev polynomial entry point */ -bcheval(PJ_UV in, Tseries *T) { - PJ_UV w2, w; - PJ_UV out; - /* scale to +-1 */ - w.u = ( in.u + in.u - T->a.u ) * T->b.u; - w.v = ( in.v + in.v - T->a.v ) * T->b.v; - if (fabs(w.u) > NEAR_ONE || fabs(w.v) > NEAR_ONE) { - out.u = out.v = HUGE_VAL; - pj_errno = PJD_ERR_TCHEBY_VAL_OUT_OF_RANGE; - } else { /* double evaluation */ - w2.u = w.u + w.u; - w2.v = w.v + w.v; - out.u = ceval(T->cu, T->mu, w, w2); - out.v = ceval(T->cv, T->mv, w, w2); - } - return out; -} - -PJ_UV /* bivariate power polynomial entry point */ -bpseval(PJ_UV in, Tseries *T) { - PJ_UV out; - double *c, row; - int i, m; - - out.u = out.v = 0.; - for (i = T->mu; i >= 0; --i) { - row = 0.; - if ((m = T->cu[i].m) != 0) { - c = T->cu[i].c + m; - while (m--) - row = *--c + in.v * row; - } - out.u = row + in.u * out.u; - } - for (i = T->mv; i >= 0; --i) { - row = 0.; - if ((m = T->cv[i].m) != 0) { - c = T->cv[i].c + m; - while (m--) - row = *--c + in.v * row; - } - out.v = row + in.u * out.v; - } - return out; -} - -PJ_UV /* general entry point selecting evaluation mode */ -biveval(PJ_UV in, Tseries *T) { - - if (T->power) { - return bpseval(in, T); - } else { - return bcheval(in, T); - } -} - diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index de11f181..240d11b3 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -34,6 +34,7 @@ #include <cstdarg> #include <cstring> #include <map> +#include <new> #include <utility> #include <vector> @@ -144,9 +145,11 @@ static PJ *pj_obj_create(PJ_CONTEXT *ctx, const IdentifiedObjectNNPtr &objIn) { // PROJ string. } } - auto pj = new PJ(); - pj->descr = "ISO-19111 object"; - pj->iso_obj = objIn; + auto pj = pj_new(); + if (pj) { + pj->descr = "ISO-19111 object"; + pj->iso_obj = objIn; + } return pj; } //! @endcond @@ -1204,9 +1207,8 @@ const char *proj_as_wkt(PJ_CONTEXT *ctx, const PJ *obj, PJ_WKT_TYPE type, * @param type PROJ String version. * @param options NULL-terminated list of strings with "KEY=VALUE" format. or * NULL. - * The currently recognized option is USE_ETMERC=YES to use - * +proj=etmerc instead of +proj=tmerc (or USE_ETMERC=NO to disable implicit - * use of etmerc by utm conversions) + * The currently recognized option is USE_APPROX_TMERC=YES to add the +approx + * flag to +proj=tmerc or +proj=utm * @return a string, or NULL in case of error. */ const char *proj_as_proj_string(PJ_CONTEXT *ctx, const PJ *obj, @@ -1240,10 +1242,8 @@ const char *proj_as_proj_string(PJ_CONTEXT *ctx, const PJ *obj, try { auto formatter = PROJStringFormatter::create(convention, dbContext); if (options != nullptr && options[0] != nullptr) { - if (ci_equal(options[0], "USE_ETMERC=YES")) { - formatter->setUseETMercForTMerc(true); - } else if (ci_equal(options[0], "USE_ETMERC=NO")) { - formatter->setUseETMercForTMerc(false); + if (ci_equal(options[0], "USE_APPROX_TMERC=YES")) { + formatter->setUseApproxTMerc(true); } } obj->lastPROJString = exportable->exportToPROJString(formatter.get()); @@ -1895,6 +1895,8 @@ PROJ_STRING_LIST proj_get_authorities_from_database(PJ_CONTEXT *ctx) { * * @return a NULL terminated list of NUL-terminated strings that must be * freed with proj_string_list_destroy(), or NULL in case of error. + * + * @see proj_get_crs_info_list_from_database() */ PROJ_STRING_LIST proj_get_codes_from_database(PJ_CONTEXT *ctx, const char *auth_name, @@ -1932,6 +1934,215 @@ void proj_string_list_destroy(PROJ_STRING_LIST list) { // --------------------------------------------------------------------------- +/** \brief Instanciate a default set of parameters to be used by + * proj_get_crs_list(). + * + * @return a new object to free with proj_get_crs_list_parameters_destroy() */ +PROJ_CRS_LIST_PARAMETERS *proj_get_crs_list_parameters_create() { + auto ret = new (std::nothrow) PROJ_CRS_LIST_PARAMETERS(); + if (ret) { + ret->types = nullptr; + ret->typesCount = 0; + ret->crs_area_of_use_contains_bbox = TRUE; + ret->bbox_valid = FALSE; + ret->west_lon_degree = 0.0; + ret->south_lat_degree = 0.0; + ret->east_lon_degree = 0.0; + ret->north_lat_degree = 0.0; + ret->allow_deprecated = FALSE; + } + return ret; +} + +// --------------------------------------------------------------------------- + +/** \brief Destroy an object returned by proj_get_crs_list_parameters_create() + */ +void proj_get_crs_list_parameters_destroy(PROJ_CRS_LIST_PARAMETERS *params) { + delete params; +} + +// --------------------------------------------------------------------------- + +/** \brief Enumerate CRS objects from the database, taking into account various + * criteria. + * + * The returned object is an array of PROJ_CRS_INFO* pointers, whose last + * entry is NULL. This array should be freed with proj_crs_info_list_destroy() + * + * When no filter parameters are set, this is functionnaly equivalent to + * proj_get_crs_info_list_from_database(), instanciating a PJ* object for each + * of the proj_create_from_database() and retrieving information with the + * various getters. However this function will be much faster. + * + * @param ctx PROJ context, or NULL for default context + * @param auth_name Authority name, used to restrict the search. + * Or NULL for all authorities. + * @param params Additional criteria, or NULL. If not-NULL, params SHOULD + * have been allocated by proj_get_crs_list_parameters_create(), as the + * PROJ_CRS_LIST_PARAMETERS structure might grow over time. + * @param out_result_count Output parameter pointing to an integer to receive + * the size of the result list. Might be NULL + * @return an array of PROJ_CRS_INFO* pointers to be freed with + * proj_crs_info_list_destroy(), or NULL in case of error. + */ +PROJ_CRS_INFO ** +proj_get_crs_info_list_from_database(PJ_CONTEXT *ctx, const char *auth_name, + const PROJ_CRS_LIST_PARAMETERS *params, + int *out_result_count) { + SANITIZE_CTX(ctx); + PROJ_CRS_INFO **ret = nullptr; + int i = 0; + try { + auto factory = AuthorityFactory::create(getDBcontext(ctx), + auth_name ? auth_name : ""); + auto list = factory->getCRSInfoList(); + ret = new PROJ_CRS_INFO *[list.size() + 1]; + GeographicBoundingBoxPtr bbox; + if (params && params->bbox_valid) { + bbox = GeographicBoundingBox::create( + params->west_lon_degree, params->south_lat_degree, + params->east_lon_degree, params->north_lat_degree) + .as_nullable(); + } + for (const auto &info : list) { + auto type = PJ_TYPE_CRS; + if (info.type == AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS) { + type = PJ_TYPE_GEOGRAPHIC_2D_CRS; + } else if (info.type == + AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS) { + type = PJ_TYPE_GEOGRAPHIC_3D_CRS; + } else if (info.type == + AuthorityFactory::ObjectType::GEOCENTRIC_CRS) { + type = PJ_TYPE_GEOCENTRIC_CRS; + } else if (info.type == + AuthorityFactory::ObjectType::PROJECTED_CRS) { + type = PJ_TYPE_PROJECTED_CRS; + } else if (info.type == + AuthorityFactory::ObjectType::VERTICAL_CRS) { + type = PJ_TYPE_VERTICAL_CRS; + } else if (info.type == + AuthorityFactory::ObjectType::COMPOUND_CRS) { + type = PJ_TYPE_COMPOUND_CRS; + } + if (params && params->typesCount) { + bool typeValid = false; + for (size_t j = 0; j < params->typesCount; j++) { + if (params->types[j] == type) { + typeValid = true; + break; + } else if (params->types[j] == PJ_TYPE_GEOGRAPHIC_CRS && + (type == PJ_TYPE_GEOGRAPHIC_2D_CRS || + type == PJ_TYPE_GEOGRAPHIC_3D_CRS)) { + typeValid = true; + break; + } else if (params->types[j] == PJ_TYPE_GEODETIC_CRS && + (type == PJ_TYPE_GEOCENTRIC_CRS || + type == PJ_TYPE_GEOGRAPHIC_2D_CRS || + type == PJ_TYPE_GEOGRAPHIC_3D_CRS)) { + typeValid = true; + break; + } + } + if (!typeValid) { + continue; + } + } + if (params && !params->allow_deprecated && info.deprecated) { + continue; + } + if (params && params->bbox_valid) { + if (!info.bbox_valid) { + continue; + } + if (info.west_lon_degree <= info.east_lon_degree && + params->west_lon_degree <= params->east_lon_degree) { + if (params->crs_area_of_use_contains_bbox) { + if (params->west_lon_degree < info.west_lon_degree || + params->east_lon_degree > info.east_lon_degree || + params->south_lat_degree < info.south_lat_degree || + params->north_lat_degree > info.north_lat_degree) { + continue; + } + } else { + if (info.east_lon_degree < params->west_lon_degree || + info.west_lon_degree > params->east_lon_degree || + info.north_lat_degree < params->south_lat_degree || + info.south_lat_degree > params->north_lat_degree) { + continue; + } + } + } else { + auto crsExtent = GeographicBoundingBox::create( + info.west_lon_degree, info.south_lat_degree, + info.east_lon_degree, info.north_lat_degree); + if (params->crs_area_of_use_contains_bbox) { + if (!crsExtent->contains(NN_NO_CHECK(bbox))) { + continue; + } + } else { + if (!bbox->intersects(crsExtent)) { + continue; + } + } + } + } + + ret[i] = new PROJ_CRS_INFO; + ret[i]->auth_name = pj_strdup(info.authName.c_str()); + ret[i]->code = pj_strdup(info.code.c_str()); + ret[i]->name = pj_strdup(info.name.c_str()); + ret[i]->type = type; + ret[i]->deprecated = info.deprecated; + ret[i]->bbox_valid = info.bbox_valid; + ret[i]->west_lon_degree = info.west_lon_degree; + ret[i]->south_lat_degree = info.south_lat_degree; + ret[i]->east_lon_degree = info.east_lon_degree; + ret[i]->north_lat_degree = info.north_lat_degree; + ret[i]->area_name = pj_strdup(info.areaName.c_str()); + ret[i]->projection_method_name = + info.projectionMethodName.empty() + ? nullptr + : pj_strdup(info.projectionMethodName.c_str()); + i++; + } + ret[i] = nullptr; + if (out_result_count) + *out_result_count = i; + return ret; + } catch (const std::exception &e) { + proj_log_error(ctx, __FUNCTION__, e.what()); + if (ret) { + ret[i + 1] = nullptr; + proj_crs_info_list_destroy(ret); + } + if (out_result_count) + *out_result_count = 0; + } + return nullptr; +} + +// --------------------------------------------------------------------------- + +/** \brief Destroy the result returned by + * proj_get_crs_info_list_from_database(). + */ +void proj_crs_info_list_destroy(PROJ_CRS_INFO **list) { + if (list) { + for (int i = 0; list[i] != nullptr; i++) { + pj_dalloc(list[i]->auth_name); + pj_dalloc(list[i]->code); + pj_dalloc(list[i]->name); + pj_dalloc(list[i]->area_name); + pj_dalloc(list[i]->projection_method_name); + delete list[i]; + } + delete[] list; + } +} + +// --------------------------------------------------------------------------- + /** \brief Return the Conversion of a DerivedCRS (such as a ProjectedCRS), * or the Transformation from the baseCRS to the hubCRS of a BoundCRS * diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 8a10bc5a..d36d901f 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -5381,13 +5381,11 @@ bool Conversion::addWKTExtensionNode(io::WKTFormatter *formatter) const { const auto &l_method = method(); const auto &methodName = l_method->nameStr(); const int methodEPSGCode = l_method->getEPSGCode(); - int zone = 0; - bool north = true; - if (l_method->getPrivate()->projMethodOverride_ == "etmerc" && - !isUTM(zone, north)) { + if (l_method->getPrivate()->projMethodOverride_ == "tmerc approx" || + l_method->getPrivate()->projMethodOverride_ == "utm approx") { auto projFormatter = io::PROJStringFormatter::create(); projFormatter->setCRSExport(true); - projFormatter->setUseETMercForTMerc(true); + projFormatter->setUseApproxTMerc(true); formatter->startNode(io::WKTConstants::EXTENSION, false); formatter->addQuotedString("PROJ4"); _exportToPROJString(projFormatter.get()); @@ -5479,17 +5477,21 @@ void Conversion::_exportToPROJString( const auto &convName = nameStr(); bool bConversionDone = false; bool bEllipsoidParametersDone = false; - bool useETMerc = false; + bool useApprox = false; if (methodEPSGCode == EPSG_CODE_METHOD_TRANSVERSE_MERCATOR) { // Check for UTM int zone = 0; bool north = true; - bool etMercSettingSet = false; - useETMerc = formatter->getUseETMercForTMerc(etMercSettingSet) || - l_method->getPrivate()->projMethodOverride_ == "etmerc"; - if (isUTM(zone, north) && !(etMercSettingSet && !useETMerc)) { + useApprox = + formatter->getUseApproxTMerc() || + l_method->getPrivate()->projMethodOverride_ == "tmerc approx" || + l_method->getPrivate()->projMethodOverride_ == "utm approx"; + if (isUTM(zone, north)) { bConversionDone = true; formatter->addStep("utm"); + if( useApprox ) { + formatter->addParam("approx"); + } formatter->addParam("zone", zone); if (!north) { formatter->addParam("south"); @@ -5662,7 +5664,10 @@ void Conversion::_exportToPROJString( if (!bConversionDone) { const MethodMapping *mapping = getMapping(l_method.get()); if (mapping && mapping->proj_name_main) { - formatter->addStep(useETMerc ? "etmerc" : mapping->proj_name_main); + formatter->addStep(mapping->proj_name_main); + if (useApprox) { + formatter->addParam("approx"); + } if (mapping->proj_name_aux) { if (internal::starts_with(mapping->proj_name_aux, "axis=")) { bAxisSpecFound = true; diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index c862cc5b..81abcdf1 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -77,9 +77,17 @@ namespace io { //! @cond Doxygen_Suppress -#define GEOG_2D "'geographic 2D'" -#define GEOG_3D "'geographic 3D'" -#define GEOCENTRIC "'geocentric'" +// CRS subtypes +#define GEOG_2D "geographic 2D" +#define GEOG_3D "geographic 3D" +#define GEOCENTRIC "geocentric" +#define PROJECTED "projected" +#define VERTICAL "vertical" +#define COMPOUND "compound" + +#define GEOG_2D_SINGLE_QUOTED "'geographic 2D'" +#define GEOG_3D_SINGLE_QUOTED "'geographic 3D'" +#define GEOCENTRIC_SINGLE_QUOTED "'geocentric'" // --------------------------------------------------------------------------- @@ -1051,7 +1059,7 @@ DatabaseContext::getAliasFromOfficialName(const std::string &officialName, sql += replaceAll(tableName, "\"", "\"\""); sql += "\" WHERE name = ?"; if (tableName == "geodetic_crs") { - sql += " AND type = " GEOG_2D; + sql += " AND type = " GEOG_2D_SINGLE_QUOTED; } auto res = d->run(sql, {officialName}); if (res.empty()) { @@ -1485,6 +1493,12 @@ AuthorityFactory::createExtent(const std::string &code) const { try { const auto &row = res.front(); const auto &name = row[0]; + if (row[1].empty()) { + auto extent = metadata::Extent::create( + util::optional<std::string>(name), {}, {}, {}); + d->context()->d->cache(code, extent); + return extent; + } double south_lat = c_locale_stod(row[1]); double north_lat = c_locale_stod(row[2]); double west_lon = c_locale_stod(row[3]); @@ -2061,7 +2075,8 @@ AuthorityFactory::createGeodeticCRS(const std::string &code, "deprecated FROM " "geodetic_crs WHERE auth_name = ? AND code = ?"); if (geographicOnly) { - sql += " AND type in (" GEOG_2D "," GEOG_3D ")"; + sql += " AND type in (" GEOG_2D_SINGLE_QUOTED "," GEOG_3D_SINGLE_QUOTED + ")"; } auto res = d->runWithCodeParam(sql, code); if (res.empty()) { @@ -2118,15 +2133,14 @@ AuthorityFactory::createGeodeticCRS(const std::string &code, auto ellipsoidalCS = util::nn_dynamic_pointer_cast<cs::EllipsoidalCS>(cs); - if ((type == "geographic 2D" || type == "geographic 3D") && - ellipsoidalCS) { + if ((type == GEOG_2D || type == GEOG_3D) && ellipsoidalCS) { auto crsRet = crs::GeographicCRS::create( props, datum, NN_NO_CHECK(ellipsoidalCS)); d->context()->d->cache(cacheKey, crsRet); return crsRet; } auto geocentricCS = util::nn_dynamic_pointer_cast<cs::CartesianCS>(cs); - if (type == "geocentric" && geocentricCS) { + if (type == GEOCENTRIC && geocentricCS) { auto crsRet = crs::GeodeticCRS::create(props, datum, NN_NO_CHECK(geocentricCS)); d->context()->d->cache(cacheKey, crsRet); @@ -2471,17 +2485,16 @@ AuthorityFactory::createCoordinateReferenceSystem(const std::string &code, code); } const auto &type = res.front()[0]; - if (type == "geographic 2D" || type == "geographic 3D" || - type == "geocentric") { + if (type == GEOG_2D || type == GEOG_3D || type == GEOCENTRIC) { return createGeodeticCRS(code); } - if (type == "vertical") { + if (type == VERTICAL) { return createVerticalCRS(code); } - if (type == "projected") { + if (type == PROJECTED) { return createProjectedCRS(code); } - if (allowCompound && type == "compound") { + if (allowCompound && type == COMPOUND) { return createCompoundCRS(code); } throw FactoryException("unhandled CRS type: " + type); @@ -3776,17 +3789,22 @@ AuthorityFactory::getAuthorityCodes(const ObjectType &type, sql = "SELECT code FROM geodetic_crs WHERE "; break; case ObjectType::GEOCENTRIC_CRS: - sql = "SELECT code FROM geodetic_crs WHERE type = " GEOCENTRIC " AND "; + sql = "SELECT code FROM geodetic_crs WHERE type " + "= " GEOCENTRIC_SINGLE_QUOTED " AND "; break; case ObjectType::GEOGRAPHIC_CRS: - sql = "SELECT code FROM geodetic_crs WHERE type IN (" GEOG_2D - "," GEOG_3D ") AND "; + sql = "SELECT code FROM geodetic_crs WHERE type IN " + "(" GEOG_2D_SINGLE_QUOTED "," GEOG_3D_SINGLE_QUOTED ") AND "; break; case ObjectType::GEOGRAPHIC_2D_CRS: - sql = "SELECT code FROM geodetic_crs WHERE type = " GEOG_2D " AND "; + sql = + "SELECT code FROM geodetic_crs WHERE type = " GEOG_2D_SINGLE_QUOTED + " AND "; break; case ObjectType::GEOGRAPHIC_3D_CRS: - sql = "SELECT code FROM geodetic_crs WHERE type = " GEOG_3D " AND "; + sql = + "SELECT code FROM geodetic_crs WHERE type = " GEOG_3D_SINGLE_QUOTED + " AND "; break; case ObjectType::VERTICAL_CRS: sql = "SELECT code FROM vertical_crs WHERE "; @@ -3852,6 +3870,107 @@ AuthorityFactory::getDescriptionText(const std::string &code) const { // --------------------------------------------------------------------------- +/** \brief Return a list of information on CRS objects + * + * This is functionnaly equivalent of listing the codes from an authority, + * instanciating + * a CRS object for each of them and getting the information from this CRS + * object, but this implementation has much less overhead. + * + * @throw FactoryException + */ +std::list<AuthorityFactory::CRSInfo> AuthorityFactory::getCRSInfoList() const { + std::string sql = "SELECT c.auth_name, c.code, c.name, c.type, " + "c.deprecated, " + "a.west_lon, a.south_lat, a.east_lon, a.north_lat, " + "a.name, NULL FROM geodetic_crs c " + "JOIN area a ON " + "c.area_of_use_auth_name = a.auth_name AND " + "c.area_of_use_code = a.code"; + ListOfParams params; + if (d->hasAuthorityRestriction()) { + sql += " WHERE c.auth_name = ?"; + params.emplace_back(d->authority()); + } + sql += " UNION ALL "; + sql += "SELECT c.auth_name, c.code, c.name, 'projected', " + "c.deprecated, " + "a.west_lon, a.south_lat, a.east_lon, a.north_lat, " + "a.name, conv.method_name FROM projected_crs c " + "JOIN area a ON " + "c.area_of_use_auth_name = a.auth_name AND " + "c.area_of_use_code = a.code " + "LEFT JOIN conversion conv ON " + "c.conversion_auth_name = conv.auth_name AND " + "c.conversion_code = conv.code"; + if (d->hasAuthorityRestriction()) { + sql += " WHERE c.auth_name = ?"; + params.emplace_back(d->authority()); + } + sql += " UNION ALL "; + sql += "SELECT c.auth_name, c.code, c.name, 'vertical', " + "c.deprecated, " + "a.west_lon, a.south_lat, a.east_lon, a.north_lat, " + "a.name, NULL FROM vertical_crs c " + "JOIN area a ON " + "c.area_of_use_auth_name = a.auth_name AND " + "c.area_of_use_code = a.code"; + if (d->hasAuthorityRestriction()) { + sql += " WHERE c.auth_name = ?"; + params.emplace_back(d->authority()); + } + sql += " UNION ALL "; + sql += "SELECT c.auth_name, c.code, c.name, 'compound', " + "c.deprecated, " + "a.west_lon, a.south_lat, a.east_lon, a.north_lat, " + "a.name, NULL FROM compound_crs c " + "JOIN area a ON " + "c.area_of_use_auth_name = a.auth_name AND " + "c.area_of_use_code = a.code"; + if (d->hasAuthorityRestriction()) { + sql += " WHERE c.auth_name = ?"; + params.emplace_back(d->authority()); + } + auto sqlRes = d->run(sql, params); + std::list<AuthorityFactory::CRSInfo> res; + for (const auto &row : sqlRes) { + AuthorityFactory::CRSInfo info; + info.authName = row[0]; + info.code = row[1]; + info.name = row[2]; + const auto &type = row[3]; + if (type == GEOG_2D) { + info.type = AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS; + } else if (type == GEOG_3D) { + info.type = AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS; + } else if (type == GEOCENTRIC) { + info.type = AuthorityFactory::ObjectType::GEOCENTRIC_CRS; + } else if (type == PROJECTED) { + info.type = AuthorityFactory::ObjectType::PROJECTED_CRS; + } else if (type == VERTICAL) { + info.type = AuthorityFactory::ObjectType::VERTICAL_CRS; + } else if (type == COMPOUND) { + info.type = AuthorityFactory::ObjectType::COMPOUND_CRS; + } + info.deprecated = row[4] == "1"; + if (row[5].empty()) { + info.bbox_valid = false; + } else { + info.bbox_valid = true; + info.west_lon_degree = c_locale_stod(row[5]); + info.south_lat_degree = c_locale_stod(row[6]); + info.east_lon_degree = c_locale_stod(row[7]); + info.north_lat_degree = c_locale_stod(row[8]); + } + info.areaName = row[9]; + info.projectionMethodName = row[10]; + res.emplace_back(info); + } + return res; +} + +// --------------------------------------------------------------------------- + /** \brief Gets the official name from a possibly alias name. * * @param aliasedName Alias name. @@ -4050,23 +4169,25 @@ AuthorityFactory::createObjectsFromName( break; case ObjectType::GEOCENTRIC_CRS: addToListStringWithOR(otherConditions, - "(table_name = " GEOCENTRIC " AND " - "type = " GEOCENTRIC ")"); + "(table_name = " GEOCENTRIC_SINGLE_QUOTED + " AND " + "type = " GEOCENTRIC_SINGLE_QUOTED ")"); break; case ObjectType::GEOGRAPHIC_CRS: addToListStringWithOR(otherConditions, "(table_name = 'geodetic_crs' AND " - "type IN (" GEOG_2D "," GEOG_3D "))"); + "type IN (" GEOG_2D_SINGLE_QUOTED + "," GEOG_3D_SINGLE_QUOTED "))"); break; case ObjectType::GEOGRAPHIC_2D_CRS: addToListStringWithOR(otherConditions, "(table_name = 'geodetic_crs' AND " - "type = " GEOG_2D ")"); + "type = " GEOG_2D_SINGLE_QUOTED ")"); break; case ObjectType::GEOGRAPHIC_3D_CRS: addToListStringWithOR(otherConditions, "(table_name = 'geodetic_crs' AND " - "type = " GEOG_3D ")"); + "type = " GEOG_3D_SINGLE_QUOTED ")"); break; case ObjectType::PROJECTED_CRS: addToListString(tableNameList, "'projected_crs'"); diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index a851b41b..078180a5 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -4746,8 +4746,7 @@ struct PROJStringFormatter::Private { bool omitProjLongLatIfPossible_ = false; bool omitZUnitConversion_ = false; DatabaseContextPtr dbContext_{}; - bool useETMercForTMerc_ = false; - bool useETMercForTMercSet_ = false; + bool useApproxTMerc_ = false; bool addNoDefs_ = true; bool coordOperationOptimizations_ = false; bool crsExport_ = false; @@ -4803,11 +4802,9 @@ PROJStringFormatter::create(Convention conventionIn, // --------------------------------------------------------------------------- -/** \brief Set whether Extended Transverse Mercator (etmerc) should be used - * instead of tmerc */ -void PROJStringFormatter::setUseETMercForTMerc(bool flag) { - d->useETMercForTMerc_ = flag; - d->useETMercForTMercSet_ = true; +/** \brief Set whether approximate Transverse Mercator or UTM should be used */ +void PROJStringFormatter::setUseApproxTMerc(bool flag) { + d->useApproxTMerc_ = flag; } // --------------------------------------------------------------------------- @@ -5256,9 +5253,8 @@ PROJStringFormatter::Convention PROJStringFormatter::convention() const { // --------------------------------------------------------------------------- -bool PROJStringFormatter::getUseETMercForTMerc(bool &settingSetOut) const { - settingSetOut = d->useETMercForTMercSet_; - return d->useETMercForTMerc_; +bool PROJStringFormatter::getUseApproxTMerc() const { + return d->useApproxTMerc_; } // --------------------------------------------------------------------------- @@ -7081,8 +7077,10 @@ CRSNNPtr PROJStringParser::Private::buildProjectedCRS( : UnitOfMeasure::NONE))); } - if (step.name == "etmerc") { - methodMap.set("proj_method", "etmerc"); + if (step.name == "tmerc" && hasParamValue(step, "approx")) { + methodMap.set("proj_method", "tmerc approx"); + } else if (step.name == "utm" && hasParamValue(step, "approx")) { + methodMap.set("proj_method", "utm approx"); } conv = Conversion::create(mapWithUnknownName, methodMap, parameters, diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index b595d87e..f28a1d68 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -117,7 +117,6 @@ SET(SRC_LIBPROJ_PROJECTIONS projections/wag7.cpp projections/lcca.cpp projections/geos.cpp - projections/etmerc.cpp projections/boggs.cpp projections/collg.cpp projections/comill.cpp @@ -201,15 +200,15 @@ SET(SRC_LIBPROJ_ISO19111 SET(SRC_LIBPROJ_CORE pj_list.h proj_internal.h proj_math.h - aasincos.cpp adjlon.cpp bch2bps.cpp bchgen.cpp - biveval.cpp dmstor.cpp mk_cheby.cpp auth.cpp + aasincos.cpp adjlon.cpp + dmstor.cpp auth.cpp deriv.cpp ell_set.cpp ellps.cpp errno.cpp factors.cpp fwd.cpp init.cpp inv.cpp list.cpp malloc.cpp mlfn.cpp msfn.cpp proj_mdist.cpp open_lib.cpp param.cpp phi2.cpp pr_list.cpp qsfn.cpp strerrno.cpp tsfn.cpp units.cpp ctx.cpp log.cpp zpoly1.cpp rtodms.cpp - vector1.cpp release.cpp gauss.cpp + release.cpp gauss.cpp fileapi.cpp gc_reader.cpp gridcatalog.cpp nad_cvt.cpp nad_init.cpp nad_intr.cpp diff --git a/src/mk_cheby.cpp b/src/mk_cheby.cpp deleted file mode 100644 index 0f3b97ed..00000000 --- a/src/mk_cheby.cpp +++ /dev/null @@ -1,193 +0,0 @@ -#include "proj.h" -#include "proj_internal.h" -static void /* sum coefficients less than res */ -eval(PJ_UV **w, int nu, int nv, double res, PJ_UV *resid) { - int i, j; - double ab; - PJ_UV *s; - - resid->u = resid->v = 0.; - for (i = 0; i < nu; ++i) { - s = w[i]; - for (j = 0; j < nv; ++j) { - if ((ab = fabs(s->u)) < res) - resid->u += ab; - if ((ab = fabs(s->v)) < res) - resid->v += ab; - ++s; - } - } -} -static Tseries * /* create power series structure */ -makeT(int nru, int nrv) { - Tseries *T; - int i; - - if (!(T = (Tseries *)pj_malloc(sizeof(Tseries)))) - return nullptr; - if (!(T->cu = (struct PW_COEF *)pj_malloc(sizeof(struct PW_COEF) * nru))) { - pj_dalloc(T); - return nullptr; - } - if (!(T->cv = (struct PW_COEF *)pj_malloc(sizeof(struct PW_COEF) * nrv))) { - pj_dalloc(T->cu); - pj_dalloc(T); - return nullptr; - } - - for (i = 0; i < nru; ++i) - T->cu[i].c = nullptr; - for (i = 0; i < nrv; ++i) - T->cv[i].c = nullptr; - return T; -} -Tseries * -mk_cheby(PJ_UV a, PJ_UV b, double res, PJ_UV *resid, PJ_UV (*func)(PJ_UV), - int nu, int nv, int power) { - int j, i, nru, nrv, *ncu, *ncv; - Tseries *T = nullptr; - PJ_UV **w; - double cutres; - - if (!(w = (PJ_UV **)vector2(nu, nv, sizeof(PJ_UV)))) - return nullptr; - if (!(ncu = (int *)vector1(nu + nv, sizeof(int)))) { - freev2((void **)w, nu); - return nullptr; - } - ncv = ncu + nu; - if (!bchgen(a, b, nu, nv, w, func)) { - PJ_UV *s; - double ab, *p; - - /* analyse coefficients and adjust until residual OK */ - cutres = res; - for (i = 4; i ; --i) { - eval(w, nu, nv, cutres, resid); - if (resid->u < res && resid->v < res) - break; - cutres *= 0.5; - } - if (i <= 0) /* warn of too many tries */ - resid->u = - resid->u; - /* apply cut resolution and set pointers */ - nru = nrv = 0; - for (j = 0; j < nu; ++j) { - ncu[j] = ncv[j] = 0; /* clear column maxes */ - s = w[j]; - for (i = 0; i < nv; ++i) { - if ((ab = fabs(s->u)) < cutres) /* < resolution ? */ - s->u = 0.; /* clear coefficient */ - else - ncu[j] = i + 1; /* update column max */ - if ((ab = fabs(s->v)) < cutres) /* same for v coef's */ - s->v = 0.; - else - ncv[j] = i + 1; - ++s; - } - if (ncu[j]) nru = j + 1; /* update row max */ - if (ncv[j]) nrv = j + 1; - } - if (power) { /* convert to bivariate power series */ - if (!bch2bps(a, b, w, nu, nv)) - goto error; - /* possible change in some row counts, so readjust */ - nru = nrv = 0; - for (j = 0; j < nu; ++j) { - ncu[j] = ncv[j] = 0; /* clear column maxes */ - s = w[j]; - for (i = 0; i < nv; ++i) { - if (s->u != 0.0) - ncu[j] = i + 1; /* update column max */ - if (s->v != 0.0) - ncv[j] = i + 1; - ++s; - } - if (ncu[j]) nru = j + 1; /* update row max */ - if (ncv[j]) nrv = j + 1; - } - if ((T = makeT(nru, nrv)) != nullptr ) { - T->a = a; - T->b = b; - T->mu = nru - 1; - T->mv = nrv - 1; - T->power = 1; - for (i = 0; i < nru; ++i) /* store coefficient rows for u */ - { - if ((T->cu[i].m = ncu[i]) != 0) - { - if ((p = T->cu[i].c = - (double *)pj_malloc(sizeof(double) * ncu[i]))) - for (j = 0; j < ncu[i]; ++j) - *p++ = (w[i] + j)->u; - else - goto error; - } - } - for (i = 0; i < nrv; ++i) /* same for v */ - { - if ((T->cv[i].m = ncv[i]) != 0) - { - if ((p = T->cv[i].c = - (double *)pj_malloc(sizeof(double) * ncv[i]))) - for (j = 0; j < ncv[i]; ++j) - *p++ = (w[i] + j)->v; - else - goto error; - } - } - } - } else if ((T = makeT(nru, nrv)) != nullptr) { - /* else make returned Chebyshev coefficient structure */ - T->mu = nru - 1; /* save row degree */ - T->mv = nrv - 1; - T->a.u = a.u + b.u; /* set argument scaling */ - T->a.v = a.v + b.v; - T->b.u = 1. / (b.u - a.u); - T->b.v = 1. / (b.v - a.v); - T->power = 0; - for (i = 0; i < nru; ++i) /* store coefficient rows for u */ - { - if ((T->cu[i].m = ncu[i]) != 0) - { - if ((p = T->cu[i].c = - (double *)pj_malloc(sizeof(double) * ncu[i]))) - for (j = 0; j < ncu[i]; ++j) - *p++ = (w[i] + j)->u; - else - goto error; - } - } - for (i = 0; i < nrv; ++i) /* same for v */ - { - if ((T->cv[i].m = ncv[i]) != 0) - { - if ((p = T->cv[i].c = - (double *)pj_malloc(sizeof(double) * ncv[i]))) - for (j = 0; j < ncv[i]; ++j) - *p++ = (w[i] + j)->v; - else - goto error; - } - } - } else - goto error; - } - goto gohome; - error: - if (T) { /* pj_dalloc up possible allocations */ - for (i = 0; i <= T->mu; ++i) - if (T->cu[i].c) - pj_dalloc(T->cu[i].c); - for (i = 0; i <= T->mv; ++i) - if (T->cv[i].c) - pj_dalloc(T->cv[i].c); - pj_dalloc(T); - } - T = nullptr; - gohome: - freev2((void **) w, nu); - pj_dalloc(ncu); - return T; -} @@ -645,6 +645,74 @@ typedef enum PJ_CS_TYPE_TEMPORALMEASURE } PJ_COORDINATE_SYSTEM_TYPE; +/** \brief Structure given overall description of a CRS. + * + * This structure may grow over time, and should not be directly allocated by + * client code. + */ +typedef struct +{ + /** Authority name. */ + char* auth_name; + /** Object code. */ + char* code; + /** Object name. */ + char* name; + /** Object type. */ + PJ_TYPE type; + /** Whether the object is deprecated */ + int deprecated; + /** Whereas the west_lon_degree, south_lat_degree, east_lon_degree and + * north_lat_degree fields are valid. */ + int bbox_valid; + /** Western-most longitude of the area of use, in degrees. */ + double west_lon_degree; + /** Southern-most latitude of the area of use, in degrees. */ + double south_lat_degree; + /** Eastern-most longitude of the area of use, in degrees. */ + double east_lon_degree; + /** Northern-most latitude of the area of use, in degrees. */ + double north_lat_degree; + /** Name of the area of use. */ + char* area_name; + /** Name of the projection method for a projected CRS. Might be NULL even + *for projected CRS in some cases. */ + char* projection_method_name; +} PROJ_CRS_INFO; + +/** \brief Structure describing optional parameters for proj_get_crs_list(); + * + * This structure may grow over time, and should not be directly allocated by + * client code. + */ +typedef struct +{ + /** Array of allowed object types. Should be NULL if all types are allowed*/ + const PJ_TYPE* types; + /** Size of types. Should be 0 if all types are allowed*/ + size_t typesCount; + + /** If TRUE and bbox_valid == TRUE, then only CRS whose area of use + * entirely contains the specified bounding box will be returned. + * If FALSE and bbox_valid == TRUE, then only CRS whose area of use + * intersects the specified bounding box will be returned. + */ + int crs_area_of_use_contains_bbox; + /** To set to TRUE so that west_lon_degree, south_lat_degree, + * east_lon_degree and north_lat_degree fields are taken into account. */ + int bbox_valid; + /** Western-most longitude of the area of use, in degrees. */ + double west_lon_degree; + /** Southern-most latitude of the area of use, in degrees. */ + double south_lat_degree; + /** Eastern-most longitude of the area of use, in degrees. */ + double east_lon_degree; + /** Northern-most latitude of the area of use, in degrees. */ + double north_lat_degree; + + /** Whether deprecated objects are allowed. Default to FALSE. */ + int allow_deprecated; +} PROJ_CRS_LIST_PARAMETERS; /**@}*/ @@ -774,6 +842,19 @@ PROJ_STRING_LIST PROJ_DLL proj_get_codes_from_database(PJ_CONTEXT *ctx, PJ_TYPE type, int allow_deprecated); +PROJ_CRS_LIST_PARAMETERS PROJ_DLL *proj_get_crs_list_parameters_create(void); + +void PROJ_DLL proj_get_crs_list_parameters_destroy( + PROJ_CRS_LIST_PARAMETERS* params); + +PROJ_CRS_INFO PROJ_DLL **proj_get_crs_info_list_from_database( + PJ_CONTEXT *ctx, + const char *auth_name, + const PROJ_CRS_LIST_PARAMETERS* params, + int *out_result_count); + +void PROJ_DLL proj_crs_info_list_destroy(PROJ_CRS_INFO** list); + /* ------------------------------------------------------------------------- */ diff --git a/src/proj_internal.h b/src/proj_internal.h index 453bd654..9ffcc2b3 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -653,7 +653,7 @@ struct FACTORS { #define PJD_ERR_LAT_0_OR_ALPHA_EQ_90 -33 #define PJD_ERR_ELLIPSOID_USE_REQUIRED -34 #define PJD_ERR_INVALID_UTM_ZONE -35 -#define PJD_ERR_TCHEBY_VAL_OUT_OF_RANGE -36 +/* -36 no longer used */ #define PJD_ERR_FAILED_TO_FIND_PROJ -37 #define PJD_ERR_FAILED_TO_LOAD_GRID -38 #define PJD_ERR_INVALID_M_OR_N -39 @@ -677,8 +677,8 @@ struct FACTORS { #define PJD_ERR_TOO_MANY_INITS -57 #define PJD_ERR_INVALID_ARG -58 #define PJD_ERR_INCONSISTENT_UNIT -59 -/* NOTE: Remember to update pj_strerrno.c and transient_error in */ -/* pj_transform.c when adding new value */ +/* NOTE: Remember to update src/strerrno.cpp, src/apps/gie.cpp and transient_error in */ +/* src/transform.cpp when adding new value */ struct projFileAPI_t; @@ -854,30 +854,6 @@ COMPLEX pj_zpolyd1(COMPLEX, const COMPLEX *, int, COMPLEX *); int pj_deriv(PJ_LP, double, const PJ *, struct DERIVS *); int pj_factors(PJ_LP, const PJ *, double, struct FACTORS *); -struct PW_COEF { /* row coefficient structure */ - int m; /* number of c coefficients (=0 for none) */ - double *c; /* power coefficients */ -}; - -/* Approximation structures and procedures */ -typedef struct { /* Chebyshev or Power series structure */ - PJ_UV a, b; /* power series range for evaluation */ - /* or Chebyshev argument shift/scaling */ - struct PW_COEF *cu, *cv; - int mu, mv; /* maximum cu and cv index (+1 for count) */ - int power; /* != 0 if power series, else Chebyshev */ -} Tseries; - -Tseries PROJ_DLL *mk_cheby(PJ_UV, PJ_UV, double, PJ_UV *, PJ_UV (*)(PJ_UV), int, int, int); -PJ_UV bpseval(PJ_UV, Tseries *); -PJ_UV bcheval(PJ_UV, Tseries *); -PJ_UV biveval(PJ_UV, Tseries *); -void *vector1(int, int); -void **vector2(int, int, int); -void freev2(void **v, int nrows); -int bchgen(PJ_UV, PJ_UV, int, int, PJ_UV **, PJ_UV(*)(PJ_UV)); -int bch2bps(PJ_UV, PJ_UV, PJ_UV **, int, int); - /* nadcon related protos */ PJ_LP nad_intr(PJ_LP, struct CTABLE *); PJ_LP nad_cvt(PJ_LP, int, struct CTABLE *); diff --git a/src/projections/etmerc.cpp b/src/projections/etmerc.cpp deleted file mode 100644 index e75bc168..00000000 --- a/src/projections/etmerc.cpp +++ /dev/null @@ -1,362 +0,0 @@ -/* -** libproj -- library of cartographic projections -** -** Copyright (c) 2008 Gerald I. Evenden -*/ - -/* -** 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. -*/ - -/* The code in this file is largly based upon procedures: - * - * Written by: Knud Poder and Karsten Engsager - * - * Based on math from: R.Koenig and K.H. Weise, "Mathematische - * Grundlagen der hoeheren Geodaesie und Kartographie, - * Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951. - * - * Modified and used here by permission of Reference Networks - * Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark - * -*/ - -#define PJ_LIB__ - -#include <errno.h> - -#include "proj.h" -#include "proj_internal.h" -#include "proj_math.h" - - -namespace { // anonymous namespace -struct pj_opaque { - double Qn; /* Merid. quad., scaled to the projection */ \ - double Zb; /* Radius vector in polar coord. systems */ \ - double cgb[6]; /* Constants for Gauss -> Geo lat */ \ - double cbg[6]; /* Constants for Geo lat -> Gauss */ \ - double utg[6]; /* Constants for transv. merc. -> geo */ \ - double gtu[6]; /* Constants for geo -> transv. merc. */ -}; -} // anonymous namespace - -PROJ_HEAD(etmerc, "Extended Transverse Mercator") - "\n\tCyl, Sph\n\tlat_ts=(0)\nlat_0=(0)"; -PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") - "\n\tCyl, Sph\n\tzone= south"; - -#define PROJ_ETMERC_ORDER 6 - -#ifdef _GNU_SOURCE - inline -#endif -static double gatg(double *p1, int len_p1, double B) { - double *p; - double h = 0, h1, h2 = 0, cos_2B; - - cos_2B = 2*cos(2*B); - p = p1 + len_p1; - h1 = *--p; - while (p - p1) { - h = -h2 + cos_2B*h1 + *--p; - h2 = h1; - h1 = h; - } - return (B + h*sin(2*B)); -} - -/* Complex Clenshaw summation */ -#ifdef _GNU_SOURCE - inline -#endif -static double clenS(double *a, int size, double arg_r, double arg_i, double *R, double *I) { - double *p, r, i, hr, hr1, hr2, hi, hi1, hi2; - double sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i; - - /* arguments */ - p = a + size; -#ifdef _GNU_SOURCE - sincos(arg_r, &sin_arg_r, &cos_arg_r); -#else - sin_arg_r = sin(arg_r); - cos_arg_r = cos(arg_r); -#endif - sinh_arg_i = sinh(arg_i); - cosh_arg_i = cosh(arg_i); - r = 2*cos_arg_r*cosh_arg_i; - i = -2*sin_arg_r*sinh_arg_i; - - /* summation loop */ - hi1 = hr1 = hi = 0; - hr = *--p; - for (; a - p;) { - hr2 = hr1; - hi2 = hi1; - hr1 = hr; - hi1 = hi; - hr = -hr2 + r*hr1 - i*hi1 + *--p; - hi = -hi2 + i*hr1 + r*hi1; - } - - r = sin_arg_r*cosh_arg_i; - i = cos_arg_r*sinh_arg_i; - *R = r*hr - i*hi; - *I = r*hi + i*hr; - return *R; -} - - -/* Real Clenshaw summation */ -static double clens(double *a, int size, double arg_r) { - double *p, r, hr, hr1, hr2, cos_arg_r; - - p = a + size; - cos_arg_r = cos(arg_r); - r = 2*cos_arg_r; - - /* summation loop */ - hr1 = 0; - hr = *--p; - for (; a - p;) { - hr2 = hr1; - hr1 = hr; - hr = -hr2 + r*hr1 + *--p; - } - return sin (arg_r)*hr; -} - - - -static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ - PJ_XY xy = {0.0,0.0}; - struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); - double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; - double Cn = lp.phi, Ce = lp.lam; - - /* ell. LAT, LNG -> Gaussian LAT, LNG */ - Cn = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn); - /* Gaussian LAT, LNG -> compl. sph. LAT */ -#ifdef _GNU_SOURCE - sincos (Cn, &sin_Cn, &cos_Cn); - sincos (Ce, &sin_Ce, &cos_Ce); -#else - sin_Cn = sin (Cn); - cos_Cn = cos (Cn); - sin_Ce = sin (Ce); - cos_Ce = cos (Ce); -#endif - - Cn = atan2 (sin_Cn, cos_Ce*cos_Cn); - Ce = atan2 (sin_Ce*cos_Cn, hypot (sin_Cn, cos_Cn*cos_Ce)); - - /* compl. sph. N, E -> ell. norm. N, E */ - Ce = asinh ( tan (Ce) ); /* Replaces: Ce = log(tan(FORTPI + Ce*0.5)); */ - Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); - Ce += dCe; - if (fabs (Ce) <= 2.623395162778) { - xy.y = Q->Qn * Cn + Q->Zb; /* Northing */ - xy.x = Q->Qn * Ce; /* Easting */ - } else - xy.x = xy.y = HUGE_VAL; - return xy; -} - - - -static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ - PJ_LP lp = {0.0,0.0}; - struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); - double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; - double Cn = xy.y, Ce = xy.x; - - /* normalize N, E */ - Cn = (Cn - Q->Zb)/Q->Qn; - Ce = Ce/Q->Qn; - - if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */ - /* norm. N, E -> compl. sph. LAT, LNG */ - Cn += clenS(Q->utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); - Ce += dCe; - Ce = atan (sinh (Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */ - /* compl. sph. LAT -> Gaussian LAT, LNG */ -#ifdef _GNU_SOURCE - sincos (Cn, &sin_Cn, &cos_Cn); - sincos (Ce, &sin_Ce, &cos_Ce); -#else - sin_Cn = sin (Cn); - cos_Cn = cos (Cn); - sin_Ce = sin (Ce); - cos_Ce = cos (Ce); -#endif - Ce = atan2 (sin_Ce, cos_Ce*cos_Cn); - Cn = atan2 (sin_Cn*cos_Ce, hypot (sin_Ce, cos_Ce*cos_Cn)); - /* Gaussian LAT, LNG -> ell. LAT, LNG */ - lp.phi = gatg (Q->cgb, PROJ_ETMERC_ORDER, Cn); - lp.lam = Ce; - } - else - lp.phi = lp.lam = HUGE_VAL; - return lp; -} - - -static PJ *setup(PJ *P) { /* general initialization */ - double f, n, np, Z; - struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); - - if (P->es <= 0) { - return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); - } - - /* flattening */ - f = P->es / (1 + sqrt (1 - P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */ - - /* third flattening */ - np = n = f/(2 - f); - - /* COEF. OF TRIG SERIES GEO <-> GAUSS */ - /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */ - /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */ - /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */ - - Q->cgb[0] = n*( 2 + n*(-2/3.0 + n*(-2 + n*(116/45.0 + n*(26/45.0 + - n*(-2854/675.0 )))))); - Q->cbg[0] = n*(-2 + n*( 2/3.0 + n*( 4/3.0 + n*(-82/45.0 + n*(32/45.0 + - n*( 4642/4725.0)))))); - np *= n; - Q->cgb[1] = np*(7/3.0 + n*( -8/5.0 + n*(-227/45.0 + n*(2704/315.0 + - n*( 2323/945.0))))); - Q->cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0 + n*( 904/315.0 + - n*(-1522/945.0))))); - np *= n; - /* n^5 coeff corrected from 1262/105 -> -1262/105 */ - Q->cgb[2] = np*( 56/15.0 + n*(-136/35.0 + n*(-1262/105.0 + - n*( 73814/2835.0)))); - Q->cbg[2] = np*(-26/15.0 + n*( 34/21.0 + n*( 8/5.0 + - n*(-12686/2835.0)))); - np *= n; - /* n^5 coeff corrected from 322/35 -> 332/35 */ - Q->cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0))); - Q->cbg[3] = np*(1237/630.0 + n*( -12/5.0 + n*( -24832/14175.0))); - np *= n; - Q->cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 )); - Q->cbg[4] = np*(-734/315.0 + n*( 109598/31185.0)); - np *= n; - Q->cgb[5] = np*(601676/22275.0 ); - Q->cbg[5] = np*(444337/155925.0); - - /* Constants of the projections */ - /* Transverse Mercator (UTM, ITM, etc) */ - np = n*n; - /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */ - Q->Qn = P->k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0))); - /* coef of trig series */ - /* utg := ell. N, E -> sph. N, E, KW p194 (65) */ - /* gtu := sph. N, E -> ell. N, E, KW p196 (69) */ - Q->utg[0] = n*(-0.5 + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 + - n*( 81/512.0 + n*(-96199/604800.0)))))); - Q->gtu[0] = n*( 0.5 + n*(-2/3.0 + n*( 5/16.0 + n*(41/180.0 + - n*(-127/288.0 + n*( 7891/37800.0 )))))); - Q->utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 + - n*( 1118711/3870720.0))))); - Q->gtu[1] = np*(13/48.0 + n*(-3/5.0 + n*(557/1440.0 + n*(281/630.0 + - n*(-1983433/1935360.0))))); - np *= n; - Q->utg[2] = np*(-17/480.0 + n*( 37/840.0 + n*( 209/4480.0 + - n*( -5569/90720.0 )))); - Q->gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 + - n*(167603/181440.0)))); - np *= n; - Q->utg[3] = np*(-4397/161280.0 + n*( 11/504.0 + n*( 830251/7257600.0))); - Q->gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0))); - np *= n; - Q->utg[4] = np*(-4583/161280.0 + n*( 108847/3991680.0)); - Q->gtu[4] = np*(34729/80640.0 + n*(-3418889/1995840.0)); - np *= n; - Q->utg[5] = np*(-20648693/638668800.0); - Q->gtu[5] = np*(212378941/319334400.0); - - /* Gaussian latitude value of the origin latitude */ - Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0); - - /* Origin northing minus true northing at the origin latitude */ - /* i.e. true northing = N - P->Zb */ - Q->Zb = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z)); - P->inv = e_inverse; - P->fwd = e_forward; - return P; -} - - - -PJ *PROJECTION(etmerc) { - struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = Q; - return setup (P); -} - - - -/* utm uses etmerc for the underlying projection */ - - -PJ *PROJECTION(utm) { - long zone; - struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = Q; - - if (P->es == 0.0) { - proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); - return pj_default_destructor(P, ENOMEM); - } - if (P->lam0 < -1000.0 || P->lam0 > 1000.0) { - return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE); - } - - P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.; - P->x0 = 500000.; - if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */ - { - zone = pj_param(P->ctx, P->params, "izone").i; - if (zone > 0 && zone <= 60) - --zone; - else { - return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE); - } - } - else /* nearest central meridian input */ - { - zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI))); - if (zone < 0) - zone = 0; - else if (zone >= 60) - zone = 59; - } - P->lam0 = (zone + .5) * M_PI / 30. - M_PI; - P->k0 = 0.9996; - P->phi0 = 0.; - - return setup (P); -} diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp index d1938116..c91c5174 100644 --- a/src/projections/tmerc.cpp +++ b/src/projections/tmerc.cpp @@ -1,3 +1,15 @@ +/* +* Transverse Mercator implementations +* +* In this file two transverse mercator implementations are found. One of Gerald +* Evenden/John Snyder origin and one of Knud Poder/Karsten Engsager origin. The +* former is regarded as "approximate" in the following and the latter is "exact". +* This word choice has been made to distinguish between the two algorithms, where +* the Evenden/Snyder implementation is the faster, less accurate implementation +* and the Poder/Engsager algorithm is a slightly slower, but more accurate +* implementation. +*/ + #define PJ_LIB__ #include <errno.h> @@ -5,18 +17,32 @@ #include "proj.h" #include "proj_internal.h" +#include "proj_math.h" -PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell"; +PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell\n\tapprox"; +PROJ_HEAD(etmerc, "Extended Transverse Mercator") "\n\tCyl, Sph"; +PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") "\n\tCyl, Sph\n\tzone= south approx"; namespace { // anonymous namespace -struct pj_opaque { +struct pj_opaque_approx { double esp; double ml0; double *en; }; + +struct pj_opaque_exact { + double Qn; /* Merid. quad., scaled to the projection */ + double Zb; /* Radius vector in polar coord. systems */ + double cgb[6]; /* Constants for Gauss -> Geo lat */ + double cbg[6]; /* Constants for Geo lat -> Gauss */ + double utg[6]; /* Constants for transv. merc. -> geo */ + double gtu[6]; /* Constants for geo -> transv. merc. */ +}; + } // anonymous namespace +/* Constants for "approximate" transverse mercator */ #define EPS10 1.e-10 #define FC1 1. #define FC2 .5 @@ -27,10 +53,18 @@ struct pj_opaque { #define FC7 .02380952380952380952 #define FC8 .01785714285714285714 +/* Constant for "exact" transverse mercator */ +#define PROJ_ETMERC_ORDER 6 -static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ +/*****************************************************************************/ +// +// Approximate Transverse Mercator functions +// +/*****************************************************************************/ + +static PJ_XY approx_e_fwd (PJ_LP lp, PJ *P) { PJ_XY xy = {0.0, 0.0}; - struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque); double al, als, n, cosphi, sinphi, t; /* @@ -70,7 +104,7 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */ } -static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ +static PJ_XY approx_s_fwd (PJ_LP lp, PJ *P) { PJ_XY xy = {0.0,0.0}; double b, cosphi; @@ -95,7 +129,7 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ return xy; } - xy.x = static_cast<struct pj_opaque*>(P->opaque)->ml0 * log ((1. + b) / (1. - b)); + xy.x = static_cast<struct pj_opaque_approx*>(P->opaque)->ml0 * log ((1. + b) / (1. - b)); xy.y = cosphi * cos (lp.lam) / sqrt (1. - b * b); b = fabs ( xy.y ); @@ -110,14 +144,14 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) { /* Spheroidal, forward */ if (lp.phi < 0.) xy.y = -xy.y; - xy.y = static_cast<struct pj_opaque*>(P->opaque)->esp * (xy.y - P->phi0); + xy.y = static_cast<struct pj_opaque_approx*>(P->opaque)->esp * (xy.y - P->phi0); return xy; } -static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ +static PJ_LP approx_e_inv (PJ_XY xy, PJ *P) { PJ_LP lp = {0.0,0.0}; - struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); + struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque); double n, con, cosphi, d, ds, sinphi, t; lp.phi = pj_inv_mlfn(P->ctx, Q->ml0 + xy.y / P->k0, P->es, Q->en); @@ -149,13 +183,13 @@ static PJ_LP e_inverse (PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */ } -static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ +static PJ_LP approx_s_inv (PJ_XY xy, PJ *P) { PJ_LP lp = {0.0, 0.0}; double h, g; - h = exp(xy.x / static_cast<struct pj_opaque*>(P->opaque)->esp); + h = exp(xy.x / static_cast<struct pj_opaque_approx*>(P->opaque)->esp); g = .5 * (h - 1. / h); - h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque*>(P->opaque)->esp); + h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque_approx*>(P->opaque)->esp); lp.phi = asin(sqrt((1. - h * h) / (1. + g * g))); /* Make sure that phi is on the correct hemisphere when false northing is used */ @@ -166,45 +200,386 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) { /* Spheroidal, inverse */ } -static PJ *destructor(PJ *P, int errlev) { /* Destructor */ +static PJ *destructor_approx(PJ *P, int errlev) { if (nullptr==P) return nullptr; if (nullptr==P->opaque) return pj_default_destructor(P, errlev); - pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->en); + pj_dealloc (static_cast<struct pj_opaque_approx*>(P->opaque)->en); return pj_default_destructor(P, errlev); } -static PJ *setup(PJ *P) { /* general initialization */ - struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque); +static PJ *setup_approx(PJ *P) { + struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(P->opaque); + + P->destructor = destructor_approx; + if (P->es != 0.0) { if (!(Q->en = pj_enfn(P->es))) return pj_default_destructor(P, ENOMEM); Q->ml0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), Q->en); Q->esp = P->es / (1. - P->es); - P->inv = e_inverse; - P->fwd = e_forward; + P->inv = approx_e_inv; + P->fwd = approx_e_fwd; } else { Q->esp = P->k0; Q->ml0 = .5 * Q->esp; - P->inv = s_inverse; - P->fwd = s_forward; + P->inv = approx_s_inv; + P->fwd = approx_s_fwd; } return P; } + +/*****************************************************************************/ +// +// Exact Transverse Mercator functions +// +// +// The code in this file is largly based upon procedures: +// +// Written by: Knud Poder and Karsten Engsager +// +// Based on math from: R.Koenig and K.H. Weise, "Mathematische +// Grundlagen der hoeheren Geodaesie und Kartographie, +// Springer-Verlag, Berlin/Goettingen" Heidelberg, 1951. +// +// Modified and used here by permission of Reference Networks +// Division, Kort og Matrikelstyrelsen (KMS), Copenhagen, Denmark +// +/*****************************************************************************/ + +/* Helper functios for "exact" transverse mercator */ +#ifdef _GNU_SOURCE + inline +#endif +static double gatg(double *p1, int len_p1, double B) { + double *p; + double h = 0, h1, h2 = 0, cos_2B; + + cos_2B = 2*cos(2*B); + p = p1 + len_p1; + h1 = *--p; + while (p - p1) { + h = -h2 + cos_2B*h1 + *--p; + h2 = h1; + h1 = h; + } + return (B + h*sin(2*B)); +} + +/* Complex Clenshaw summation */ +#ifdef _GNU_SOURCE + inline +#endif +static double clenS(double *a, int size, double arg_r, double arg_i, double *R, double *I) { + double *p, r, i, hr, hr1, hr2, hi, hi1, hi2; + double sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i; + + /* arguments */ + p = a + size; +#ifdef _GNU_SOURCE + sincos(arg_r, &sin_arg_r, &cos_arg_r); +#else + sin_arg_r = sin(arg_r); + cos_arg_r = cos(arg_r); +#endif + sinh_arg_i = sinh(arg_i); + cosh_arg_i = cosh(arg_i); + r = 2*cos_arg_r*cosh_arg_i; + i = -2*sin_arg_r*sinh_arg_i; + + /* summation loop */ + hi1 = hr1 = hi = 0; + hr = *--p; + for (; a - p;) { + hr2 = hr1; + hi2 = hi1; + hr1 = hr; + hi1 = hi; + hr = -hr2 + r*hr1 - i*hi1 + *--p; + hi = -hi2 + i*hr1 + r*hi1; + } + + r = sin_arg_r*cosh_arg_i; + i = cos_arg_r*sinh_arg_i; + *R = r*hr - i*hi; + *I = r*hi + i*hr; + return *R; +} + + +/* Real Clenshaw summation */ +static double clens(double *a, int size, double arg_r) { + double *p, r, hr, hr1, hr2, cos_arg_r; + + p = a + size; + cos_arg_r = cos(arg_r); + r = 2*cos_arg_r; + + /* summation loop */ + hr1 = 0; + hr = *--p; + for (; a - p;) { + hr2 = hr1; + hr1 = hr; + hr = -hr2 + r*hr1 + *--p; + } + return sin (arg_r)*hr; +} + +/* Ellipsoidal, forward */ +static PJ_XY exact_e_fwd (PJ_LP lp, PJ *P) { + PJ_XY xy = {0.0,0.0}; + struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque); + double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; + double Cn = lp.phi, Ce = lp.lam; + + /* ell. LAT, LNG -> Gaussian LAT, LNG */ + Cn = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn); + /* Gaussian LAT, LNG -> compl. sph. LAT */ +#ifdef _GNU_SOURCE + sincos (Cn, &sin_Cn, &cos_Cn); + sincos (Ce, &sin_Ce, &cos_Ce); +#else + sin_Cn = sin (Cn); + cos_Cn = cos (Cn); + sin_Ce = sin (Ce); + cos_Ce = cos (Ce); +#endif + + Cn = atan2 (sin_Cn, cos_Ce*cos_Cn); + Ce = atan2 (sin_Ce*cos_Cn, hypot (sin_Cn, cos_Cn*cos_Ce)); + + /* compl. sph. N, E -> ell. norm. N, E */ + Ce = asinh ( tan (Ce) ); /* Replaces: Ce = log(tan(FORTPI + Ce*0.5)); */ + Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); + Ce += dCe; + if (fabs (Ce) <= 2.623395162778) { + xy.y = Q->Qn * Cn + Q->Zb; /* Northing */ + xy.x = Q->Qn * Ce; /* Easting */ + } else + xy.x = xy.y = HUGE_VAL; + return xy; +} + + +/* Ellipsoidal, inverse */ +static PJ_LP exact_e_inv (PJ_XY xy, PJ *P) { + PJ_LP lp = {0.0,0.0}; + struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque); + double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe; + double Cn = xy.y, Ce = xy.x; + + /* normalize N, E */ + Cn = (Cn - Q->Zb)/Q->Qn; + Ce = Ce/Q->Qn; + + if (fabs(Ce) <= 2.623395162778) { /* 150 degrees */ + /* norm. N, E -> compl. sph. LAT, LNG */ + Cn += clenS(Q->utg, PROJ_ETMERC_ORDER, 2*Cn, 2*Ce, &dCn, &dCe); + Ce += dCe; + Ce = atan (sinh (Ce)); /* Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI); */ + /* compl. sph. LAT -> Gaussian LAT, LNG */ +#ifdef _GNU_SOURCE + sincos (Cn, &sin_Cn, &cos_Cn); + sincos (Ce, &sin_Ce, &cos_Ce); +#else + sin_Cn = sin (Cn); + cos_Cn = cos (Cn); + sin_Ce = sin (Ce); + cos_Ce = cos (Ce); +#endif + Ce = atan2 (sin_Ce, cos_Ce*cos_Cn); + Cn = atan2 (sin_Cn*cos_Ce, hypot (sin_Ce, cos_Ce*cos_Cn)); + /* Gaussian LAT, LNG -> ell. LAT, LNG */ + lp.phi = gatg (Q->cgb, PROJ_ETMERC_ORDER, Cn); + lp.lam = Ce; + } + else + lp.phi = lp.lam = HUGE_VAL; + return lp; +} + +static PJ *setup_exact(PJ *P) { + double f, n, np, Z; + struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque); + + if (P->es <= 0) { + return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + } + + /* flattening */ + f = P->es / (1 + sqrt (1 - P->es)); /* Replaces: f = 1 - sqrt(1-P->es); */ + + /* third flattening */ + np = n = f/(2 - f); + + /* COEF. OF TRIG SERIES GEO <-> GAUSS */ + /* cgb := Gaussian -> Geodetic, KW p190 - 191 (61) - (62) */ + /* cbg := Geodetic -> Gaussian, KW p186 - 187 (51) - (52) */ + /* PROJ_ETMERC_ORDER = 6th degree : Engsager and Poder: ICC2007 */ + + Q->cgb[0] = n*( 2 + n*(-2/3.0 + n*(-2 + n*(116/45.0 + n*(26/45.0 + + n*(-2854/675.0 )))))); + Q->cbg[0] = n*(-2 + n*( 2/3.0 + n*( 4/3.0 + n*(-82/45.0 + n*(32/45.0 + + n*( 4642/4725.0)))))); + np *= n; + Q->cgb[1] = np*(7/3.0 + n*( -8/5.0 + n*(-227/45.0 + n*(2704/315.0 + + n*( 2323/945.0))))); + Q->cbg[1] = np*(5/3.0 + n*(-16/15.0 + n*( -13/9.0 + n*( 904/315.0 + + n*(-1522/945.0))))); + np *= n; + /* n^5 coeff corrected from 1262/105 -> -1262/105 */ + Q->cgb[2] = np*( 56/15.0 + n*(-136/35.0 + n*(-1262/105.0 + + n*( 73814/2835.0)))); + Q->cbg[2] = np*(-26/15.0 + n*( 34/21.0 + n*( 8/5.0 + + n*(-12686/2835.0)))); + np *= n; + /* n^5 coeff corrected from 322/35 -> 332/35 */ + Q->cgb[3] = np*(4279/630.0 + n*(-332/35.0 + n*(-399572/14175.0))); + Q->cbg[3] = np*(1237/630.0 + n*( -12/5.0 + n*( -24832/14175.0))); + np *= n; + Q->cgb[4] = np*(4174/315.0 + n*(-144838/6237.0 )); + Q->cbg[4] = np*(-734/315.0 + n*( 109598/31185.0)); + np *= n; + Q->cgb[5] = np*(601676/22275.0 ); + Q->cbg[5] = np*(444337/155925.0); + + /* Constants of the projections */ + /* Transverse Mercator (UTM, ITM, etc) */ + np = n*n; + /* Norm. mer. quad, K&W p.50 (96), p.19 (38b), p.5 (2) */ + Q->Qn = P->k0/(1 + n) * (1 + np*(1/4.0 + np*(1/64.0 + np/256.0))); + /* coef of trig series */ + /* utg := ell. N, E -> sph. N, E, KW p194 (65) */ + /* gtu := sph. N, E -> ell. N, E, KW p196 (69) */ + Q->utg[0] = n*(-0.5 + n*( 2/3.0 + n*(-37/96.0 + n*( 1/360.0 + + n*( 81/512.0 + n*(-96199/604800.0)))))); + Q->gtu[0] = n*( 0.5 + n*(-2/3.0 + n*( 5/16.0 + n*(41/180.0 + + n*(-127/288.0 + n*( 7891/37800.0 )))))); + Q->utg[1] = np*(-1/48.0 + n*(-1/15.0 + n*(437/1440.0 + n*(-46/105.0 + + n*( 1118711/3870720.0))))); + Q->gtu[1] = np*(13/48.0 + n*(-3/5.0 + n*(557/1440.0 + n*(281/630.0 + + n*(-1983433/1935360.0))))); + np *= n; + Q->utg[2] = np*(-17/480.0 + n*( 37/840.0 + n*( 209/4480.0 + + n*( -5569/90720.0 )))); + Q->gtu[2] = np*( 61/240.0 + n*(-103/140.0 + n*(15061/26880.0 + + n*(167603/181440.0)))); + np *= n; + Q->utg[3] = np*(-4397/161280.0 + n*( 11/504.0 + n*( 830251/7257600.0))); + Q->gtu[3] = np*(49561/161280.0 + n*(-179/168.0 + n*(6601661/7257600.0))); + np *= n; + Q->utg[4] = np*(-4583/161280.0 + n*( 108847/3991680.0)); + Q->gtu[4] = np*(34729/80640.0 + n*(-3418889/1995840.0)); + np *= n; + Q->utg[5] = np*(-20648693/638668800.0); + Q->gtu[5] = np*(212378941/319334400.0); + + /* Gaussian latitude value of the origin latitude */ + Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0); + + /* Origin northing minus true northing at the origin latitude */ + /* i.e. true northing = N - P->Zb */ + Q->Zb = - Q->Qn*(Z + clens(Q->gtu, PROJ_ETMERC_ORDER, 2*Z)); + P->inv = exact_e_inv; + P->fwd = exact_e_fwd; + return P; +} + + + + +/*****************************************************************************/ +// +// Operation Setups +// +/*****************************************************************************/ + PJ *PROJECTION(tmerc) { - struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque))); + /* exact transverse mercator only exists in ellipsoidal form, */ + /* use approximate version if +a sphere is requested */ + if (pj_param (P->ctx, P->params, "bapprox").i || P->es <= 0) { + struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(pj_calloc (1, sizeof (struct pj_opaque_approx))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + + P->opaque = Q; + + return setup_approx(P); + } else { + struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + return setup_exact (P); + } +} + + +PJ *PROJECTION(etmerc) { + struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact))); if (nullptr==Q) return pj_default_destructor (P, ENOMEM); - P->opaque = Q; - P->destructor = destructor; + return setup_exact (P); +} + - return setup(P); +/* UTM uses the Poder/Engsager implementation for the underlying projection */ +/* UNLESS +approx is set in which case the Evenden/Snyder implemenation is used. */ +PJ *PROJECTION(utm) { + long zone; + if (P->es == 0.0) { + proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED); + return pj_default_destructor(P, ENOMEM); + } + if (P->lam0 < -1000.0 || P->lam0 > 1000.0) { + return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE); + } + + P->y0 = pj_param (P->ctx, P->params, "bsouth").i ? 10000000. : 0.; + P->x0 = 500000.; + if (pj_param (P->ctx, P->params, "tzone").i) /* zone input ? */ + { + zone = pj_param(P->ctx, P->params, "izone").i; + if (zone > 0 && zone <= 60) + --zone; + else { + return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE); + } + } + else /* nearest central meridian input */ + { + zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI))); + if (zone < 0) + zone = 0; + else if (zone >= 60) + zone = 59; + } + P->lam0 = (zone + .5) * M_PI / 30. - M_PI; + P->k0 = 0.9996; + P->phi0 = 0.; + + if (pj_param(P->ctx, P->params, "bapprox").i) { + struct pj_opaque_approx *Q = static_cast<struct pj_opaque_approx*>(pj_calloc (1, sizeof (struct pj_opaque_approx))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + return setup_approx(P); + } else { + struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(pj_calloc (1, sizeof (struct pj_opaque_exact))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + return setup_exact(P); + } } diff --git a/src/strerrno.cpp b/src/strerrno.cpp index 9f690041..13e9c757 100644 --- a/src/strerrno.cpp +++ b/src/strerrno.cpp @@ -44,7 +44,7 @@ pj_err_list[] = { "lat_1=lat_2 or lat_1=0 or lat_2=90", /* -33 */ "elliptical usage required", /* -34 */ "invalid UTM zone number", /* -35 */ - "arg(s) out of range for Tcheby eval", /* -36 */ + "", /* no longer used */ /* -36 */ "failed to find projection to be rotated", /* -37 */ "failed to load datum shift file", /* -38 */ "both n & m must be spec'd and > 0", /* -39 */ diff --git a/src/vector1.cpp b/src/vector1.cpp deleted file mode 100644 index fc69f5c3..00000000 --- a/src/vector1.cpp +++ /dev/null @@ -1,30 +0,0 @@ -/* make storage for one and two dimensional matricies */ -#include <stdlib.h> -#include "proj.h" -#include "proj_internal.h" - void * /* one dimension array */ -vector1(int nvals, int size) { return((void *)pj_malloc(size * nvals)); } - void /* free 2D array */ -freev2(void **v, int nrows) { - if (v) { - for (v += nrows; nrows > 0; --nrows) - pj_dalloc(*--v); - pj_dalloc(v); - } -} - void ** /* two dimension array */ -vector2(int nrows, int ncols, int size) { - void **s; - - if ((s = (void **)pj_malloc(sizeof(void *) * nrows)) != nullptr) { - int rsize, i; - - rsize = size * ncols; - for (i = 0; i < nrows; ++i) - if (!(s[i] = pj_malloc(rsize))) { - freev2(s, i); - return (void **)nullptr; - } - } - return s; -} diff --git a/test/cli/CMakeLists.txt b/test/cli/CMakeLists.txt index d197b2aa..0c4ccf1b 100644 --- a/test/cli/CMakeLists.txt +++ b/test/cli/CMakeLists.txt @@ -3,10 +3,13 @@ # set(CS2CS_BIN "cs2cs") set(PROJ_BIN "proj") +set(PROJINFO_BIN "projinfo") +set(CCT_BIN "cct") proj_add_test_script_sh("test27" PROJ_BIN ) proj_add_test_script_sh("test83" PROJ_BIN ) proj_add_test_script_sh("testvarious" CS2CS_BIN ) proj_add_test_script_sh("testdatumfile" CS2CS_BIN "connu") proj_add_test_script_sh("testIGNF" CS2CS_BIN "ntf_r93.gsb") proj_add_test_script_sh("testntv2" CS2CS_BIN "ntv2_0.gsb") - +proj_add_test_script_sh("testprojinfo" PROJINFO_BIN ) +proj_add_test_script_sh("testcct" CCT_BIN ) diff --git a/test/cli/Makefile.am b/test/cli/Makefile.am index 50654968..c0bc0871 100644 --- a/test/cli/Makefile.am +++ b/test/cli/Makefile.am @@ -5,6 +5,7 @@ EXEPATH = ../../src PROJEXE = $(EXEPATH)/proj CS2CSEXE = $(EXEPATH)/cs2cs PROJINFOEXE = $(EXEPATH)/projinfo +CCTEXE = $(EXEPATH)/cct # PROJ.4 test scripts TEST27 = $(THIS_DIR)/test27 @@ -15,12 +16,14 @@ TESTFLAKY = $(THIS_DIR)/testflaky TESTDATUMFILE = $(THIS_DIR)/testdatumfile TESTIGN = $(THIS_DIR)/testIGNF TESTPROJINFO = $(THIS_DIR)/testprojinfo +TESTCCT = $(THIS_DIR)/testcct EXTRA_DIST = pj_out27.dist pj_out83.dist td_out.dist \ test27 test83 tv_out.dist tf_out.dist \ testflaky testvarious testdatumfile testntv2 ntv2_out.dist \ testIGNF proj_outIGNF.dist \ testprojinfo testprojinfo_out.dist \ + testcct testcct_out.dist \ CMakeLists.txt testprojinfo-check: @@ -50,4 +53,7 @@ testntv2-check: PROJ_LIB=$(DATAPATH) $(TESTNTV2) $(CS2CSEXE) ; \ fi -check-local: testprojinfo-check test27-check test83-check testvarious-check testdatumfile-check testign-check testntv2-check +testcct-check: + PROJ_LIB=$(DATAPATH) $(TESTCCT) $(CCTEXE) + +check-local: testprojinfo-check test27-check test83-check testvarious-check testdatumfile-check testign-check testntv2-check testcct-check diff --git a/test/cli/testcct b/test/cli/testcct new file mode 100755 index 00000000..93749052 --- /dev/null +++ b/test/cli/testcct @@ -0,0 +1,51 @@ +: +# Test cct + +TEST_CLI_DIR=`dirname $0` +DATA_DIR=`dirname $0`/../../data +EXE=$1 + +usage() +{ + echo "Usage: ${0} <path to 'cct' program>" + echo + exit 1 +} + +if test -z "${EXE}"; then + EXE=../../src/cct +fi + +if test ! -x ${EXE}; then + echo "*** ERROR: Can not find '${EXE}' program!" + exit 1 +fi + +echo "============================================" +echo "Running ${0} using ${EXE}:" +echo "============================================" + +OUT=testcct_out + +rm -f ${OUT} + +echo "Testing cct -d 8 +proj=merc +R=1" >> ${OUT} +echo "90 45" 0 | $EXE -d 8 +proj=merc +R=1 >>${OUT} +echo "" >>${OUT} + +# do 'diff' with distribution results +echo "diff ${OUT} with testcct_out.dist" +diff -u ${OUT} ${TEST_CLI_DIR}/testcct_out.dist +if [ $? -ne 0 ] ; then + echo "" + echo "PROBLEMS HAVE OCCURRED" + echo "test file ${OUT} saved" + echo + exit 100 +else + echo "TEST OK" + echo "test file ${OUT} removed" + echo + /bin/rm -f ${OUT} + exit 0 +fi diff --git a/test/cli/testcct_out.dist b/test/cli/testcct_out.dist new file mode 100644 index 00000000..44dd6964 --- /dev/null +++ b/test/cli/testcct_out.dist @@ -0,0 +1,3 @@ +Testing cct -d 8 +proj=merc +R=1 + 1.57079633 0.88137359 0.00000000 inf + diff --git a/test/cli/testvarious b/test/cli/testvarious index 7ec50bb3..c1fa61df 100755 --- a/test/cli/testvarious +++ b/test/cli/testvarious @@ -239,7 +239,7 @@ EOF echo "##############################################################" >> ${OUT} echo "Test transverse mercator (#97)" >> ${OUT} # -$EXE +proj=tmerc +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \ +$EXE +proj=tmerc +approx +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \ +to +proj=latlong +datum=WGS84 \ -E >>${OUT} <<EOF 10000 20000 @@ -253,7 +253,7 @@ echo "##############################################################" >> ${OUT} echo "Test transverse mercator inverse (#97)" >> ${OUT} # $EXE +proj=latlong +datum=WGS84 \ - +to +proj=tmerc +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \ + +to +proj=tmerc +approx +k=0.998 +lon_0=-20 +datum=WGS84 +x_0=10000 +y_0=20000 \ -E >>${OUT} <<EOF 0dN 0.000 15d22'16.108"W 17d52'53.478"N 0.000 diff --git a/test/cli/tv_out.dist b/test/cli/tv_out.dist index 148d413d..72e95634 100644 --- a/test/cli/tv_out.dist +++ b/test/cli/tv_out.dist @@ -352,7 +352,7 @@ Test inverse handling 10 20 -1384841.19 7581707.88 0.00 ############################################################## Test MGI datum gives expected results (#207) -16.33 48.20 595710.3732102 5357598.4645755 -44.4951085 +16.33 48.20 595710.3731286 5357598.4645652 -44.4951085 ############################################################## Test omerc sensitivity with locations 90d from origin(#114) 56.958381652832 72.8798 -9985.16336453 -227.67701050 0.00000000 diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie index 82222210..930ca71f 100644 --- a/test/gie/builtins.gie +++ b/test/gie/builtins.gie @@ -5218,6 +5218,16 @@ accept -200 -100 expect -7.490535682 -0.000901935 +operation +proj=utm +zone=32 +tolerance 0.001 mm +accept 12 56 0 2000 +expect 687071.43910944 6210141.32674801 0.00000000 2000.0000 + +operation +proj=utm +zone=32 +approx +tolerance 0.001 mm +accept 12 56 0 2000 +expect 687071.43911000 6210141.32675053 0.00000000 2000.0000 + =============================================================================== van der Grinten (I) Misc Sph diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp index ff29185f..d535e412 100644 --- a/test/unit/test_c_api.cpp +++ b/test/unit/test_c_api.cpp @@ -507,36 +507,21 @@ TEST_F(CApi, proj_as_proj_string_incompatible_WKT1) { // --------------------------------------------------------------------------- -TEST_F(CApi, proj_as_proj_string_etmerc_option_yes) { +TEST_F(CApi, proj_as_proj_string_approx_tmerc_option_yes) { auto obj = proj_create(m_ctxt, "+proj=tmerc +type=crs"); ObjectKeeper keeper(obj); ASSERT_NE(obj, nullptr); - const char *options[] = {"USE_ETMERC=YES", nullptr}; + const char *options[] = {"USE_APPROX_TMERC=YES", nullptr}; auto str = proj_as_proj_string(m_ctxt, obj, PJ_PROJ_4, options); ASSERT_NE(str, nullptr); EXPECT_EQ(str, - std::string("+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 " + std::string("+proj=tmerc +approx +lat_0=0 +lon_0=0 +k=1 +x_0=0 " "+y_0=0 +datum=WGS84 +units=m +no_defs +type=crs")); } // --------------------------------------------------------------------------- -TEST_F(CApi, proj_as_proj_string_etmerc_option_no) { - auto obj = proj_create(m_ctxt, "+proj=utm +zone=31 +type=crs"); - ObjectKeeper keeper(obj); - ASSERT_NE(obj, nullptr); - - const char *options[] = {"USE_ETMERC=NO", nullptr}; - auto str = proj_as_proj_string(m_ctxt, obj, PJ_PROJ_4, options); - ASSERT_NE(str, nullptr); - EXPECT_EQ(str, std::string("+proj=tmerc +lat_0=0 +lon_0=3 +k=0.9996 " - "+x_0=500000 +y_0=0 +datum=WGS84 +units=m " - "+no_defs +type=crs")); -} - -// --------------------------------------------------------------------------- - TEST_F(CApi, proj_crs_create_bound_crs_to_WGS84) { auto crs = proj_create_from_database(m_ctxt, "EPSG", "3844", PJ_CATEGORY_CRS, false, nullptr); @@ -3022,4 +3007,210 @@ TEST_F(CApi, proj_create_cartesian_2D_cs) { } } +// --------------------------------------------------------------------------- + +TEST_F(CApi, proj_get_crs_info_list_from_database) { + { proj_crs_info_list_destroy(nullptr); } + + { proj_get_crs_list_parameters_destroy(nullptr); } + + // All null parameters + { + auto list = proj_get_crs_info_list_from_database(nullptr, nullptr, + nullptr, nullptr); + ASSERT_NE(list, nullptr); + ASSERT_NE(list[0], nullptr); + EXPECT_NE(list[0]->auth_name, nullptr); + EXPECT_NE(list[0]->code, nullptr); + EXPECT_NE(list[0]->name, nullptr); + proj_crs_info_list_destroy(list); + } + + // Default parameters + { + int result_count = 0; + auto params = proj_get_crs_list_parameters_create(); + auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params, + &result_count); + proj_get_crs_list_parameters_destroy(params); + ASSERT_NE(list, nullptr); + EXPECT_GT(result_count, 1); + EXPECT_EQ(list[result_count], nullptr); + bool found4326 = false; + bool found4978 = false; + bool found4979 = false; + bool found32631 = false; + bool found3855 = false; + bool found3901 = false; + for (int i = 0; i < result_count; i++) { + auto code = std::string(list[i]->code); + if (code == "4326") { + found4326 = true; + EXPECT_EQ(std::string(list[i]->auth_name), "EPSG"); + EXPECT_EQ(std::string(list[i]->name), "WGS 84"); + EXPECT_EQ(list[i]->type, PJ_TYPE_GEOGRAPHIC_2D_CRS); + EXPECT_EQ(list[i]->deprecated, 0); + EXPECT_EQ(list[i]->bbox_valid, 1); + EXPECT_EQ(list[i]->west_lon_degree, -180.0); + EXPECT_EQ(list[i]->south_lat_degree, -90.0); + EXPECT_EQ(list[i]->east_lon_degree, 180.0); + EXPECT_EQ(list[i]->north_lat_degree, 90.0); + EXPECT_EQ(std::string(list[i]->area_name), "World"); + EXPECT_EQ(list[i]->projection_method_name, nullptr); + } else if (code == "4978") { + found4978 = true; + EXPECT_EQ(list[i]->type, PJ_TYPE_GEOCENTRIC_CRS); + } else if (code == "4979") { + found4979 = true; + EXPECT_EQ(list[i]->type, PJ_TYPE_GEOGRAPHIC_3D_CRS); + } else if (code == "32631") { + found32631 = true; + EXPECT_EQ(list[i]->type, PJ_TYPE_PROJECTED_CRS); + EXPECT_EQ(std::string(list[i]->projection_method_name), + "Transverse Mercator"); + } else if (code == "3855") { + found3855 = true; + EXPECT_EQ(list[i]->type, PJ_TYPE_VERTICAL_CRS); + } else if (code == "3901") { + found3901 = true; + EXPECT_EQ(list[i]->type, PJ_TYPE_COMPOUND_CRS); + } + EXPECT_EQ(list[i]->deprecated, 0); + } + EXPECT_TRUE(found4326); + EXPECT_TRUE(found4978); + EXPECT_TRUE(found4979); + EXPECT_TRUE(found32631); + EXPECT_TRUE(found3855); + EXPECT_TRUE(found3901); + proj_crs_info_list_destroy(list); + } + + // Filter on only geodetic crs + { + int result_count = 0; + auto params = proj_get_crs_list_parameters_create(); + params->typesCount = 1; + auto type = PJ_TYPE_GEODETIC_CRS; + params->types = &type; + auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params, + &result_count); + bool foundGeog2D = false; + bool foundGeog3D = false; + bool foundGeocentric = false; + for (int i = 0; i < result_count; i++) { + foundGeog2D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS; + foundGeog3D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS; + foundGeocentric |= list[i]->type == PJ_TYPE_GEOCENTRIC_CRS; + EXPECT_TRUE(list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS || + list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS || + list[i]->type == PJ_TYPE_GEOCENTRIC_CRS); + } + EXPECT_TRUE(foundGeog2D); + EXPECT_TRUE(foundGeog3D); + EXPECT_TRUE(foundGeocentric); + proj_get_crs_list_parameters_destroy(params); + proj_crs_info_list_destroy(list); + } + + // Filter on only geographic crs + { + int result_count = 0; + auto params = proj_get_crs_list_parameters_create(); + params->typesCount = 1; + auto type = PJ_TYPE_GEOGRAPHIC_CRS; + params->types = &type; + auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params, + &result_count); + bool foundGeog2D = false; + bool foundGeog3D = false; + for (int i = 0; i < result_count; i++) { + foundGeog2D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS; + foundGeog3D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS; + EXPECT_TRUE(list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS || + list[i]->type == PJ_TYPE_GEOGRAPHIC_3D_CRS); + } + EXPECT_TRUE(foundGeog2D); + EXPECT_TRUE(foundGeog3D); + proj_get_crs_list_parameters_destroy(params); + proj_crs_info_list_destroy(list); + } + + // Filter on only geographic 2D crs and projected CRS + { + int result_count = 0; + auto params = proj_get_crs_list_parameters_create(); + params->typesCount = 2; + const PJ_TYPE types[] = {PJ_TYPE_GEOGRAPHIC_2D_CRS, + PJ_TYPE_PROJECTED_CRS}; + params->types = types; + auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params, + &result_count); + bool foundGeog2D = false; + bool foundProjected = false; + for (int i = 0; i < result_count; i++) { + foundGeog2D |= list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS; + foundProjected |= list[i]->type == PJ_TYPE_PROJECTED_CRS; + EXPECT_TRUE(list[i]->type == PJ_TYPE_GEOGRAPHIC_2D_CRS || + list[i]->type == PJ_TYPE_PROJECTED_CRS); + } + EXPECT_TRUE(foundGeog2D); + EXPECT_TRUE(foundProjected); + proj_get_crs_list_parameters_destroy(params); + proj_crs_info_list_destroy(list); + } + + // Filter on bbox (inclusion) + { + int result_count = 0; + auto params = proj_get_crs_list_parameters_create(); + params->bbox_valid = 1; + params->west_lon_degree = 2; + params->south_lat_degree = 49; + params->east_lon_degree = 2.1; + params->north_lat_degree = 49.1; + params->typesCount = 1; + auto type = PJ_TYPE_PROJECTED_CRS; + params->types = &type; + auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params, + &result_count); + ASSERT_NE(list, nullptr); + EXPECT_GT(result_count, 1); + for (int i = 0; i < result_count; i++) { + EXPECT_LE(list[i]->west_lon_degree, params->west_lon_degree); + EXPECT_LE(list[i]->south_lat_degree, params->south_lat_degree); + EXPECT_GE(list[i]->east_lon_degree, params->east_lon_degree); + EXPECT_GE(list[i]->north_lat_degree, params->north_lat_degree); + } + proj_get_crs_list_parameters_destroy(params); + proj_crs_info_list_destroy(list); + } + + // Filter on bbox (intersection) + { + int result_count = 0; + auto params = proj_get_crs_list_parameters_create(); + params->bbox_valid = 1; + params->west_lon_degree = 2; + params->south_lat_degree = 49; + params->east_lon_degree = 2.1; + params->north_lat_degree = 49.1; + params->crs_area_of_use_contains_bbox = 0; + params->typesCount = 1; + auto type = PJ_TYPE_PROJECTED_CRS; + params->types = &type; + auto list = proj_get_crs_info_list_from_database(m_ctxt, "EPSG", params, + &result_count); + ASSERT_NE(list, nullptr); + EXPECT_GT(result_count, 1); + for (int i = 0; i < result_count; i++) { + EXPECT_LE(list[i]->west_lon_degree, params->east_lon_degree); + EXPECT_LE(list[i]->south_lat_degree, params->north_lat_degree); + EXPECT_GE(list[i]->east_lon_degree, params->west_lon_degree); + EXPECT_GE(list[i]->north_lat_degree, params->south_lat_degree); + } + proj_get_crs_list_parameters_destroy(params); + proj_crs_info_list_destroy(list); + } +} } // namespace diff --git a/test/unit/test_factory.cpp b/test/unit/test_factory.cpp index 81f745f5..ce019079 100644 --- a/test/unit/test_factory.cpp +++ b/test/unit/test_factory.cpp @@ -222,6 +222,16 @@ TEST(factory, AuthorityFactory_createExtent) { // --------------------------------------------------------------------------- +TEST(factory, AuthorityFactory_createExtent_no_bbox) { + auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto extent = factory->createExtent("1361"); // Sudan - south. Deprecated + EXPECT_EQ(*(extent->description()), "Sudan - south"); + const auto &geogElts = extent->geographicElements(); + EXPECT_TRUE(geogElts.empty()); +} + +// --------------------------------------------------------------------------- + TEST(factory, AuthorityFactory_createGeodeticDatum) { auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); EXPECT_THROW(factory->createGeodeticDatum("-1"), @@ -414,6 +424,17 @@ TEST(factory, AuthorityFactory_createGeodeticCRS_geographic2D) { // --------------------------------------------------------------------------- +TEST(factory, AuthorityFactory_createGeodeticCRS_geographic2D_area_no_bbox) { + auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); + auto crs = factory->createGeodeticCRS("4296"); // Sudan - deprecated + auto domain = crs->domains()[0]; + auto extent = domain->domainOfValidity(); + ASSERT_TRUE(extent != nullptr); + EXPECT_TRUE(extent->isEquivalentTo(factory->createExtent("1361").get())); +} + +// --------------------------------------------------------------------------- + TEST(factory, AuthorityFactory_createGeodeticCRS_geographic3D) { auto factory = AuthorityFactory::create(DatabaseContext::create(), "EPSG"); auto crs = factory->createGeodeticCRS("4979"); @@ -2789,4 +2810,103 @@ TEST(factory, listAreaOfUseFromName) { } } +// --------------------------------------------------------------------------- + +TEST(factory, getCRSInfoList) { + auto ctxt = DatabaseContext::create(); + { + auto factory = AuthorityFactory::create(ctxt, std::string()); + auto list = factory->getCRSInfoList(); + EXPECT_GT(list.size(), 1U); + bool foundEPSG = false; + bool foundIGNF = true; + bool found4326 = false; + for (const auto &info : list) { + foundEPSG |= info.authName == "EPSG"; + foundIGNF |= info.authName == "IGNF"; + if (info.authName == "EPSG" && info.code == "4326") { + found4326 = true; + } + } + EXPECT_TRUE(foundEPSG); + EXPECT_TRUE(foundIGNF); + EXPECT_TRUE(found4326); + } + { + auto factory = AuthorityFactory::create(ctxt, "EPSG"); + auto list = factory->getCRSInfoList(); + EXPECT_GT(list.size(), 1U); + bool found4326 = false; + bool found4978 = false; + bool found4979 = false; + bool found32631 = false; + bool found3855 = false; + bool found6871 = false; + for (const auto &info : list) { + EXPECT_EQ(info.authName, "EPSG"); + if (info.code == "4326") { + EXPECT_EQ(info.name, "WGS 84"); + EXPECT_EQ(info.type, + AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS); + EXPECT_EQ(info.deprecated, false); + EXPECT_EQ(info.bbox_valid, true); + EXPECT_EQ(info.west_lon_degree, -180.0); + EXPECT_EQ(info.south_lat_degree, -90.0); + EXPECT_EQ(info.east_lon_degree, 180.0); + EXPECT_EQ(info.north_lat_degree, 90.0); + EXPECT_EQ(info.areaName, "World"); + EXPECT_TRUE(info.projectionMethodName.empty()); + found4326 = true; + } else if (info.code == "4296") { // Soudan - deprecated + EXPECT_EQ(info.bbox_valid, false); + EXPECT_EQ(info.west_lon_degree, 0.0); + EXPECT_EQ(info.south_lat_degree, 0.0); + EXPECT_EQ(info.east_lon_degree, 0.0); + EXPECT_EQ(info.north_lat_degree, 0.0); + } else if (info.code == "4978") { + EXPECT_EQ(info.name, "WGS 84"); + EXPECT_EQ(info.type, + AuthorityFactory::ObjectType::GEOCENTRIC_CRS); + found4978 = true; + } else if (info.code == "4979") { + EXPECT_EQ(info.name, "WGS 84"); + EXPECT_EQ(info.type, + AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS); + found4979 = true; + } else if (info.code == "32631") { + EXPECT_EQ(info.name, "WGS 84 / UTM zone 31N"); + EXPECT_EQ(info.type, + AuthorityFactory::ObjectType::PROJECTED_CRS); + EXPECT_EQ(info.deprecated, false); + EXPECT_EQ(info.bbox_valid, true); + EXPECT_EQ(info.west_lon_degree, 0.0); + EXPECT_EQ(info.south_lat_degree, 0.0); + EXPECT_EQ(info.east_lon_degree, 6.0); + EXPECT_EQ(info.north_lat_degree, 84.0); + EXPECT_EQ(info.areaName, "World - N hemisphere - 0\xC2\xB0" + "E to 6\xC2\xB0" + "E - by country"); + EXPECT_EQ(info.projectionMethodName, "Transverse Mercator"); + found32631 = true; + } else if (info.code == "3855") { + EXPECT_EQ(info.name, "EGM2008 height"); + EXPECT_EQ(info.type, + AuthorityFactory::ObjectType::VERTICAL_CRS); + found3855 = true; + } else if (info.code == "6871") { + EXPECT_EQ(info.name, + "WGS 84 / Pseudo-Mercator + EGM2008 geoid height"); + EXPECT_EQ(info.type, + AuthorityFactory::ObjectType::COMPOUND_CRS); + found6871 = true; + } + } + EXPECT_TRUE(found4326); + EXPECT_TRUE(found4978); + EXPECT_TRUE(found4979); + EXPECT_TRUE(found32631); + EXPECT_TRUE(found3855); + EXPECT_TRUE(found6871); + } +} } // namespace diff --git a/test/unit/test_io.cpp b/test/unit/test_io.cpp index 69ef9073..3cb50352 100644 --- a/test/unit/test_io.cpp +++ b/test/unit/test_io.cpp @@ -7654,12 +7654,39 @@ TEST(io, projparse_etmerc) { crs->exportToPROJString( PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) .get()), - "+proj=etmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 " + "+proj=tmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 " "+datum=WGS84 +units=m +no_defs +type=crs"); auto wkt1 = crs->exportToWKT( WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); - EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=etmerc +lat_0=0 " + EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\"") == std::string::npos) + << wkt1; +} + +// --------------------------------------------------------------------------- + +TEST(io, projparse_tmerc_approx) { + auto obj = + PROJStringParser().createFromPROJString("+proj=tmerc +approx +type=crs"); + auto crs = nn_dynamic_pointer_cast<ProjectedCRS>(obj); + ASSERT_TRUE(crs != nullptr); + auto wkt2 = crs->exportToWKT( + &WKTFormatter::create()->simulCurNodeHasId().setMultiLine(false)); + EXPECT_TRUE( + wkt2.find("METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]]") != + std::string::npos) + << wkt2; + + EXPECT_EQ( + crs->exportToPROJString( + PROJStringFormatter::create(PROJStringFormatter::Convention::PROJ_4) + .get()), + "+proj=tmerc +approx +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 " + "+datum=WGS84 +units=m +no_defs +type=crs"); + + auto wkt1 = crs->exportToWKT( + WKTFormatter::create(WKTFormatter::Convention::WKT1_GDAL).get()); + EXPECT_TRUE(wkt1.find("EXTENSION[\"PROJ4\",\"+proj=tmerc +approx +lat_0=0 " "+lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m " "+no_defs\"]") != std::string::npos) << wkt1; diff --git a/test/unit/test_operation.cpp b/test/unit/test_operation.cpp index 6c1ecf1c..586226e2 100644 --- a/test/unit/test_operation.cpp +++ b/test/unit/test_operation.cpp @@ -1237,9 +1237,9 @@ TEST(operation, tmerc_export) { { auto formatter = PROJStringFormatter::create(); - formatter->setUseETMercForTMerc(true); + formatter->setUseApproxTMerc(true); EXPECT_EQ(conv->exportToPROJString(formatter.get()), - "+proj=etmerc +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5"); + "+proj=tmerc +approx +lat_0=1 +lon_0=2 +k=3 +x_0=4 +y_0=5"); } EXPECT_EQ(conv->exportToWKT(WKTFormatter::create().get()), @@ -6061,21 +6061,21 @@ TEST(operation, compoundCRS_to_compoundCRS_with_vertical_transform) { "+ellps=WGS84"); { auto formatter = PROJStringFormatter::create(); - formatter->setUseETMercForTMerc(true); + formatter->setUseApproxTMerc(true); EXPECT_EQ(op->exportToPROJString(formatter.get()), - "+proj=pipeline +step +inv +proj=etmerc +lat_0=1 +lon_0=2 " + "+proj=pipeline +step +inv +proj=tmerc +approx +lat_0=1 +lon_0=2 " "+k=3 +x_0=4 +y_0=5 +ellps=WGS84 +step " "+proj=vgridshift +grids=bla.gtx +multiplier=0.001 +step " - "+proj=utm +zone=32 " + "+proj=utm +approx +zone=32 " "+ellps=WGS84"); } { auto formatter = PROJStringFormatter::create(); - formatter->setUseETMercForTMerc(true); + formatter->setUseApproxTMerc(true); EXPECT_EQ(op->inverse()->exportToPROJString(formatter.get()), - "+proj=pipeline +step +inv +proj=utm +zone=32 +ellps=WGS84 " + "+proj=pipeline +step +inv +proj=utm +approx +zone=32 +ellps=WGS84 " "+step +inv +proj=vgridshift +grids=bla.gtx " - "+multiplier=0.001 +step +proj=etmerc +lat_0=1 +lon_0=2 " + "+multiplier=0.001 +step +proj=tmerc +approx +lat_0=1 +lon_0=2 " "+k=3 +x_0=4 +y_0=5 +ellps=WGS84"); } |
