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