From 1c60f4cc408e85aff78482659a80fe974ee5d57b Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 2 Feb 2019 09:44:44 +0100 Subject: PROJStringSyntaxParser: avoid assertion on illegal input. Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12837. Credit to OSS Fuzz --- src/iso19111/io.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 3517c225..f854e21a 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -5285,6 +5285,8 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, const char *c_str = projString.c_str(); std::vector tokens; + bool hasProj = false; + bool hasInit = false; { size_t i = 0; while (true) { @@ -5313,6 +5315,13 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, if (token.empty()) { break; } + if (!hasProj && + (starts_with(token, "proj=") || starts_with(token, "+proj="))) { + hasProj = true; + } else if (!hasInit && (starts_with(token, "init=") || + starts_with(token, "+init="))) { + hasInit = true; + } tokens.emplace_back(token); } } @@ -5320,14 +5329,6 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, bool prevWasTitle = false; if (projString.find("proj=pipeline") == std::string::npos) { - const bool hasProj = projString.find("proj=") == 0 || - projString.find("+proj=") == 0 || - projString.find(" proj=") != std::string::npos || - projString.find(" +proj=") != std::string::npos; - const bool hasInit = projString.find("init=") == 0 || - projString.find("+init=") == 0 || - projString.find(" init=") != std::string::npos || - projString.find(" +init=") != std::string::npos; if (hasProj || hasInit) { steps.push_back(Step()); } -- cgit v1.2.3 From 3ed6153400299db04163ffe8fd9fd785eb9ca12e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 2 Feb 2019 11:13:12 +0100 Subject: createFromUserInput(): fix infinite recursion. Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12842. Credit to OSS Fuzz --- src/iso19111/io.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index f854e21a..1f4a7c8b 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -5287,6 +5287,7 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, bool hasProj = false; bool hasInit = false; + bool hasPipeline = false; { size_t i = 0; while (true) { @@ -5315,8 +5316,11 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, if (token.empty()) { break; } - if (!hasProj && - (starts_with(token, "proj=") || starts_with(token, "+proj="))) { + if (!hasPipeline && + (token == "proj=pipeline" || token == "+proj=pipeline")) { + hasPipeline = true; + } else if (!hasProj && (starts_with(token, "proj=") || + starts_with(token, "+proj="))) { hasProj = true; } else if (!hasInit && (starts_with(token, "init=") || starts_with(token, "+init="))) { @@ -5328,7 +5332,7 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, bool prevWasTitle = false; - if (projString.find("proj=pipeline") == std::string::npos) { + if (!hasPipeline) { if (hasProj || hasInit) { steps.push_back(Step()); } -- cgit v1.2.3 From abad23412f0920276c32567c8f237be23aa94941 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sat, 2 Feb 2019 11:25:36 +0100 Subject: pj_ellipsoid(): avoid division by zero when R=0. Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12843. Credit to OSS Fuzz --- src/ell_set.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/ell_set.cpp b/src/ell_set.cpp index f4228be8..4c9fc892 100644 --- a/src/ell_set.cpp +++ b/src/ell_set.cpp @@ -88,7 +88,8 @@ int pj_ellipsoid (PJ *P) { /* Specifying R overrules everything */ if (pj_get_param (P->params, "R")) { - ellps_size (P); + if (0 != ellps_size (P)) + return 1; pj_calc_ellipsoid_params (P, P->a, 0); if (proj_errno (P)) return 1; -- cgit v1.2.3 From b30ed5cba50688cea9ac2af90855177f84efea5a Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Feb 2019 12:23:50 +0100 Subject: labrd: avoid floating point division by zero. Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12849. Credit to OSS Fuzz --- src/projections/labrd.cpp | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'src') diff --git a/src/projections/labrd.cpp b/src/projections/labrd.cpp index 330c105f..85ab3ddd 100644 --- a/src/projections/labrd.cpp +++ b/src/projections/labrd.cpp @@ -108,6 +108,10 @@ PJ *PROJECTION(labrd) { return pj_default_destructor (P, ENOMEM); P->opaque = Q; + if (P->phi0 == 0.) { + return pj_default_destructor(P, PJD_ERR_LAT_0_IS_ZERO); + } + Az = pj_param(P->ctx, P->params, "razi").f; sinp = sin(P->phi0); t = 1. - P->es * sinp * sinp; -- cgit v1.2.3 From 24c1272cf4397614fb7ea5502bba9258e6a8f972 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Feb 2019 14:48:17 +0100 Subject: WKT1 export: avoid division by zero. Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12854. Credit to OSS Fuzz --- src/iso19111/common.cpp | 15 ++++++++------- src/iso19111/coordinateoperation.cpp | 16 ++++++++++++---- 2 files changed, 20 insertions(+), 11 deletions(-) (limited to 'src') diff --git a/src/iso19111/common.cpp b/src/iso19111/common.cpp index 57654d84..bdd836e1 100644 --- a/src/iso19111/common.cpp +++ b/src/iso19111/common.cpp @@ -172,21 +172,22 @@ void UnitOfMeasure::_exportToWKT( { const bool isWKT2 = formatter->version() == WKTFormatter::Version::WKT2; - if (formatter->forceUNITKeyword() && type() != Type::PARAMETRIC) { + const auto l_type = type(); + if (formatter->forceUNITKeyword() && l_type != Type::PARAMETRIC) { formatter->startNode(WKTConstants::UNIT, !codeSpace().empty()); } else if (!unitType.empty()) { formatter->startNode(unitType, !codeSpace().empty()); } else { - if (isWKT2 && type() == Type::LINEAR) { + if (isWKT2 && l_type == Type::LINEAR) { formatter->startNode(WKTConstants::LENGTHUNIT, !codeSpace().empty()); - } else if (isWKT2 && type() == Type::ANGULAR) { + } else if (isWKT2 && l_type == Type::ANGULAR) { formatter->startNode(WKTConstants::ANGLEUNIT, !codeSpace().empty()); - } else if (isWKT2 && type() == Type::SCALE) { + } else if (isWKT2 && l_type == Type::SCALE) { formatter->startNode(WKTConstants::SCALEUNIT, !codeSpace().empty()); - } else if (isWKT2 && type() == Type::TIME) { + } else if (isWKT2 && l_type == Type::TIME) { formatter->startNode(WKTConstants::TIMEUNIT, !codeSpace().empty()); - } else if (isWKT2 && type() == Type::PARAMETRIC) { + } else if (isWKT2 && l_type == Type::PARAMETRIC) { formatter->startNode(WKTConstants::PARAMETRICUNIT, !codeSpace().empty()); } else { @@ -211,7 +212,7 @@ void UnitOfMeasure::_exportToWKT( formatter->addQuotedString(l_name); } const auto &factor = conversionToSI(); - if (!isWKT2 || factor != 0.0) { + if (!isWKT2 || l_type != Type::TIME || factor != 0.0) { // Some TIMEUNIT do not have a conversion factor formatter->add(factor); } diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp index 6f9b6283..8a10bc5a 100644 --- a/src/iso19111/coordinateoperation.cpp +++ b/src/iso19111/coordinateoperation.cpp @@ -2147,11 +2147,19 @@ void ParameterValue::_exportToWKT(io::WKTFormatter *formatter) const { // registered linear / angular unit. const auto &unitType = unit.type(); if (unitType == common::UnitOfMeasure::Type::LINEAR) { - formatter->add( - l_value.convertToUnit(*(formatter->axisLinearUnit()))); + const auto &targetUnit = *(formatter->axisLinearUnit()); + if (targetUnit.conversionToSI() == 0.0) { + throw io::FormattingException( + "cannot convert value to target linear unit"); + } + formatter->add(l_value.convertToUnit(targetUnit)); } else if (unitType == common::UnitOfMeasure::Type::ANGULAR) { - formatter->add( - l_value.convertToUnit(*(formatter->axisAngularUnit()))); + const auto &targetUnit = *(formatter->axisAngularUnit()); + if (targetUnit.conversionToSI() == 0.0) { + throw io::FormattingException( + "cannot convert value to target angular unit"); + } + formatter->add(l_value.convertToUnit(targetUnit)); } else { formatter->add(l_value.getSIValue()); } -- cgit v1.2.3 From e121dc35abafd9a2eda01b05aa5fa13b13e6196e Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Feb 2019 15:16:01 +0100 Subject: Avoid division by zero in Ellipsoid::computeSemiMinorAxis(). Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12867. Credit to OSS Fuzz. master only --- src/iso19111/datum.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/iso19111/datum.cpp b/src/iso19111/datum.cpp index 4df1b319..66717f70 100644 --- a/src/iso19111/datum.cpp +++ b/src/iso19111/datum.cpp @@ -646,7 +646,8 @@ EllipsoidNNPtr Ellipsoid::createSphere(const util::PropertyMap &properties, * @param properties See \ref general_properties. * At minimum the name should be defined. * @param semiMajorAxisIn the semi-major axis. - * @param invFlattening the inverse/reverse flattening. + * @param invFlattening the inverse/reverse flattening. If set to 0, this will + * be considered as a sphere. * @param celestialBody Name of the celestial body on which the ellipsoid refers * to. * @return new Ellipsoid. @@ -654,8 +655,11 @@ EllipsoidNNPtr Ellipsoid::createSphere(const util::PropertyMap &properties, EllipsoidNNPtr Ellipsoid::createFlattenedSphere( const util::PropertyMap &properties, const common::Length &semiMajorAxisIn, const common::Scale &invFlattening, const std::string &celestialBody) { - auto ellipsoid(Ellipsoid::nn_make_shared( - semiMajorAxisIn, invFlattening, celestialBody)); + auto ellipsoid(invFlattening.value() == 0 + ? Ellipsoid::nn_make_shared(semiMajorAxisIn, + celestialBody) + : Ellipsoid::nn_make_shared( + semiMajorAxisIn, invFlattening, celestialBody)); ellipsoid->setProperties(properties); return ellipsoid; } -- cgit v1.2.3 From 4abfe10ab7a64552db242aa240bbcf2f92598604 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 3 Feb 2019 16:07:10 +0100 Subject: init(): repair to_meter=num/denom that was broken in the general case in PROJ 5; repair vto_meter=num/denom that was broken, and avoid division by zero, which fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12869. Credit to OSS Fuzz --- src/init.cpp | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) (limited to 'src') diff --git a/src/init.cpp b/src/init.cpp index 482e398d..2961bcca 100644 --- a/src/init.cpp +++ b/src/init.cpp @@ -764,20 +764,18 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i s = units[i].to_meter; } if (s || (s = pj_param(ctx, start, "sto_meter").s)) { - double factor; - int ratio = 0; - - /* ratio number? */ - if (strlen (s) > 1 && s[0] == '1' && s[1]=='/') { - ratio = 1; - s += 2; + char* end_ptr = const_cast(s); + PIN->to_meter = pj_strtod(s, &end_ptr); + s = end_ptr; + if (*s == '/') { /* ratio number */ + ++s; + double denom = pj_strtod(s, nullptr); + if (denom == 0.0) + return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0); + PIN->to_meter /= denom; } - - factor = pj_strtod(s, nullptr); - if ((factor <= 0.0) || (1/factor==0)) + if (PIN->to_meter <= 0.0) return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0); - - PIN->to_meter = ratio? 1 / factor: factor; PIN->fr_meter = 1 / PIN->to_meter; } else @@ -792,9 +790,16 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i s = units[i].to_meter; } if (s || (s = pj_param(ctx, start, "svto_meter").s)) { - PIN->vto_meter = pj_strtod(s, nullptr); - if (*s == '/') /* ratio number */ - PIN->vto_meter /= pj_strtod(++s, nullptr); + char* end_ptr = const_cast(s); + PIN->vto_meter = pj_strtod(s, &end_ptr); + s = end_ptr; + if (*s == '/') { /* ratio number */ + ++s; + double denom = pj_strtod(s, nullptr); + if (denom == 0.0) + return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0); + PIN->vto_meter /= denom; + } if (PIN->vto_meter <= 0.0) return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0); PIN->vfr_meter = 1. / PIN->vto_meter; -- cgit v1.2.3 From a1ea80e2b2d94bea9c1020a9651f7a6eda461767 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fischer?= Date: Sun, 3 Feb 2019 17:31:06 +0100 Subject: typo fix --- src/iso19111/factory.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index 54fc3fd3..c862cc5b 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -1733,7 +1733,7 @@ AuthorityFactory::createEllipsoid(const std::string &code) const { common::Length(c_locale_stod(semi_minor_axis_str), uom), body); } } catch (const std::exception &ex) { - throw buildFactoryException("elllipsoid", code, ex); + throw buildFactoryException("ellipsoid", code, ex); } } -- cgit v1.2.3 From 322493235696097b94294f5f466a32f015a95626 Mon Sep 17 00:00:00 2001 From: Elliott Sales de Andrade Date: Sun, 3 Feb 2019 01:33:39 -0500 Subject: Use pkgconfig to find sqlite3. --- src/Makefile.am | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index fabea4bf..e60cd3c5 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -7,7 +7,7 @@ TESTS = geodtest check_PROGRAMS = geodtest AM_CPPFLAGS = -DPROJ_LIB=\"$(pkgdatadir)\" \ - -DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ -I$(top_srcdir)/include @SQLITE3_FLAGS@ + -DMUTEX_@MUTEX_SETTING@ @JNI_INCLUDE@ -I$(top_srcdir)/include @SQLITE3_CFLAGS@ AM_CXXFLAGS = @CXX_WFLAGS@ @FLTO_FLAG@ -DPROJ_COMPILATION include_HEADERS = proj.h proj_experimental.h proj_constants.h proj_api.h geodesic.h \ @@ -45,7 +45,7 @@ geodtest_LDADD = libproj.la lib_LTLIBRARIES = libproj.la libproj_la_LDFLAGS = -no-undefined -version-info 14:1:1 -libproj_la_LIBADD = @SQLITE3_LDFLAGS@ +libproj_la_LIBADD = @SQLITE3_LIBS@ libproj_la_SOURCES = \ pj_list.h proj_internal.h proj_math.h \ -- cgit v1.2.3 From d6aaddc5e0cb2bd9080a8fc4b77de947172c81d8 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 4 Feb 2019 14:48:52 +0100 Subject: Remove all traces of nad2bin and nad2nad The source material for the default grids used by PROJ has been moved to the proj-datumgrid repository. For that reason it is no longer necessary to include the nad2bin program in the PROJ repository and source distribution. From now on the nad2bin application will be kept in the proj-datumgrid repo. Previously the null grid was generated by running nad2bin on the null.lla file. Since nad2bin is no longer available null.lla has been replaced by its binary counterpart null. This file will be distributed and installed alongside PROJ. Build scripts and documenation has been adjusted so that nad2bin is not mentioned anywhere. Additionally all references to nad2nad has been removed as well. nad2nad has not been part of the PROJ distribution for quite some time so this has been long overdue. --- src/CMakeLists.txt | 6 - src/Makefile.am | 6 +- src/apps/nad2bin.cpp | 382 -------------------------------------------------- src/bin_nad2bin.cmake | 17 --- src/nad_list.h | 6 - 5 files changed, 2 insertions(+), 415 deletions(-) delete mode 100644 src/apps/nad2bin.cpp delete mode 100644 src/bin_nad2bin.cmake delete mode 100644 src/nad_list.h (limited to 'src') diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index de8ca3e5..122227bf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,7 +7,6 @@ option(BUILD_CCT "Build cct (coordinate conversion and transformation tool)" option(BUILD_CS2CS "Build cs2cs (coordinate systems to coordinate systems translation tool)" ON) option(BUILD_GEOD "Build geod (computation of geodesic lines)" ON) option(BUILD_GIE "Build gie (geospatial integrity investigation environment - a PROJ.4 test tool)" ON) -option(BUILD_NAD2BIN "Build nad2bin (format conversion tool)" ON) option(BUILD_PROJ "Build proj (cartographic projection tool : latlong <-> projected coordinates)" ON) option(BUILD_PROJINFO "Build projinfo (SRS and coordinate operation metadata/query tool)" ON) @@ -50,11 +49,6 @@ if(BUILD_GEOD) set(BIN_TARGETS ${BIN_TARGETS} geod) endif(BUILD_GEOD) -if(BUILD_NAD2BIN) - include(bin_nad2bin.cmake) - set(BIN_TARGETS ${BIN_TARGETS} nad2bin) -endif(BUILD_NAD2BIN) - if(BUILD_PROJ) include(bin_proj.cmake) set(BIN_TARGETS ${BIN_TARGETS} binproj) diff --git a/src/Makefile.am b/src/Makefile.am index fabea4bf..bf1cd32c 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ AM_CFLAGS = @C_WFLAGS@ -bin_PROGRAMS = proj nad2bin geod cs2cs gie cct projinfo +bin_PROGRAMS = proj geod cs2cs gie cct projinfo EXTRA_PROGRAMS = multistresstest test228 TESTS = geodtest @@ -14,7 +14,7 @@ include_HEADERS = proj.h proj_experimental.h proj_constants.h proj_api.h geodesi org_proj4_PJ.h proj_symbol_rename.h EXTRA_DIST = bin_cct.cmake bin_gie.cmake bin_cs2cs.cmake \ - bin_geod.cmake bin_nad2bin.cmake bin_proj.cmake bin_projinfo.cmake \ + bin_geod.cmake bin_proj.cmake bin_projinfo.cmake \ lib_proj.cmake CMakeLists.txt bin_geodtest.cmake tests/geodtest.cpp \ wkt1_grammar.y wkt2_grammar.y apps/emess.h @@ -22,7 +22,6 @@ proj_SOURCES = apps/proj.cpp apps/emess.cpp projinfo_SOURCES = apps/projinfo.cpp cs2cs_SOURCES = apps/cs2cs.cpp apps/emess.cpp cct_SOURCES = apps/cct.cpp apps/proj_strtod.cpp apps/proj_strtod.h apps/optargpm.h -nad2bin_SOURCES = apps/nad2bin.cpp geod_SOURCES = apps/geod.cpp apps/geod_set.cpp apps/geod_interface.cpp apps/geod_interface.h apps/emess.cpp gie_SOURCES = apps/gie.cpp apps/proj_strtod.cpp apps/proj_strtod.h apps/optargpm.h @@ -35,7 +34,6 @@ cs2cs_LDADD = libproj.la geod_LDADD = libproj.la proj_LDADD = libproj.la projinfo_LDADD = libproj.la -nad2bin_LDADD = libproj.la gie_LDADD = libproj.la multistresstest_LDADD = libproj.la @THREAD_LIB@ diff --git a/src/apps/nad2bin.cpp b/src/apps/nad2bin.cpp deleted file mode 100644 index a684b087..00000000 --- a/src/apps/nad2bin.cpp +++ /dev/null @@ -1,382 +0,0 @@ -/* Convert bivariate ASCII NAD27 to NAD83 tables to NTv2 binary structure */ -#include -#include - -#define PJ_LIB__ -#include "proj_internal.h" -#include "proj_internal.h" -#define U_SEC_TO_RAD 4.848136811095359935899141023e-12 - -/************************************************************************/ -/* swap_words() */ -/* */ -/* Convert the byte order of the given word(s) in place. */ -/************************************************************************/ - -static const int byte_order_test = 1; -#define IS_LSB (((const unsigned char *) (&byte_order_test))[0] == 1) - -static void swap_words( void *data_in, int word_size, int word_count ) - -{ - int word; - unsigned char *data = (unsigned char *) data_in; - - for( word = 0; word < word_count; word++ ) - { - int i; - - for( i = 0; i < word_size/2; i++ ) - { - unsigned char t; - - t = data[i]; - data[i] = data[word_size-i-1]; - data[word_size-i-1] = t; - } - - data += word_size; - } -} - -/************************************************************************/ -/* Usage() */ -/************************************************************************/ - -static void Usage() -{ - fprintf(stderr, - "usage: nad2bin [-f ctable/ctable2/ntv2] binary_output < ascii_source\n" ); - exit(1); -} - -/************************************************************************/ -/* main() */ -/************************************************************************/ -int main(int argc, char **argv) { - struct CTABLE ct; - FLP *p, t; - size_t tsize; - int i, j, ichk; - long lam, laml, phi, phil; - FILE *fp; - - const char *output_file = nullptr; - - const char *format = "ctable2"; - const char *GS_TYPE = "SECONDS"; - const char *VERSION = ""; - const char *SYSTEM_F = "NAD27"; - const char *SYSTEM_T = "NAD83"; - const char *SUB_NAME = ""; - const char *CREATED = ""; - const char *UPDATED = ""; - -/* ==================================================================== */ -/* Process arguments. */ -/* ==================================================================== */ - for( i = 1; i < argc; i++ ) - { - if( i < argc-1 && strcmp(argv[i],"-f") == 0 ) - { - format = argv[++i]; - } - else if( output_file == nullptr ) - { - output_file = argv[i]; - } - else - Usage(); - } - - if( output_file == nullptr ) - Usage(); - - fprintf( stdout, "Output Binary File Format: %s\n", format ); - -/* ==================================================================== */ -/* Read the ASCII Table */ -/* ==================================================================== */ - - memset(ct.id,0,MAX_TAB_ID); - if ( nullptr == fgets(ct.id, MAX_TAB_ID, stdin) ) { - perror("fgets"); - exit(1); - } - /* cppcheck-suppress invalidscanf */ - if ( EOF == scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi, - &ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi) ) { - perror("scanf"); - exit(1); - } - if (!(ct.cvs = (FLP *)malloc(tsize = ct.lim.lam * ct.lim.phi * - sizeof(FLP)))) { - perror("mem. alloc"); - exit(1); - } - ct.ll.lam *= DEG_TO_RAD; - ct.ll.phi *= DEG_TO_RAD; - ct.del.lam *= DEG_TO_RAD; - ct.del.phi *= DEG_TO_RAD; - /* load table */ - p = ct.cvs; - for (i = 0; i < ct.lim.phi; ++i) { - /* cppcheck-suppress invalidscanf */ - if ( EOF == scanf("%d:%ld %ld", &ichk, &laml, &phil) ) { - perror("scanf on row"); - exit(1); - } - if (ichk != i) { - fprintf(stderr,"format check on row\n"); - exit(1); - } - t.lam = (float) (laml * U_SEC_TO_RAD); - t.phi = (float) (phil * U_SEC_TO_RAD); - *p++ = t; - for (j = 1; j < ct.lim.lam; ++j) { - /* cppcheck-suppress invalidscanf */ - if ( EOF == scanf("%ld %ld", &lam, &phi) ) { - perror("scanf on column"); - exit(1); - } - t.lam = (float) ((laml += lam) * U_SEC_TO_RAD); - t.phi = (float) ((phil += phi) * U_SEC_TO_RAD); - *p++ = t; - } - } - if (feof(stdin)) { - fprintf(stderr, "premature EOF\n"); - exit(1); - } - -/* ==================================================================== */ -/* Write out the old ctable format - this is machine and byte */ -/* order specific. */ -/* ==================================================================== */ - if( strcmp(format,"ctable") == 0 ) - { - if (!(fp = fopen(output_file, "wb"))) { - perror(output_file); - exit(2); - } - if (fwrite(&ct, sizeof(ct), 1, fp) != 1 || - fwrite(ct.cvs, tsize, 1, fp) != 1) { - fprintf(stderr, "output failure\n"); - exit(2); - } - fclose( fp ); - exit(0); /* normal completion */ - } - -/* ==================================================================== */ -/* Write out the old ctable format - this is machine and byte */ -/* order specific. */ -/* ==================================================================== */ - if( strcmp(format,"ctable2") == 0 ) - { - char header[160]; - - if (!(fp = fopen(output_file, "wb"))) { - perror(output_file); - exit(2); - } - - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( MAX_TAB_ID == 80 ); - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(pj_int32) == 4 ); /* for ct.lim.lam/phi */ - - memset( header, 0, sizeof(header) ); - - memcpy( header + 0, "CTABLE V2.0 ", 16 ); - memcpy( header + 16, ct.id, 80 ); - memcpy( header + 96, &ct.ll.lam, 8 ); - memcpy( header + 104, &ct.ll.phi, 8 ); - memcpy( header + 112, &ct.del.lam, 8 ); - memcpy( header + 120, &ct.del.phi, 8 ); - memcpy( header + 128, &ct.lim.lam, 4 ); - memcpy( header + 132, &ct.lim.phi, 4 ); - - /* force into LSB format */ - if( !IS_LSB ) - { - swap_words( header + 96, 8, 4 ); - swap_words( header + 128, 4, 2 ); - swap_words( ct.cvs, 4, ct.lim.lam * 2 * ct.lim.phi ); - } - - if( fwrite( header, sizeof(header), 1, fp ) != 1 ) { - perror( "fwrite" ); - exit( 2 ); - } - - if (fwrite(ct.cvs, tsize, 1, fp) != 1) { - perror( "fwrite" ); - exit(2); - } - - fclose( fp ); - exit(0); /* normal completion */ - } - -/* ==================================================================== */ -/* Write out the NTv2 format grid shift file. */ -/* ==================================================================== */ - if( strcmp(format,"ntv2") == 0 ) - { - if (!(fp = fopen(output_file, "wb"))) - { - perror(output_file); - exit(2); - } - -/* -------------------------------------------------------------------- */ -/* Write the file header. */ -/* -------------------------------------------------------------------- */ - { - char achHeader[11*16]; - - memset( achHeader, 0, sizeof(achHeader) ); - - memcpy( achHeader + 0*16, "NUM_OREC", 8 ); - achHeader[ 0*16 + 8] = 0xb; - - memcpy( achHeader + 1*16, "NUM_SREC", 8 ); - achHeader[ 1*16 + 8] = 0xb; - - memcpy( achHeader + 2*16, "NUM_FILE", 8 ); - achHeader[ 2*16 + 8] = 0x1; - - memcpy( achHeader + 3*16, "GS_TYPE ", 16 ); - memcpy( achHeader + 3*16+8, GS_TYPE, MIN(16,strlen(GS_TYPE)) ); - - memcpy( achHeader + 4*16, "VERSION ", 16 ); - memcpy( achHeader + 4*16+8, VERSION, MIN(16,strlen(VERSION)) ); - - memcpy( achHeader + 5*16, "SYSTEM_F ", 16 ); - memcpy( achHeader + 5*16+8, SYSTEM_F, MIN(16,strlen(SYSTEM_F)) ); - - memcpy( achHeader + 6*16, "SYSTEM_T ", 16 ); - memcpy( achHeader + 6*16+8, SYSTEM_T, MIN(16,strlen(SYSTEM_T)) ); - - memcpy( achHeader + 7*16, "MAJOR_F ", 8); - memcpy( achHeader + 8*16, "MINOR_F ", 8 ); - memcpy( achHeader + 9*16, "MAJOR_T ", 8 ); - memcpy( achHeader + 10*16, "MINOR_T ", 8 ); - - fwrite( achHeader, 1, sizeof(achHeader), fp ); - } - -/* -------------------------------------------------------------------- */ -/* Write the grid header. */ -/* -------------------------------------------------------------------- */ - { - unsigned char achHeader[11*16]; - double dfValue; - pj_int32 nGSCount = ct.lim.lam * ct.lim.phi; - PJ_LP ur; - - ur.lam = ct.ll.lam + (ct.lim.lam-1) * ct.del.lam; - ur.phi = ct.ll.phi + (ct.lim.phi-1) * ct.del.phi; - - /* cppcheck-suppress sizeofCalculation */ - STATIC_ASSERT( sizeof(nGSCount) == 4 ); - - memset( achHeader, 0, sizeof(achHeader) ); - - memcpy( achHeader + 0*16, "SUB_NAME ", 16 ); - memcpy( achHeader + 0*16+8, SUB_NAME, MIN(16,strlen(SUB_NAME)) ); - - memcpy( achHeader + 1*16, "PARENT ", 16 ); - memcpy( achHeader + 1*16+8, "NONE", MIN(16,strlen("NONE")) ); - - memcpy( achHeader + 2*16, "CREATED ", 16 ); - memcpy( achHeader + 2*16+8, CREATED, MIN(16,strlen(CREATED)) ); - - memcpy( achHeader + 3*16, "UPDATED ", 16 ); - memcpy( achHeader + 3*16+8, UPDATED, MIN(16,strlen(UPDATED)) ); - - memcpy( achHeader + 4*16, "S_LAT ", 8 ); - dfValue = ct.ll.phi * 3600.0 / DEG_TO_RAD; - memcpy( achHeader + 4*16 + 8, &dfValue, 8 ); - - memcpy( achHeader + 5*16, "N_LAT ", 8 ); - dfValue = ur.phi * 3600.0 / DEG_TO_RAD; - memcpy( achHeader + 5*16 + 8, &dfValue, 8 ); - - memcpy( achHeader + 6*16, "E_LONG ", 8 ); - dfValue = -1 * ur.lam * 3600.0 / DEG_TO_RAD; - memcpy( achHeader + 6*16 + 8, &dfValue, 8 ); - - memcpy( achHeader + 7*16, "W_LONG ", 8 ); - dfValue = -1 * ct.ll.lam * 3600.0 / DEG_TO_RAD; - memcpy( achHeader + 7*16 + 8, &dfValue, 8 ); - - memcpy( achHeader + 8*16, "LAT_INC ", 8 ); - dfValue = ct.del.phi * 3600.0 / DEG_TO_RAD; - memcpy( achHeader + 8*16 + 8, &dfValue, 8 ); - - memcpy( achHeader + 9*16, "LONG_INC", 8 ); - dfValue = ct.del.lam * 3600.0 / DEG_TO_RAD; - memcpy( achHeader + 9*16 + 8, &dfValue, 8 ); - - memcpy( achHeader + 10*16, "GS_COUNT", 8 ); - memcpy( achHeader + 10*16+8, &nGSCount, 4 ); - - if( !IS_LSB ) - { - swap_words( achHeader + 4*16 + 8, 8, 1 ); - swap_words( achHeader + 5*16 + 8, 8, 1 ); - swap_words( achHeader + 6*16 + 8, 8, 1 ); - swap_words( achHeader + 7*16 + 8, 8, 1 ); - swap_words( achHeader + 8*16 + 8, 8, 1 ); - swap_words( achHeader + 9*16 + 8, 8, 1 ); - swap_words( achHeader + 10*16 + 8, 4, 1 ); - } - - fwrite( achHeader, 1, sizeof(achHeader), fp ); - } - -/* -------------------------------------------------------------------- */ -/* Write the actual grid cells. */ -/* -------------------------------------------------------------------- */ - { - float *row_buf; - int row; - - row_buf = (float *) pj_malloc(ct.lim.lam * sizeof(float) * 4); - memset( row_buf, 0, sizeof(float)*4 ); - - for( row = 0; row < ct.lim.phi; row++ ) - { - for( i = 0; i < ct.lim.lam; i++ ) - { - FLP *cvs = ct.cvs + (row) * ct.lim.lam - + (ct.lim.lam - i - 1); - - /* convert radians to seconds */ - row_buf[i*4+0] = (float) (cvs->phi * (3600.0 / (M_PI/180.0))); - row_buf[i*4+1] = (float) (cvs->lam * (3600.0 / (M_PI/180.0))); - - /* We leave the accuracy values as zero */ - } - - if( !IS_LSB ) - swap_words( row_buf, 4, ct.lim.lam * 4 ); - - if( fwrite( row_buf, sizeof(float), ct.lim.lam*4, fp ) - != (size_t)( 4 * ct.lim.lam ) ) - { - perror( "write()" ); - exit( 2 ); - } - } - } - - fclose( fp ); - exit(0); /* normal completion */ - } - - fprintf( stderr, "Unsupported format, nothing written.\n" ); - exit( 3 ); -} diff --git a/src/bin_nad2bin.cmake b/src/bin_nad2bin.cmake deleted file mode 100644 index 271aac93..00000000 --- a/src/bin_nad2bin.cmake +++ /dev/null @@ -1,17 +0,0 @@ -if(WIN32 AND BUILD_LIBPROJ_SHARED) - message(warning " nad2nad can't be build with a DLL proj library you need a static lib") -endif(WIN32 AND BUILD_LIBPROJ_SHARED) - - -set(NAD2BIN_SRC apps/nad2bin.cpp) -source_group("Source Files\\Bin" FILES ${NAD2BIN_SRC}) - -#Executable -add_executable(nad2bin ${NAD2BIN_SRC}) -target_link_libraries(nad2bin ${PROJ_LIBRARIES}) -install(TARGETS nad2bin - RUNTIME DESTINATION ${BINDIR}) - -if(MSVC AND BUILD_LIBPROJ_SHARED) - target_compile_definitions(nad2bin PRIVATE PROJ_MSVC_DLL_IMPORT=1) -endif() diff --git a/src/nad_list.h b/src/nad_list.h deleted file mode 100644 index f82a2ab7..00000000 --- a/src/nad_list.h +++ /dev/null @@ -1,6 +0,0 @@ -/* projection list for program nad2nad */ -PROJ_HEAD(lcc, "Lambert Conformal Conic") -PROJ_HEAD(omerc, "Oblique Mercator") -PROJ_HEAD(poly, "Polyconic (American)") -PROJ_HEAD(tmerc, "Transverse Mercator") -PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") -- cgit v1.2.3 From 149bd81691e309d4ab22bed944ea69fbeaec450f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Tue, 5 Feb 2019 23:48:32 +0100 Subject: PROJStringParser::createFromPROJString(): avoid potential infinite loop. Fixes https://bugs.chromium.org/p/oss-fuzz/issues/detail?id=12913. Credit to OSS Fuzz --- src/iso19111/io.cpp | 16 ++++++++++++++++ src/proj_internal.h | 2 ++ 2 files changed, 18 insertions(+) (limited to 'src') diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp index 1f4a7c8b..431c75af 100644 --- a/src/iso19111/io.cpp +++ b/src/iso19111/io.cpp @@ -5313,6 +5313,9 @@ PROJStringSyntaxParser(const std::string &projString, std::vector &steps, } token += c_str[i]; } + if (in_string) { + throw ParsingException("Unbalanced double quote"); + } if (token.empty()) { break; } @@ -7452,6 +7455,13 @@ static const metadata::ExtentPtr &getExtent(const crs::CRS *crs) { */ BaseObjectNNPtr PROJStringParser::createFromPROJString(const std::string &projString) { + + // In some abnormal situations involving init=epsg:XXXX syntax, we could + // have infinite loop + if (d->ctx_ && d->ctx_->curStringInCreateFromPROJString == projString) { + throw ParsingException("invalid PROJ string"); + } + d->steps_.clear(); d->title_.clear(); d->globalParamValues_.clear(); @@ -7760,11 +7770,17 @@ PROJStringParser::createFromPROJString(const std::string &projString) { proj_log_func(pj_context, &logger, Logger::log); proj_context_use_proj4_init_rules(pj_context, d->usePROJ4InitRules_); } + if (d->ctx_) { + d->ctx_->curStringInCreateFromPROJString = projString; + } auto pj = pj_create_internal( pj_context, (projString.find("type=crs") != std::string::npos ? projString + " +disable_grid_presence_check" : projString) .c_str()); + if (d->ctx_) { + d->ctx_->curStringInCreateFromPROJString.clear(); + } valid = pj != nullptr; // Remove parameters not understood by PROJ. diff --git a/src/proj_internal.h b/src/proj_internal.h index f5196939..453bd654 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -702,6 +702,8 @@ struct projCtx_t { const char* (*file_finder) (PJ_CONTEXT *, const char*, void* user_data) = nullptr; void* file_finder_user_data = nullptr; + std::string curStringInCreateFromPROJString{}; + projCtx_t() = default; projCtx_t(const projCtx_t&); ~projCtx_t(); -- cgit v1.2.3 From 457b173cbe8fdb790b011d1828a0fd9f8f6221a4 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 6 Feb 2019 14:26:06 +0100 Subject: ISO19111: Handle database area objects with no bounding box --- src/iso19111/factory.cpp | 6 ++++++ 1 file changed, 6 insertions(+) (limited to 'src') diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp index c862cc5b..a3047f8b 100644 --- a/src/iso19111/factory.cpp +++ b/src/iso19111/factory.cpp @@ -1485,6 +1485,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(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]); -- cgit v1.2.3 From 02efe72181814097284196de9b9b984db0fb3d26 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 6 Feb 2019 18:32:58 +0100 Subject: Add proj_get_crs_info_list_from_database() This method is intended to be used typically by GUI that lists all possible CRS. What is does could be done by previously existing functions, but it is much faster. It typically runs in less than 0.1s (hot run) versus ~0.5s with the method that consists in enumerating all codes and instanciating a PJ object for each of them. --- src/iso19111/c_api.cpp | 212 +++++++++++++++++++++++++++++++++++++++++++++++ src/iso19111/factory.cpp | 161 ++++++++++++++++++++++++++++++----- src/proj.h | 81 ++++++++++++++++++ 3 files changed, 431 insertions(+), 23 deletions(-) (limited to 'src') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index de11f181..5982da47 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -1895,6 +1896,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 +1935,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_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_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_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_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/factory.cpp b/src/iso19111/factory.cpp index a3047f8b..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()) { @@ -2067,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()) { @@ -2124,15 +2133,14 @@ AuthorityFactory::createGeodeticCRS(const std::string &code, auto ellipsoidalCS = util::nn_dynamic_pointer_cast(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); - 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); @@ -2477,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); @@ -3782,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 "; @@ -3858,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::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 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. @@ -4056,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/proj.h b/src/proj.h index 552f4b26..4f5fb1d1 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_list_destroy(PROJ_CRS_INFO** list); + /* ------------------------------------------------------------------------- */ -- cgit v1.2.3 From 96fe92361156c0b242301dd8cb56115cec243af4 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 7 Feb 2019 15:05:31 +0100 Subject: Rename proj_crs_list_destroy() to proj_crs_info_list_destroy() --- src/iso19111/c_api.cpp | 8 ++++---- src/proj.h | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 5982da47..c205d039 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -1969,7 +1969,7 @@ void proj_get_crs_list_parameters_destroy(PROJ_CRS_LIST_PARAMETERS *params) { * 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_list_destroy() + * 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 @@ -1985,7 +1985,7 @@ void proj_get_crs_list_parameters_destroy(PROJ_CRS_LIST_PARAMETERS *params) { * @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_list_destroy(), or NULL in case of error. + * 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, @@ -2115,7 +2115,7 @@ proj_get_crs_info_list_from_database(PJ_CONTEXT *ctx, const char *auth_name, proj_log_error(ctx, __FUNCTION__, e.what()); if (ret) { ret[i + 1] = nullptr; - proj_crs_list_destroy(ret); + proj_crs_info_list_destroy(ret); } if (out_result_count) *out_result_count = 0; @@ -2128,7 +2128,7 @@ proj_get_crs_info_list_from_database(PJ_CONTEXT *ctx, const char *auth_name, /** \brief Destroy the result returned by * proj_get_crs_info_list_from_database(). */ -void proj_crs_list_destroy(PROJ_CRS_INFO **list) { +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); diff --git a/src/proj.h b/src/proj.h index 4f5fb1d1..f25c3228 100644 --- a/src/proj.h +++ b/src/proj.h @@ -853,7 +853,7 @@ PROJ_CRS_INFO PROJ_DLL **proj_get_crs_info_list_from_database( const PROJ_CRS_LIST_PARAMETERS* params, int *out_result_count); -void PROJ_DLL proj_crs_list_destroy(PROJ_CRS_INFO** list); +void PROJ_DLL proj_crs_info_list_destroy(PROJ_CRS_INFO** list); /* ------------------------------------------------------------------------- */ -- cgit v1.2.3 From c45d7693ee799fe7317550ffb3dd9a40c4b0a17c Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Sun, 10 Feb 2019 12:06:43 +0100 Subject: Completely remove Chebychev remains from codebase c6ab83f5742bc5ac6f9cb9a8b2a4f1ea241b6f63 already removed their availability in user facing application, but the library code remained, and appeared to be unused by the rest of the library, and not available to library users, the API being only in proj_internal.h. So remove all remains. --- src/Makefile.am | 6 +- src/apps/gie.cpp | 6 +- src/bch2bps.cpp | 154 ----------------------------------------- src/bchgen.cpp | 59 ---------------- src/biveval.cpp | 88 ------------------------ src/lib_proj.cmake | 6 +- src/mk_cheby.cpp | 193 ---------------------------------------------------- src/proj_internal.h | 30 +------- src/strerrno.cpp | 2 +- src/vector1.cpp | 30 -------- 10 files changed, 14 insertions(+), 560 deletions(-) delete mode 100644 src/bch2bps.cpp delete mode 100644 src/bchgen.cpp delete mode 100644 src/biveval.cpp delete mode 100644 src/mk_cheby.cpp delete mode 100644 src/vector1.cpp (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index 43bd91df..2af53e30 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -184,15 +184,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/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/lib_proj.cmake b/src/lib_proj.cmake index b595d87e..cd380b24 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -201,15 +201,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_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/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 -#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; -} -- cgit v1.2.3 From 048cd7f3fc2bee07b28ca1b9418a05ff9f35d35f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 11 Feb 2019 21:11:28 +0100 Subject: Use pj_new() --- src/iso19111/c_api.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index c205d039..6fafa2c8 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -145,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 -- cgit v1.2.3 From dd1075784fecfbfa929787dbed271e876ca1693f Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Mon, 11 Feb 2019 22:11:01 +0100 Subject: cct: fix memory leak --- src/apps/cct.cpp | 2 ++ 1 file changed, 2 insertions(+) (limited to 'src') 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); -- cgit v1.2.3 From 5141b3908e59a26c9fe66de94bb7388bff741b58 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 11 Feb 2019 23:58:16 +0100 Subject: Make tmerc an alias for etmerc. (#1234) * Make tmerc an alias for etmerc This switches the algorithm used in tmerc to the Poder/Engsager tmerc algorithm. The original tmerc algorithm of Evenden/Snyder origin can still be accessed by adding the +approx flag when initializing a tmerc projection. The +approx flag can also be used when initializing UTM projections, in which case the Evenden/Snyder algorithm is used as well. If a tmerc projection is instantiated on a spherical earth the Evenden/Snyder algorithm is used as well since the Poder/Engsager algorithm is only defined on the ellipsoid. +proj=etmerc can still be instantiated for backwards compatibility reasons. Co-authored-by: Kristian Evers Co-authored-by: Even Rouault --- src/Makefile.am | 1 - src/iso19111/c_api.cpp | 11 +- src/iso19111/coordinateoperation.cpp | 27 ++- src/iso19111/io.cpp | 22 +- src/lib_proj.cmake | 1 - src/projections/etmerc.cpp | 362 ------------------------------ src/projections/tmerc.cpp | 423 +++++++++++++++++++++++++++++++++-- 7 files changed, 429 insertions(+), 418 deletions(-) delete mode 100644 src/projections/etmerc.cpp (limited to 'src') diff --git a/src/Makefile.am b/src/Makefile.am index 2af53e30..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 \ diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp index 6fafa2c8..240d11b3 100644 --- a/src/iso19111/c_api.cpp +++ b/src/iso19111/c_api.cpp @@ -1207,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, @@ -1243,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()); 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/io.cpp b/src/iso19111/io.cpp index 431c75af..60c28201 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 cd380b24..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 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 - -#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(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(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(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(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(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 @@ -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(P->opaque); + struct pj_opaque_approx *Q = static_cast(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(P->opaque)->ml0 * log ((1. + b) / (1. - b)); + xy.x = static_cast(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(P->opaque)->esp * (xy.y - P->phi0); + xy.y = static_cast(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(P->opaque); + struct pj_opaque_approx *Q = static_cast(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(P->opaque)->esp); + h = exp(xy.x / static_cast(P->opaque)->esp); g = .5 * (h - 1. / h); - h = cos (P->phi0 + xy.y / static_cast(P->opaque)->esp); + h = cos (P->phi0 + xy.y / static_cast(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(P->opaque)->en); + pj_dealloc (static_cast(P->opaque)->en); return pj_default_destructor(P, errlev); } -static PJ *setup(PJ *P) { /* general initialization */ - struct pj_opaque *Q = static_cast(P->opaque); +static PJ *setup_approx(PJ *P) { + struct pj_opaque_approx *Q = static_cast(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(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(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(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(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(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(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(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(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(pj_calloc (1, sizeof (struct pj_opaque_exact))); + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = Q; + + return setup_exact(P); + } } -- cgit v1.2.3