aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am7
-rw-r--r--src/apps/cct.cpp2
-rw-r--r--src/apps/gie.cpp6
-rw-r--r--src/bch2bps.cpp154
-rw-r--r--src/bchgen.cpp59
-rw-r--r--src/biveval.cpp88
-rw-r--r--src/iso19111/c_api.cpp231
-rw-r--r--src/iso19111/coordinateoperation.cpp27
-rw-r--r--src/iso19111/factory.cpp167
-rw-r--r--src/iso19111/io.cpp22
-rw-r--r--src/lib_proj.cmake7
-rw-r--r--src/mk_cheby.cpp193
-rw-r--r--src/proj.h81
-rw-r--r--src/proj_internal.h30
-rw-r--r--src/projections/etmerc.cpp362
-rw-r--r--src/projections/tmerc.cpp423
-rw-r--r--src/strerrno.cpp2
-rw-r--r--src/vector1.cpp30
18 files changed, 887 insertions, 1004 deletions
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;
-}
diff --git a/src/proj.h b/src/proj.h
index 552f4b26..f25c3228 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -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;
-}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-29 07:39:55 +0000