aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/4D_api.cpp309
-rw-r--r--src/apps/cs2cs.cpp44
-rw-r--r--src/apps/geod.cpp17
-rw-r--r--src/apps/geod_set.cpp24
-rw-r--r--src/apps/proj.cpp17
-rw-r--r--src/conversions/unitconvert.cpp6
-rw-r--r--src/grids.cpp13
-rw-r--r--src/init.cpp2
-rw-r--r--src/iso19111/c_api.cpp243
-rw-r--r--src/iso19111/common.cpp5
-rw-r--r--src/iso19111/coordinateoperation.cpp2
-rw-r--r--src/iso19111/crs.cpp43
-rw-r--r--src/iso19111/factory.cpp58
-rw-r--r--src/iso19111/io.cpp39
-rw-r--r--src/lib_proj.cmake2
-rw-r--r--src/proj.h69
-rw-r--r--src/proj_internal.h109
-rw-r--r--src/projections/stere.cpp5
-rw-r--r--src/projections/tmerc.cpp184
-rw-r--r--src/units.cpp12
20 files changed, 852 insertions, 351 deletions
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 52bd7f4f..6e494793 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -179,7 +179,58 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord) {
return proj_xyz_dist (org, t);
}
+/**************************************************************************************/
+int pj_get_suggested_operation(PJ_CONTEXT*,
+ const std::vector<CoordOperation>& opList,
+ const int iExcluded[2],
+ PJ_DIRECTION direction,
+ PJ_COORD coord)
+/**************************************************************************************/
+{
+ // Select the operations that match the area of use
+ // and has the best accuracy.
+ int iBest = -1;
+ double bestAccuracy = std::numeric_limits<double>::max();
+ const int nOperations = static_cast<int>(opList.size());
+ for( int i = 0; i < nOperations; i++ ) {
+ if( i == iExcluded[0] || i == iExcluded[1] ) {
+ continue;
+ }
+ const auto &alt = opList[i];
+ bool spatialCriterionOK = false;
+ if( direction == PJ_FWD ) {
+ if( coord.xyzt.x >= alt.minxSrc &&
+ coord.xyzt.y >= alt.minySrc &&
+ coord.xyzt.x <= alt.maxxSrc &&
+ coord.xyzt.y <= alt.maxySrc) {
+ spatialCriterionOK = true;
+ }
+ } else {
+ if( coord.xyzt.x >= alt.minxDst &&
+ coord.xyzt.y >= alt.minyDst &&
+ coord.xyzt.x <= alt.maxxDst &&
+ coord.xyzt.y <= alt.maxyDst ) {
+ spatialCriterionOK = true;
+ }
+ }
+ if( spatialCriterionOK ) {
+ // The offshore test is for the "Test bug 245 (use +datum=carthage)"
+ // of testvarious. The long=10 lat=34 point belongs both to the
+ // onshore and offshore Tunisia area of uses, but is slightly
+ // onshore. So in a general way, prefer a onshore area to a
+ // offshore one.
+ if( iBest < 0 ||
+ (alt.accuracy >= 0 && alt.accuracy < bestAccuracy &&
+ !alt.isOffshore) ) {
+ iBest = i;
+ bestAccuracy = alt.accuracy;
+ }
+ }
+ }
+
+ return iBest;
+}
/**************************************************************************************/
PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coord) {
@@ -211,45 +262,11 @@ similarly, but prefers the 2D resp. 3D interfaces if available.
{
// Do a first pass and select the operations that match the area of use
// and has the best accuracy.
- int iBest = -1;
- double bestAccuracy = std::numeric_limits<double>::max();
- for( int i = 0; i < nOperations; i++ ) {
- if( i == iExcluded[0] || i == iExcluded[1] ) {
- continue;
- }
- const auto &alt = P->alternativeCoordinateOperations[i];
- bool spatialCriterionOK = false;
- if( direction == PJ_FWD ) {
- if( coord.xyzt.x >= alt.minxSrc &&
- coord.xyzt.y >= alt.minySrc &&
- coord.xyzt.x <= alt.maxxSrc &&
- coord.xyzt.y <= alt.maxySrc) {
- spatialCriterionOK = true;
- }
- } else {
- if( coord.xyzt.x >= alt.minxDst &&
- coord.xyzt.y >= alt.minyDst &&
- coord.xyzt.x <= alt.maxxDst &&
- coord.xyzt.y <= alt.maxyDst ) {
- spatialCriterionOK = true;
- }
- }
-
- if( spatialCriterionOK ) {
- // The offshore test is for the "Test bug 245 (use +datum=carthage)"
- // of testvarious. The long=10 lat=34 point belongs both to the
- // onshore and offshore Tunisia area of uses, but is slightly
- // onshore. So in a general way, prefer a onshore area to a
- // offshore one.
- if( iBest < 0 ||
- (alt.accuracy >= 0 && alt.accuracy < bestAccuracy &&
- !alt.isOffshore) ) {
- iBest = i;
- bestAccuracy = alt.accuracy;
- }
- }
- }
-
+ int iBest = pj_get_suggested_operation(P->ctx,
+ P->alternativeCoordinateOperations,
+ iExcluded,
+ direction,
+ coord);
if( iBest < 0 ) {
break;
}
@@ -968,13 +985,15 @@ static void reproject_bbox(PJ* pjGeogToCrs,
/*****************************************************************************/
-static PJ* add_coord_op_to_list(PJ* op,
+static PJ* add_coord_op_to_list(
+ int idxInOriginalList,
+ PJ* op,
double west_lon, double south_lat,
double east_lon, double north_lat,
PJ* pjGeogToSrc,
PJ* pjGeogToDst,
bool isOffshore,
- std::vector<PJconsts::CoordOperation>& altCoordOps) {
+ std::vector<CoordOperation>& altCoordOps) {
/*****************************************************************************/
double minxSrc;
@@ -997,9 +1016,10 @@ static PJ* add_coord_op_to_list(PJ* op,
std::string name(c_name ? c_name : "");
const double accuracy = proj_coordoperation_get_accuracy(op->ctx, op);
- altCoordOps.emplace_back(minxSrc, minySrc, maxxSrc, maxySrc,
- minxDst, minyDst, maxxDst, maxyDst,
- op, name, accuracy, isOffshore);
+ altCoordOps.emplace_back(idxInOriginalList,
+ minxSrc, minySrc, maxxSrc, maxySrc,
+ minxDst, minyDst, maxxDst, maxyDst,
+ op, name, accuracy, isOffshore);
op = nullptr;
}
return op;
@@ -1140,6 +1160,91 @@ PJ *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *source_crs, const char
return ret;
}
+
+/*****************************************************************************/
+std::vector<CoordOperation> pj_create_prepared_operations(PJ_CONTEXT *ctx,
+ const PJ *source_crs,
+ const PJ *target_crs,
+ PJ_OBJ_LIST* op_list)
+/*****************************************************************************/
+{
+ auto pjGeogToSrc = create_operation_to_geog_crs(ctx, source_crs);
+ if( !pjGeogToSrc )
+ {
+ proj_context_log_debug(ctx,
+ "Cannot create transformation from geographic CRS of source CRS to source CRS");
+ return {};
+ }
+
+ auto pjGeogToDst = create_operation_to_geog_crs(ctx, target_crs);
+ if( !pjGeogToDst )
+ {
+ proj_context_log_debug(ctx,
+ "Cannot create transformation from geographic CRS of target CRS to target CRS");
+ proj_destroy(pjGeogToSrc);
+ return {};
+ }
+
+ try
+ {
+ std::vector<CoordOperation> preparedOpList;
+
+ // Iterate over source->target candidate transformations and reproject
+ // their long-lat bounding box into the source CRS.
+ const auto op_count = proj_list_get_count(op_list);
+ for( int i = 0; i < op_count; i++ )
+ {
+ auto op = proj_list_get(ctx, op_list, i);
+ assert(op);
+ double west_lon = 0.0;
+ double south_lat = 0.0;
+ double east_lon = 0.0;
+ double north_lat = 0.0;
+
+ const char* areaName = nullptr;
+ if( proj_get_area_of_use(ctx, op,
+ &west_lon, &south_lat, &east_lon, &north_lat, &areaName) )
+ {
+ const bool isOffshore =
+ areaName && strstr(areaName, "offshore");
+ if( west_lon <= east_lon )
+ {
+ op = add_coord_op_to_list(i, op,
+ west_lon, south_lat, east_lon, north_lat,
+ pjGeogToSrc, pjGeogToDst, isOffshore,
+ preparedOpList);
+ }
+ else
+ {
+ auto op_clone = proj_clone(ctx, op);
+
+ op = add_coord_op_to_list(i, op,
+ west_lon, south_lat, 180, north_lat,
+ pjGeogToSrc, pjGeogToDst, isOffshore,
+ preparedOpList);
+ op_clone = add_coord_op_to_list(i, op_clone,
+ -180, south_lat, east_lon, north_lat,
+ pjGeogToSrc, pjGeogToDst, isOffshore,
+ preparedOpList);
+ proj_destroy(op_clone);
+ }
+ }
+
+ proj_destroy(op);
+ }
+
+ proj_destroy(pjGeogToSrc);
+ proj_destroy(pjGeogToDst);
+ return preparedOpList;
+ }
+ catch( const std::exception& )
+ {
+ proj_destroy(pjGeogToSrc);
+ proj_destroy(pjGeogToDst);
+ return {};
+ }
+}
+
/*****************************************************************************/
PJ *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, const PJ *target_crs, PJ_AREA *area, const char* const *) {
/******************************************************************************
@@ -1177,16 +1282,15 @@ PJ *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
auto op_list = proj_create_operations(ctx, source_crs, target_crs, operation_ctx);
+ proj_operation_factory_context_destroy(operation_ctx);
if( !op_list ) {
- proj_operation_factory_context_destroy(operation_ctx);
return nullptr;
}
auto op_count = proj_list_get_count(op_list);
if( op_count == 0 ) {
proj_list_destroy(op_list);
- proj_operation_factory_context_destroy(operation_ctx);
proj_context_log_debug(ctx, "No operation found matching criteria");
return nullptr;
@@ -1199,112 +1303,39 @@ PJ *proj_create_crs_to_crs_from_pj (PJ_CONTEXT *ctx, const PJ *source_crs, cons
proj_get_type(source_crs) == PJ_TYPE_GEOCENTRIC_CRS ||
proj_get_type(target_crs) == PJ_TYPE_GEOCENTRIC_CRS ) {
proj_list_destroy(op_list);
- proj_operation_factory_context_destroy(operation_ctx);
return P;
}
- auto pjGeogToSrc = create_operation_to_geog_crs(ctx, source_crs);
- if( !pjGeogToSrc )
+ auto preparedOpList = pj_create_prepared_operations(ctx, source_crs, target_crs,
+ op_list);
+ proj_list_destroy(op_list);
+
+ if( preparedOpList.empty() )
{
- proj_list_destroy(op_list);
- proj_operation_factory_context_destroy(operation_ctx);
- proj_context_log_debug(ctx,
- "Cannot create transformation from geographic CRS of source CRS to source CRS");
proj_destroy(P);
return nullptr;
}
- auto pjGeogToDst = create_operation_to_geog_crs(ctx, target_crs);
- if( !pjGeogToDst )
+ // If there's finally juste a single result, return it directly
+ if( preparedOpList.size() == 1 )
{
- proj_list_destroy(op_list);
- proj_operation_factory_context_destroy(operation_ctx);
- proj_context_log_debug(ctx,
- "Cannot create transformation from geographic CRS of target CRS to target CRS");
+ auto retP = preparedOpList[0].pj;
+ preparedOpList[0].pj = nullptr;
proj_destroy(P);
- proj_destroy(pjGeogToSrc);
- return nullptr;
+ return retP;
}
- try
- {
- // Iterate over source->target candidate transformations and reproject
- // their long-lat bounding box into the source CRS.
- for( int i = 0; i < op_count; i++ )
- {
- auto op = proj_list_get(ctx, op_list, i);
- assert(op);
- double west_lon = 0.0;
- double south_lat = 0.0;
- double east_lon = 0.0;
- double north_lat = 0.0;
+ P->alternativeCoordinateOperations = std::move(preparedOpList);
+ // The returned P is rather dummy
+ P->iso_obj = nullptr;
+ P->fwd = nullptr;
+ P->inv = nullptr;
+ P->fwd3d = nullptr;
+ P->inv3d = nullptr;
+ P->fwd4d = nullptr;
+ P->inv4d = nullptr;
- const char* areaName = nullptr;
- if( proj_get_area_of_use(ctx, op,
- &west_lon, &south_lat, &east_lon, &north_lat, &areaName) )
- {
- const bool isOffshore =
- areaName && strstr(areaName, "offshore");
- if( west_lon <= east_lon )
- {
- op = add_coord_op_to_list(op,
- west_lon, south_lat, east_lon, north_lat,
- pjGeogToSrc, pjGeogToDst, isOffshore,
- P->alternativeCoordinateOperations);
- }
- else
- {
- auto op_clone = proj_clone(ctx, op);
-
- op = add_coord_op_to_list(op,
- west_lon, south_lat, 180, north_lat,
- pjGeogToSrc, pjGeogToDst, isOffshore,
- P->alternativeCoordinateOperations);
- op_clone = add_coord_op_to_list(op_clone,
- -180, south_lat, east_lon, north_lat,
- pjGeogToSrc, pjGeogToDst, isOffshore,
- P->alternativeCoordinateOperations);
- proj_destroy(op_clone);
- }
- }
-
- proj_destroy(op);
- }
-
- proj_list_destroy(op_list);
-
- proj_operation_factory_context_destroy(operation_ctx);
- proj_destroy(pjGeogToSrc);
- proj_destroy(pjGeogToDst);
-
- // If there's finally juste a single result, return it directly
- if( P->alternativeCoordinateOperations.size() == 1 ) {
- auto retP = P->alternativeCoordinateOperations[0].pj;
- P->alternativeCoordinateOperations[0].pj = nullptr;
- proj_destroy(P);
- P = retP;
- } else {
- // The returned P is rather dummy
- P->iso_obj = nullptr;
- P->fwd = nullptr;
- P->inv = nullptr;
- P->fwd3d = nullptr;
- P->inv3d = nullptr;
- P->fwd4d = nullptr;
- P->inv4d = nullptr;
- }
-
- return P;
- }
- catch( const std::exception& )
- {
- proj_list_destroy(op_list);
- proj_operation_factory_context_destroy(operation_ctx);
- proj_destroy(pjGeogToSrc);
- proj_destroy(pjGeogToDst);
- proj_destroy(P);
- return nullptr;
- }
+ return P;
}
PJ *proj_destroy (PJ *P) {
diff --git a/src/apps/cs2cs.cpp b/src/apps/cs2cs.cpp
index 6c85d4aa..c99e16e0 100644
--- a/src/apps/cs2cs.cpp
+++ b/src/apps/cs2cs.cpp
@@ -349,23 +349,14 @@ int main(int argc, char **argv) {
exit(0);
}
- // First pass to check if we have "cs2cs [-bla]* <SRC> <DEST>" syntax
- int countNonOptionArg = 0;
+ // First pass to check if we have "cs2cs [-bla]* <SRC> <DEST> [<filename>]" syntax
+ bool isProj4StyleSyntax = false;
for (int i = 1; i < argc; i++) {
- if (argv[i][0] == '-') {
- if (argv[i][1] == 'f' || argv[i][1] == 'e' || argv[i][1] == 'd' ||
- argv[i][1] == 'D' ) {
- i++;
- }
- } else {
- if (strcmp(argv[i], "+to") == 0) {
- countNonOptionArg = -1;
- break;
- }
- countNonOptionArg++;
+ if (argv[i][0] == '+') {
+ isProj4StyleSyntax = true;
+ break;
}
}
- const bool isSrcDestSyntax = (countNonOptionArg == 2);
/* process run line arguments */
while (--argc > 0) { /* collect run line arguments */
@@ -425,11 +416,18 @@ int main(int argc, char **argv) {
(void)printf("%9s %-16s %-16s %s\n", le->id,
le->major, le->ell, le->name);
} else if (arg[1] == 'u') { /* list units */
- const struct PJ_UNITS *lu;
-
- for (lu = proj_list_units(); lu->id; ++lu)
- (void)printf("%12s %-20s %s\n", lu->id,
- lu->to_meter, lu->name);
+ auto units = proj_get_units_from_database(nullptr, nullptr, "linear", false, nullptr);
+ for( int i = 0; units && units[i]; i++ )
+ {
+ if( units[i]->proj_short_name )
+ {
+ (void)printf("%12s %-20.15g %s\n",
+ units[i]->proj_short_name,
+ units[i]->conv_factor,
+ units[i]->name);
+ }
+ }
+ proj_unit_list_destroy(units);
} else if (arg[1] == 'm') { /* list prime meridians */
const struct PJ_PRIME_MERIDIANS *lpm;
@@ -485,11 +483,15 @@ int main(int argc, char **argv) {
}
break;
}
- } else if (isSrcDestSyntax) {
+ } else if (!isProj4StyleSyntax) {
if (fromStr.empty())
fromStr = *argv;
- else
+ else if( toStr.empty() )
toStr = *argv;
+ else {
+ /* assumed to be input file name(s) */
+ eargv[eargc++] = *argv;
+ }
} else if (strcmp(*argv, "+to") == 0) {
have_to_flag = 1;
diff --git a/src/apps/geod.cpp b/src/apps/geod.cpp
index b46188d3..919430ca 100644
--- a/src/apps/geod.cpp
+++ b/src/apps/geod.cpp
@@ -185,11 +185,18 @@ noargument: emess(1,"missing argument for -%c",*arg);
(void)printf("%9s %-16s %-16s %s\n",
le->id, le->major, le->ell, le->name);
} else if (arg[1] == 'u') { /* list of units */
- const struct PJ_UNITS *lu;
-
- for (lu = proj_list_units();lu->id ; ++lu)
- (void)printf("%12s %-20s %s\n",
- lu->id, lu->to_meter, lu->name);
+ auto units = proj_get_units_from_database(nullptr, nullptr, "linear", false, nullptr);
+ for( int i = 0; units && units[i]; i++ )
+ {
+ if( units[i]->proj_short_name )
+ {
+ (void)printf("%12s %-20.15g %s\n",
+ units[i]->proj_short_name,
+ units[i]->conv_factor,
+ units[i]->name);
+ }
+ }
+ proj_unit_list_destroy(units);
} else
emess(1,"invalid list option: l%c",arg[1]);
exit( 0 );
diff --git a/src/apps/geod_set.cpp b/src/apps/geod_set.cpp
index ed7edeb9..603f0d95 100644
--- a/src/apps/geod_set.cpp
+++ b/src/apps/geod_set.cpp
@@ -14,7 +14,6 @@ geod_set(int argc, char **argv) {
paralist *start = nullptr, *curr;
double es;
char *name;
- int i;
/* put arguments into internal linked list */
if (argc <= 0)
@@ -22,7 +21,7 @@ geod_set(int argc, char **argv) {
start = curr = pj_mkparam(argv[0]);
if (!curr)
emess(1, "memory allocation failed");
- for (i = 1; curr != nullptr && i < argc; ++i) {
+ for (int i = 1; curr != nullptr && i < argc; ++i) {
curr->next = pj_mkparam(argv[i]);
if (!curr->next)
emess(1, "memory allocation failed");
@@ -32,13 +31,20 @@ geod_set(int argc, char **argv) {
if (pj_ell_set(pj_get_default_ctx(),start, &geod_a, &es)) emess(1,"ellipse setup failure");
/* set units */
if ((name = pj_param(nullptr,start, "sunits").s) != nullptr) {
- const char *s;
- const struct PJ_UNITS *unit_list = proj_list_units();
- for (i = 0; (s = unit_list[i].id) && strcmp(name, s) ; ++i) ;
- if (!s)
- emess(1,"%s unknown unit conversion id", name);
- to_meter = unit_list[i].factor;
- fr_meter = 1 / to_meter;
+ bool unit_found = false;
+ auto units = proj_get_units_from_database(nullptr, nullptr, "linear", false, nullptr);
+ for( int i = 0; units && units[i]; i++ )
+ {
+ if( units[i]->proj_short_name &&
+ strcmp(units[i]->proj_short_name, name) == 0 ) {
+ unit_found = true;
+ to_meter = units[i]->conv_factor;
+ fr_meter = 1 / to_meter;
+ }
+ }
+ proj_unit_list_destroy(units);
+ if( !unit_found )
+ emess(1,"%s unknown unit conversion id", name);
} else
to_meter = fr_meter = 1;
geod_f = es/(1 + sqrt(1 - es));
diff --git a/src/apps/proj.cpp b/src/apps/proj.cpp
index 852cea04..0bf98b3a 100644
--- a/src/apps/proj.cpp
+++ b/src/apps/proj.cpp
@@ -380,11 +380,18 @@ int main(int argc, char **argv) {
(void)printf("%9s %-16s %-16s %s\n",
le->id, le->major, le->ell, le->name);
} else if (arg[1] == 'u') { /* list units */
- const struct PJ_UNITS *lu;
-
- for (lu = proj_list_units(); lu->id ; ++lu)
- (void)printf("%12s %-20s %s\n",
- lu->id, lu->to_meter, lu->name);
+ auto units = proj_get_units_from_database(nullptr, nullptr, "linear", false, nullptr);
+ for( int i = 0; units && units[i]; i++ )
+ {
+ if( units[i]->proj_short_name )
+ {
+ (void)printf("%12s %-20.15g %s\n",
+ units[i]->proj_short_name,
+ units[i]->conv_factor,
+ units[i]->name);
+ }
+ }
+ proj_unit_list_destroy(units);
} else
emess(1,"invalid list option: l%c",arg[1]);
exit(0);
diff --git a/src/conversions/unitconvert.cpp b/src/conversions/unitconvert.cpp
index f8439aee..609b30e0 100644
--- a/src/conversions/unitconvert.cpp
+++ b/src/conversions/unitconvert.cpp
@@ -393,9 +393,7 @@ static double get_unit_conversion_factor(const char* name,
/***********************************************************************/
int i;
const char* s;
- const PJ_UNITS *units;
-
- units = proj_list_units();
+ const PJ_UNITS *units = pj_list_linear_units();
/* Try first with linear units */
for (i = 0; (s = units[i].id) ; ++i) {
@@ -411,7 +409,7 @@ static double get_unit_conversion_factor(const char* name,
}
/* And then angular units */
- units = proj_list_angular_units();
+ units = pj_list_angular_units();
for (i = 0; (s = units[i].id) ; ++i) {
if ( strcmp(s, name) == 0 ) {
if( p_normalized_name ) {
diff --git a/src/grids.cpp b/src/grids.cpp
index b9da072c..c5fc6c74 100644
--- a/src/grids.cpp
+++ b/src/grids.cpp
@@ -309,8 +309,11 @@ constexpr uint16 TIFFTAG_GEOTRANSMATRIX = 34264;
constexpr uint16 TIFFTAG_GEOKEYDIRECTORY = 34735;
constexpr uint16 TIFFTAG_GEODOUBLEPARAMS = 34736;
constexpr uint16 TIFFTAG_GEOASCIIPARAMS = 34737;
+#ifndef TIFFTAG_GDAL_METADATA
+// Starting with libtiff > 4.1.0, those symbolic names are #define in tiff.h
constexpr uint16 TIFFTAG_GDAL_METADATA = 42112;
constexpr uint16 TIFFTAG_GDAL_NODATA = 42113;
+#endif
// ---------------------------------------------------------------------------
@@ -1373,7 +1376,8 @@ VerticalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
if (IsTIFF(header_size, header)) {
#ifdef TIFF_ENABLED
- auto set = GTiffVGridShiftSet::open(ctx, std::move(fp), actualName);
+ auto set = std::unique_ptr<VerticalShiftGridSet>(
+ GTiffVGridShiftSet::open(ctx, std::move(fp), actualName));
if (!set)
pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
return set;
@@ -2351,7 +2355,8 @@ HorizontalShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
} else if (IsTIFF(header_size,
reinterpret_cast<const unsigned char *>(header))) {
#ifdef TIFF_ENABLED
- auto set = GTiffHGridShiftSet::open(ctx, std::move(fp), actualName);
+ auto set = std::unique_ptr<HorizontalShiftGridSet>(
+ GTiffHGridShiftSet::open(ctx, std::move(fp), actualName));
if (!set)
pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
return set;
@@ -2686,8 +2691,8 @@ GenericShiftGridSet::open(PJ_CONTEXT *ctx, const std::string &filename) {
if (IsTIFF(header_size, header)) {
#ifdef TIFF_ENABLED
- auto set =
- GTiffGenericGridShiftSet::open(ctx, std::move(fp), actualName);
+ auto set = std::unique_ptr<GenericShiftGridSet>(
+ GTiffGenericGridShiftSet::open(ctx, std::move(fp), actualName));
if (!set)
pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
return set;
diff --git a/src/init.cpp b/src/init.cpp
index a25d1ccd..101fc8ad 100644
--- a/src/init.cpp
+++ b/src/init.cpp
@@ -738,7 +738,7 @@ pj_init_ctx_with_allow_init_epsg(projCtx ctx, int argc, char **argv, int allow_i
return pj_default_destructor (PIN, PJD_ERR_K_LESS_THAN_ZERO);
/* Set units */
- units = proj_list_units();
+ units = pj_list_linear_units();
s = nullptr;
if ((name = pj_param(ctx, start, "sunits").s) != nullptr) {
for (i = 0; (s = units[i].id) && strcmp(name, s) ; ++i) ;
diff --git a/src/iso19111/c_api.cpp b/src/iso19111/c_api.cpp
index d5f299c2..17fbd079 100644
--- a/src/iso19111/c_api.cpp
+++ b/src/iso19111/c_api.cpp
@@ -216,12 +216,17 @@ struct PJ_OBJ_LIST {
explicit PJ_OBJ_LIST(std::vector<IdentifiedObjectNNPtr> &&objectsIn)
: objects(std::move(objectsIn)) {}
+ virtual ~PJ_OBJ_LIST();
PJ_OBJ_LIST(const PJ_OBJ_LIST &) = delete;
PJ_OBJ_LIST &operator=(const PJ_OBJ_LIST &) = delete;
//! @endcond
};
+//! @cond Doxygen_Suppress
+PJ_OBJ_LIST::~PJ_OBJ_LIST() = default;
+//! @endcond
+
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
@@ -662,7 +667,8 @@ PJ *proj_create_from_database(PJ_CONTEXT *ctx, const char *auth_name,
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
-static const char *get_unit_category(UnitOfMeasure::Type type) {
+static const char *get_unit_category(const std::string &unit_name,
+ UnitOfMeasure::Type type) {
const char *ret = nullptr;
switch (type) {
case UnitOfMeasure::Type::UNKNOWN:
@@ -672,19 +678,26 @@ static const char *get_unit_category(UnitOfMeasure::Type type) {
ret = "none";
break;
case UnitOfMeasure::Type::ANGULAR:
- ret = "angular";
+ ret = unit_name.find(" per ") != std::string::npos ? "angular_per_time"
+ : "angular";
break;
case UnitOfMeasure::Type::LINEAR:
- ret = "linear";
+ ret = unit_name.find(" per ") != std::string::npos ? "linear_per_time"
+ : "linear";
break;
case UnitOfMeasure::Type::SCALE:
- ret = "scale";
+ ret = unit_name.find(" per year") != std::string::npos ||
+ unit_name.find(" per second") != std::string::npos
+ ? "scale_per_time"
+ : "scale";
break;
case UnitOfMeasure::Type::TIME:
ret = "time";
break;
case UnitOfMeasure::Type::PARAMETRIC:
- ret = "parametric";
+ ret = unit_name.find(" per ") != std::string::npos
+ ? "parametric_per_time"
+ : "parametric";
break;
}
return ret;
@@ -704,8 +717,9 @@ static const char *get_unit_category(UnitOfMeasure::Type type) {
* @param out_conv_factor Pointer to a value to store the conversion
* factor of the prime meridian longitude unit to radian. or NULL
* @param out_category Pointer to a string value to store the parameter name. or
- * NULL. This value might be "unknown", "none", "linear", "angular", "scale",
- * "time" or "parametric";
+ * NULL. This value might be "unknown", "none", "linear", "linear_per_time",
+ * "angular", "angular_per_time", "scale", "scale_per_time", "time",
+ * "parametric" or "parametric_per_time"
* @return TRUE in case of success
*/
int proj_uom_get_info_from_database(PJ_CONTEXT *ctx, const char *auth_name,
@@ -726,7 +740,7 @@ int proj_uom_get_info_from_database(PJ_CONTEXT *ctx, const char *auth_name,
*out_conv_factor = obj->conversionToSI();
}
if (out_category) {
- *out_category = get_unit_category(obj->type());
+ *out_category = get_unit_category(obj->name(), obj->type());
}
ctx->cpp_context->autoCloseDbIfNeeded();
return true;
@@ -2585,6 +2599,100 @@ void proj_crs_info_list_destroy(PROJ_CRS_INFO **list) {
// ---------------------------------------------------------------------------
+/** \brief Enumerate units from the database, taking into account various
+ * criteria.
+ *
+ * The returned object is an array of PROJ_UNIT_INFO* pointers, whose last
+ * entry is NULL. This array should be freed with proj_unit_list_destroy()
+ *
+ * @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 category Filter by category, if this parameter is not NULL. Category
+ * is one of "linear", "linear_per_time", "angular", "angular_per_time",
+ * "scale", "scale_per_time" or "time"
+ * @param allow_deprecated whether we should return deprecated objects as well.
+ * @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_UNIT_INFO* pointers to be freed with
+ * proj_unit_list_destroy(), or NULL in case of error.
+ *
+ * @since 7.1
+ */
+PROJ_UNIT_INFO **proj_get_units_from_database(PJ_CONTEXT *ctx,
+ const char *auth_name,
+ const char *category,
+ int allow_deprecated,
+ int *out_result_count) {
+ SANITIZE_CTX(ctx);
+ PROJ_UNIT_INFO **ret = nullptr;
+ int i = 0;
+ try {
+ auto factory = AuthorityFactory::create(getDBcontext(ctx),
+ auth_name ? auth_name : "");
+ auto list = factory->getUnitList();
+ ret = new PROJ_UNIT_INFO *[list.size() + 1];
+ for (const auto &info : list) {
+ if (category && info.category != category) {
+ continue;
+ }
+ if (!allow_deprecated && info.deprecated) {
+ continue;
+ }
+ ret[i] = new PROJ_UNIT_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]->category = pj_strdup(info.category.c_str());
+ ret[i]->conv_factor = info.convFactor;
+ ret[i]->proj_short_name =
+ info.projShortName.empty()
+ ? nullptr
+ : pj_strdup(info.projShortName.c_str());
+ ret[i]->deprecated = info.deprecated;
+ i++;
+ }
+ ret[i] = nullptr;
+ if (out_result_count)
+ *out_result_count = i;
+ ctx->cpp_context->autoCloseDbIfNeeded();
+ return ret;
+ } catch (const std::exception &e) {
+ proj_log_error(ctx, __FUNCTION__, e.what());
+ if (ret) {
+ ret[i + 1] = nullptr;
+ proj_unit_list_destroy(ret);
+ }
+ if (out_result_count)
+ *out_result_count = 0;
+ }
+ ctx->cpp_context->autoCloseDbIfNeeded();
+ return nullptr;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Destroy the result returned by
+ * proj_get_units_from_database().
+ *
+ * @since 7.1
+ */
+void proj_unit_list_destroy(PROJ_UNIT_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]->category);
+ pj_dalloc(list[i]->proj_short_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
*
@@ -6749,8 +6857,9 @@ int proj_coordoperation_get_param_index(PJ_CONTEXT *ctx,
* unit code. or NULL
* @param out_unit_category Pointer to a string value to store the parameter
* name. or
- * NULL. This value might be "unknown", "none", "linear", "angular", "scale",
- * "time" or "parametric";
+ * NULL. This value might be "unknown", "none", "linear", "linear_per_time",
+ * "angular", "angular_per_time", "scale", "scale_per_time", "time",
+ * "parametric" or "parametric_per_time"
* @return TRUE in case of success.
*/
@@ -6852,7 +6961,8 @@ int proj_coordoperation_get_param(
*out_unit_code = unit.code().c_str();
}
if (out_unit_category) {
- *out_unit_category = get_unit_category(unit.type());
+ *out_unit_category =
+ get_unit_category(unit.name(), unit.type());
}
}
}
@@ -7406,6 +7516,62 @@ void PROJ_DLL proj_operation_factory_context_set_discard_superseded(
// ---------------------------------------------------------------------------
+//! @cond Doxygen_Suppress
+/** \brief Opaque object representing a set of operation results. */
+struct PJ_OPERATION_LIST : PJ_OBJ_LIST {
+
+ PJ *source_crs;
+ PJ *target_crs;
+ bool hasPreparedOperation = false;
+ std::vector<CoordOperation> preparedOperations{};
+
+ explicit PJ_OPERATION_LIST(PJ_CONTEXT *ctx, const PJ *source_crsIn,
+ const PJ *target_crsIn,
+ std::vector<IdentifiedObjectNNPtr> &&objectsIn);
+ ~PJ_OPERATION_LIST() override;
+
+ PJ_OPERATION_LIST(const PJ_OPERATION_LIST &) = delete;
+ PJ_OPERATION_LIST &operator=(const PJ_OPERATION_LIST &) = delete;
+
+ const std::vector<CoordOperation> &getPreparedOperations(PJ_CONTEXT *ctx);
+};
+
+// ---------------------------------------------------------------------------
+
+PJ_OPERATION_LIST::PJ_OPERATION_LIST(
+ PJ_CONTEXT *ctx, const PJ *source_crsIn, const PJ *target_crsIn,
+ std::vector<IdentifiedObjectNNPtr> &&objectsIn)
+ : PJ_OBJ_LIST(std::move(objectsIn)),
+ source_crs(proj_clone(ctx, source_crsIn)),
+ target_crs(proj_clone(ctx, target_crsIn)) {}
+
+// ---------------------------------------------------------------------------
+
+PJ_OPERATION_LIST::~PJ_OPERATION_LIST() {
+ auto tmpCtxt = proj_context_create();
+ proj_assign_context(source_crs, tmpCtxt);
+ proj_assign_context(target_crs, tmpCtxt);
+ proj_destroy(source_crs);
+ proj_destroy(target_crs);
+ proj_context_destroy(tmpCtxt);
+}
+
+// ---------------------------------------------------------------------------
+
+const std::vector<CoordOperation> &
+PJ_OPERATION_LIST::getPreparedOperations(PJ_CONTEXT *ctx) {
+ if (!hasPreparedOperation) {
+ hasPreparedOperation = true;
+ preparedOperations =
+ pj_create_prepared_operations(ctx, source_crs, target_crs, this);
+ }
+ return preparedOperations;
+}
+
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
/** \brief Find a list of CoordinateOperation from source_crs to target_crs.
*
* The operations are sorted with the most relevant ones first: by
@@ -7457,7 +7623,8 @@ proj_create_operations(PJ_CONTEXT *ctx, const PJ *source_crs,
for (const auto &op : ops) {
objects.emplace_back(op);
}
- return new PJ_OBJ_LIST(std::move(objects));
+ return new PJ_OPERATION_LIST(ctx, source_crs, target_crs,
+ std::move(objects));
} catch (const std::exception &e) {
proj_log_error(ctx, __FUNCTION__, e.what());
return nullptr;
@@ -7466,6 +7633,54 @@ proj_create_operations(PJ_CONTEXT *ctx, const PJ *source_crs,
// ---------------------------------------------------------------------------
+/** Return the index of the operation that would be the most appropriate to
+ * transform the specified coordinates.
+ *
+ * This operation may use resources that are not locally available, depending
+ * on the search criteria used by proj_create_operations().
+ *
+ * This could be done by using proj_create_operations() with a punctual bounding
+ * box, but this function is faster when one needs to evaluate on many points
+ * with the same (source_crs, target_crs) tuple.
+ *
+ * @param ctx PROJ context, or NULL for default context
+ * @param operations List of operations returned by proj_create_operations()
+ * @param direction Direction into which to transform the point.
+ * @param coord Coordinate to transform
+ * @return the index in operations that would be used to transform coord. Or -1
+ * in case of error, or no match.
+ *
+ * @since 7.1
+ */
+int proj_get_suggested_operation(PJ_CONTEXT *ctx, PJ_OBJ_LIST *operations,
+ PJ_DIRECTION direction, PJ_COORD coord) {
+ SANITIZE_CTX(ctx);
+ auto opList = dynamic_cast<PJ_OPERATION_LIST *>(operations);
+ if (opList == nullptr) {
+ proj_log_error(ctx, __FUNCTION__,
+ "operations is not a list of operations");
+ return -1;
+ }
+
+ // Special case:
+ // proj_create_crs_to_crs_from_pj() always use the unique operation
+ // if there's a single one
+ if (opList->objects.size() == 1) {
+ return 0;
+ }
+
+ int iExcluded[2] = {-1, -1};
+ const auto &preparedOps = opList->getPreparedOperations(ctx);
+ int idx = pj_get_suggested_operation(ctx, preparedOps, iExcluded, direction,
+ coord);
+ if (idx >= 0) {
+ idx = preparedOps[idx].idxInOriginalList;
+ }
+ return idx;
+}
+
+// ---------------------------------------------------------------------------
+
/** \brief Return the number of objects in the result set
*
* @param result Object of type PJ_OBJ_LIST (must not be NULL)
@@ -7787,8 +8002,8 @@ PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) {
}
}
pjNew->alternativeCoordinateOperations.emplace_back(
- minxSrc, minySrc, maxxSrc, maxySrc, minxDst, minyDst,
- maxxDst, maxyDst,
+ alt.idxInOriginalList, minxSrc, minySrc, maxxSrc,
+ maxySrc, minxDst, minyDst, maxxDst, maxyDst,
pj_obj_create(ctx, co->normalizeForVisualization()),
co->nameStr(), alt.accuracy, alt.isOffshore);
}
diff --git a/src/iso19111/common.cpp b/src/iso19111/common.cpp
index f2e4de4c..e3c747b9 100644
--- a/src/iso19111/common.cpp
+++ b/src/iso19111/common.cpp
@@ -39,6 +39,7 @@
#include "proj/internal/io_internal.hpp"
#include "proj.h"
+#include "proj_internal.h"
#include <cmath> // M_PI
#include <cstdlib>
@@ -312,7 +313,7 @@ bool UnitOfMeasure::operator!=(const UnitOfMeasure &other) PROJ_PURE_DEFN {
//! @cond Doxygen_Suppress
std::string UnitOfMeasure::exportToPROJString() const {
if (type() == Type::LINEAR) {
- auto proj_units = proj_list_units();
+ auto proj_units = pj_list_linear_units();
for (int i = 0; proj_units[i].id != nullptr; i++) {
if (::fabs(proj_units[i].factor - conversionToSI()) <
1e-10 * conversionToSI()) {
@@ -320,7 +321,7 @@ std::string UnitOfMeasure::exportToPROJString() const {
}
}
} else if (type() == Type::ANGULAR) {
- auto proj_angular_units = proj_list_angular_units();
+ auto proj_angular_units = pj_list_angular_units();
for (int i = 0; proj_angular_units[i].id != nullptr; i++) {
if (::fabs(proj_angular_units[i].factor - conversionToSI()) <
1e-10 * conversionToSI()) {
diff --git a/src/iso19111/coordinateoperation.cpp b/src/iso19111/coordinateoperation.cpp
index 63b11761..33210d24 100644
--- a/src/iso19111/coordinateoperation.cpp
+++ b/src/iso19111/coordinateoperation.cpp
@@ -2443,7 +2443,7 @@ bool ParameterValue::_isEquivalentTo(const util::IComparable *other,
}
switch (type()) {
case Type::MEASURE: {
- return value()._isEquivalentTo(otherPV->value(), criterion);
+ return value()._isEquivalentTo(otherPV->value(), criterion, 2e-10);
}
case Type::STRING:
diff --git a/src/iso19111/crs.cpp b/src/iso19111/crs.cpp
index e1181d43..b0f5bcf9 100644
--- a/src/iso19111/crs.cpp
+++ b/src/iso19111/crs.cpp
@@ -310,17 +310,15 @@ CRSNNPtr CRS::alterCSLinearUnit(const common::UnitOfMeasure &unit) const {
auto cartCS = util::nn_dynamic_pointer_cast<cs::CartesianCS>(
engCRS->coordinateSystem());
if (cartCS) {
- auto props = createPropertyMap(this);
- props.set("FORCE_OUTPUT_CS", true);
- return EngineeringCRS::create(props, engCRS->datum(),
+ return EngineeringCRS::create(createPropertyMap(this),
+ engCRS->datum(),
cartCS->alterUnit(unit));
} else {
auto vertCS = util::nn_dynamic_pointer_cast<cs::VerticalCS>(
engCRS->coordinateSystem());
if (vertCS) {
- auto props = createPropertyMap(this);
- props.set("FORCE_OUTPUT_CS", true);
- return EngineeringCRS::create(props, engCRS->datum(),
+ return EngineeringCRS::create(createPropertyMap(this),
+ engCRS->datum(),
vertCS->alterUnit(unit));
}
}
@@ -3743,7 +3741,7 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
return res;
}
res.emplace_back(crsNN, eqName ? 90 : 70);
- } else if (crs->nameStr() == thisName &&
+ } else if (objects.size() == 1 &&
CRS::getPrivate()->implicitCS_ &&
coordinateSystem()
->axisList()[0]
@@ -3762,10 +3760,11 @@ ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
derivingConversionRef()->_isEquivalentTo(
crs->derivingConversionRef().get(),
util::IComparable::Criterion::EQUIVALENT,
- dbContext) &&
- objects.size() == 1) {
+ dbContext)) {
res.clear();
- res.emplace_back(crsNN, 100);
+ res.emplace_back(crsNN, crs->nameStr() == thisName
+ ? 100
+ : eqName ? 90 : 70);
return res;
} else {
res.emplace_back(crsNN, 25);
@@ -5433,9 +5432,7 @@ bool TemporalCRS::_isEquivalentTo(
// ---------------------------------------------------------------------------
//! @cond Doxygen_Suppress
-struct EngineeringCRS::Private {
- bool forceOutputCS_ = false;
-};
+struct EngineeringCRS::Private {};
//! @endcond
// ---------------------------------------------------------------------------
@@ -5493,17 +5490,6 @@ EngineeringCRS::create(const util::PropertyMap &properties,
crs->assignSelf(crs);
crs->setProperties(properties);
- const auto pVal = properties.get("FORCE_OUTPUT_CS");
- if (pVal) {
- if (const auto genVal =
- dynamic_cast<const util::BoxedValue *>(pVal->get())) {
- if (genVal->type() == util::BoxedValue::Type::BOOLEAN &&
- genVal->booleanValue()) {
- crs->d->forceOutputCS_ = true;
- }
- }
- }
-
return crs;
}
@@ -5518,11 +5504,16 @@ void EngineeringCRS::_exportToWKT(io::WKTFormatter *formatter) const {
formatter->addQuotedString(nameStr());
if (isWKT2 || !datum()->nameStr().empty()) {
datum()->_exportToWKT(formatter);
- coordinateSystem()->_exportToWKT(formatter);
}
- if (!isWKT2 && d->forceOutputCS_) {
+ if (!isWKT2) {
coordinateSystem()->axisList()[0]->unit()._exportToWKT(formatter);
}
+
+ const auto oldAxisOutputRule = formatter->outputAxis();
+ formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES);
+ coordinateSystem()->_exportToWKT(formatter);
+ formatter->setOutputAxis(oldAxisOutputRule);
+
ObjectUsage::baseExportToWKT(formatter);
formatter->endNode();
}
diff --git a/src/iso19111/factory.cpp b/src/iso19111/factory.cpp
index 1c836367..6a707e0d 100644
--- a/src/iso19111/factory.cpp
+++ b/src/iso19111/factory.cpp
@@ -5214,6 +5214,64 @@ std::list<AuthorityFactory::CRSInfo> AuthorityFactory::getCRSInfoList() const {
// ---------------------------------------------------------------------------
+//! @cond Doxygen_Suppress
+AuthorityFactory::UnitInfo::UnitInfo()
+ : authName{}, code{}, name{}, category{}, convFactor{}, projShortName{},
+ deprecated{} {}
+//! @endcond
+
+// ---------------------------------------------------------------------------
+
+/** \brief Return the list of units.
+ * @throw FactoryException
+ *
+ * @since 7.1
+ */
+std::list<AuthorityFactory::UnitInfo> AuthorityFactory::getUnitList() const {
+ std::string sql = "SELECT auth_name, code, name, type, conv_factor, "
+ "proj_short_name, deprecated FROM unit_of_measure";
+ ListOfParams params;
+ if (d->hasAuthorityRestriction()) {
+ sql += " WHERE auth_name = ?";
+ params.emplace_back(d->authority());
+ }
+ sql += " ORDER BY auth_name, code";
+
+ auto sqlRes = d->run(sql, params);
+ std::list<AuthorityFactory::UnitInfo> res;
+ for (const auto &row : sqlRes) {
+ AuthorityFactory::UnitInfo info;
+ info.authName = row[0];
+ info.code = row[1];
+ info.name = row[2];
+ const std::string &raw_category(row[3]);
+ if (raw_category == "length") {
+ info.category = info.name.find(" per ") != std::string::npos
+ ? "linear_per_time"
+ : "linear";
+ } else if (raw_category == "angle") {
+ info.category = info.name.find(" per ") != std::string::npos
+ ? "angular_per_time"
+ : "angular";
+ } else if (raw_category == "scale") {
+ info.category =
+ info.name.find(" per year") != std::string::npos ||
+ info.name.find(" per second") != std::string::npos
+ ? "scale_per_time"
+ : "scale";
+ } else {
+ info.category = raw_category;
+ }
+ info.convFactor = row[4].empty() ? 0 : c_locale_stod(row[4]);
+ info.projShortName = row[5];
+ info.deprecated = row[6] == "1";
+ res.emplace_back(info);
+ }
+ return res;
+}
+
+// ---------------------------------------------------------------------------
+
/** \brief Gets the official name from a possibly alias name.
*
* @param aliasedName Alias name.
diff --git a/src/iso19111/io.cpp b/src/iso19111/io.cpp
index c57014d9..6a1d7e32 100644
--- a/src/iso19111/io.cpp
+++ b/src/iso19111/io.cpp
@@ -3689,6 +3689,15 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
auto props = buildProperties(node);
+ auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
+ const auto &nodeValue = nodeP->value();
+ if (isNull(csNode) && !ci_equal(nodeValue, WKTConstants::PROJCS) &&
+ !ci_equal(nodeValue, WKTConstants::BASEPROJCRS)) {
+ ThrowMissing(WKTConstants::CS_);
+ }
+ auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
+ auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
+
const std::string projCRSName = stripQuotes(nodeP->children()[0]);
if (esriStyle_ && dbContext_) {
// It is likely that the ESRI definition of EPSG:32661 (UPS North) &
@@ -3709,6 +3718,27 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
projCRSName, "projected_crs", "ESRI", false, outTableName,
authNameFromAlias, codeFromAlias);
if (!officialName.empty()) {
+ // Special case for https://github.com/OSGeo/PROJ/issues/2086
+ // The name of the CRS to identify is
+ // NAD_1983_HARN_StatePlane_Colorado_North_FIPS_0501
+ // whereas it should be
+ // NAD_1983_HARN_StatePlane_Colorado_North_FIPS_0501_Feet
+ constexpr double US_FOOT_CONV_FACTOR = 12.0 / 39.37;
+ if (projCRSName.find("_FIPS_") != std::string::npos &&
+ projCRSName.find("_Feet") == std::string::npos &&
+ std::fabs(
+ cartesianCS->axisList()[0]->unit().conversionToSI() -
+ US_FOOT_CONV_FACTOR) < 1e-10 * US_FOOT_CONV_FACTOR) {
+ auto officialNameFromFeet =
+ authFactory->getOfficialNameFromAlias(
+ projCRSName + "_Feet", "projected_crs", "ESRI",
+ false, outTableName, authNameFromAlias,
+ codeFromAlias);
+ if (!officialNameFromFeet.empty()) {
+ officialName = officialNameFromFeet;
+ }
+ }
+
props.set(IdentifiedObject::NAME_KEY, officialName);
}
}
@@ -3739,15 +3769,6 @@ WKTParser::Private::buildProjectedCRS(const WKTNodeNNPtr &node) {
? buildConversion(conversionNode, linearUnit, angularUnit)
: buildProjection(node, projectionNode, linearUnit, angularUnit);
- auto &csNode = nodeP->lookForChild(WKTConstants::CS_);
- const auto &nodeValue = nodeP->value();
- if (isNull(csNode) && !ci_equal(nodeValue, WKTConstants::PROJCS) &&
- !ci_equal(nodeValue, WKTConstants::BASEPROJCRS)) {
- ThrowMissing(WKTConstants::CS_);
- }
- auto cs = buildCS(csNode, node, UnitOfMeasure::NONE);
- auto cartesianCS = nn_dynamic_pointer_cast<CartesianCS>(cs);
-
// No explicit AXIS node ? (WKT1)
if (isNull(nodeP->lookForChild(WKTConstants::AXIS))) {
props.set("IMPLICIT_CS", true);
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index d08a0ddf..8978e20f 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -292,7 +292,7 @@ source_group("Source Files\\Transformations"
source_group("Source Files\\ISO19111"
FILES ${SRC_LIBPROJ_ISO19111})
-include_directories(${CMAKE_SOURCE_DIR}/include)
+include_directories(${PROJ_SOURCE_DIR}/include)
include_directories(${CMAKE_CURRENT_BINARY_DIR})
source_group("CMake Files" FILES CMakeLists.txt)
diff --git a/src/proj.h b/src/proj.h
index 90a11739..d2464d69 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -151,6 +151,24 @@ extern "C" {
#endif
#endif
+#ifdef PROJ_SUPPRESS_DEPRECATION_MESSAGE
+ #define PROJ_DEPRECATED(decl, msg) decl
+#elif defined(__has_extension)
+ #if __has_extension(attribute_deprecated_with_message)
+ #define PROJ_DEPRECATED(decl, msg) decl __attribute__ ((deprecated(msg)))
+ #elif defined(__GNUC__)
+ #define PROJ_DEPRECATED(decl, msg) decl __attribute__ ((deprecated))
+ #else
+ #define PROJ_DEPRECATED(decl, msg) decl
+ #endif
+#elif defined(__GNUC__)
+ #define PROJ_DEPRECATED(decl, msg) decl __attribute__ ((deprecated))
+#elif defined(_MSVC_VER)
+ #define PROJ_DEPRECATED(decl, msg) __declspec(deprecated(msg)) decl
+#else
+ #define PROJ_DEPRECATED(decl, msg) decl
+#endif
+
/* The version numbers should be updated with every release! **/
#define PROJ_VERSION_MAJOR 7
#define PROJ_VERSION_MINOR 1
@@ -608,8 +626,8 @@ PJ_INIT_INFO PROJ_DLL proj_init_info(const char *initname);
/* Get lists of operations, ellipsoids, units and prime meridians. */
const PJ_OPERATIONS PROJ_DLL *proj_list_operations(void);
const PJ_ELLPS PROJ_DLL *proj_list_ellps(void);
-const PJ_UNITS PROJ_DLL *proj_list_units(void);
-const PJ_UNITS PROJ_DLL *proj_list_angular_units(void);
+PROJ_DEPRECATED(const PJ_UNITS PROJ_DLL *proj_list_units(void), "Deprecated by proj_get_units_from_database");
+PROJ_DEPRECATED(const PJ_UNITS PROJ_DLL *proj_list_angular_units(void), "Deprecated by proj_get_units_from_database");
const PJ_PRIME_MERIDIANS PROJ_DLL *proj_list_prime_meridians(void);
/* These are trivial, and while occasionally useful in real code, primarily here to */
@@ -912,6 +930,39 @@ typedef struct
int allow_deprecated;
} PROJ_CRS_LIST_PARAMETERS;
+/** \brief Structure given description of a unit.
+ *
+ * This structure may grow over time, and should not be directly allocated by
+ * client code.
+ * @since 7.1
+ */
+typedef struct
+{
+ /** Authority name. */
+ char* auth_name;
+
+ /** Object code. */
+ char* code;
+
+ /** Object name. For example "metre", "US survey foot", etc. */
+ char* name;
+
+ /** Category of the unit: one of "linear", "linear_per_time", "angular",
+ * "angular_per_time", "scale", "scale_per_time" or "time" */
+ char* category;
+
+ /** Conversion factor to apply to transform from that unit to the
+ * corresponding SI unit (metre for "linear", radian for "angular", etc.).
+ * It might be 0 in some cases to indicate no known conversion factor. */
+ double conv_factor;
+
+ /** PROJ short name, like "m", "ft", "us-ft", etc... Might be NULL */
+ char* proj_short_name;
+
+ /** Whether the object is deprecated */
+ int deprecated;
+} PROJ_UNIT_INFO;
+
/**@}*/
@@ -1077,6 +1128,15 @@ PROJ_CRS_INFO PROJ_DLL **proj_get_crs_info_list_from_database(
void PROJ_DLL proj_crs_info_list_destroy(PROJ_CRS_INFO** list);
+PROJ_UNIT_INFO PROJ_DLL **proj_get_units_from_database(
+ PJ_CONTEXT *ctx,
+ const char *auth_name,
+ const char *category,
+ int allow_deprecated,
+ int *out_result_count);
+
+void PROJ_DLL proj_unit_list_destroy(PROJ_UNIT_INFO** list);
+
/* ------------------------------------------------------------------------- */
@@ -1157,6 +1217,11 @@ PJ PROJ_DLL *proj_list_get(PJ_CONTEXT *ctx,
void PROJ_DLL proj_list_destroy(PJ_OBJ_LIST *result);
+int PROJ_DLL proj_get_suggested_operation(PJ_CONTEXT *ctx,
+ PJ_OBJ_LIST *operations,
+ PJ_DIRECTION direction,
+ PJ_COORD coord);
+
/* ------------------------------------------------------------------------- */
PJ PROJ_DLL *proj_crs_get_geodetic_crs(PJ_CONTEXT *ctx, const PJ *crs);
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 8f73200d..a3bc9f17 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -287,6 +287,55 @@ typedef PJ_COORD (* PJ_OPERATOR) (PJ_COORD, PJ *);
#define PJD_GRIDSHIFT 3
#define PJD_WGS84 4 /* WGS84 (or anything considered equivalent) */
+struct CoordOperation
+{
+ int idxInOriginalList;
+ double minxSrc = 0.0;
+ double minySrc = 0.0;
+ double maxxSrc = 0.0;
+ double maxySrc = 0.0;
+ double minxDst = 0.0;
+ double minyDst = 0.0;
+ double maxxDst = 0.0;
+ double maxyDst = 0.0;
+ PJ* pj = nullptr;
+ std::string name{};
+ double accuracy = -1.0;
+ bool isOffshore = false;
+
+ CoordOperation(int idxInOriginalListIn,
+ double minxSrcIn, double minySrcIn, double maxxSrcIn, double maxySrcIn,
+ double minxDstIn, double minyDstIn, double maxxDstIn, double maxyDstIn,
+ PJ* pjIn, const std::string& nameIn, double accuracyIn, bool isOffshoreIn):
+ idxInOriginalList(idxInOriginalListIn),
+ minxSrc(minxSrcIn), minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn),
+ minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn), maxyDst(maxyDstIn),
+ pj(pjIn), name(nameIn),
+ accuracy(accuracyIn),
+ isOffshore(isOffshoreIn)
+ {
+ }
+
+ CoordOperation(const CoordOperation&) = delete;
+
+ CoordOperation(CoordOperation&& other):
+ idxInOriginalList(other.idxInOriginalList),
+ minxSrc(other.minxSrc), minySrc(other.minySrc), maxxSrc(other.maxxSrc), maxySrc(other.maxySrc),
+ minxDst(other.minxDst), minyDst(other.minyDst), maxxDst(other.maxxDst), maxyDst(other.maxyDst),
+ name(std::move(other.name)),
+ accuracy(other.accuracy),
+ isOffshore(other.isOffshore) {
+ pj = other.pj;
+ other.pj = nullptr;
+ }
+
+ CoordOperation& operator=(const CoordOperation&) = delete;
+
+ ~CoordOperation()
+ {
+ proj_destroy(pj);
+ }
+};
/* base projection data structure */
struct PJconsts {
@@ -493,52 +542,6 @@ struct PJconsts {
/*************************************************************************************
proj_create_crs_to_crs() alternative coordinate operations
**************************************************************************************/
-
- struct CoordOperation
- {
- double minxSrc = 0.0;
- double minySrc = 0.0;
- double maxxSrc = 0.0;
- double maxySrc = 0.0;
- double minxDst = 0.0;
- double minyDst = 0.0;
- double maxxDst = 0.0;
- double maxyDst = 0.0;
- PJ* pj = nullptr;
- std::string name{};
- double accuracy = -1.0;
- bool isOffshore = false;
-
- CoordOperation(double minxSrcIn, double minySrcIn, double maxxSrcIn, double maxySrcIn,
- double minxDstIn, double minyDstIn, double maxxDstIn, double maxyDstIn,
- PJ* pjIn, const std::string& nameIn, double accuracyIn, bool isOffshoreIn):
- minxSrc(minxSrcIn), minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn),
- minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn), maxyDst(maxyDstIn),
- pj(pjIn), name(nameIn),
- accuracy(accuracyIn),
- isOffshore(isOffshoreIn)
- {
- }
-
- CoordOperation(const CoordOperation&) = delete;
-
- CoordOperation(CoordOperation&& other):
- minxSrc(other.minxSrc), minySrc(other.minySrc), maxxSrc(other.maxxSrc), maxySrc(other.maxySrc),
- minxDst(other.minxDst), minyDst(other.minyDst), maxxDst(other.maxxDst), maxyDst(other.maxyDst),
- name(std::move(other.name)),
- accuracy(other.accuracy),
- isOffshore(other.isOffshore) {
- pj = other.pj;
- other.pj = nullptr;
- }
-
- CoordOperation& operator=(const CoordOperation&) = delete;
-
- ~CoordOperation()
- {
- proj_destroy(pj);
- }
- };
std::vector<CoordOperation> alternativeCoordinateOperations{};
int iCurCoordOp = -1;
@@ -873,6 +876,20 @@ std::string PROJ_DLL pj_context_get_user_writable_directory(PJ_CONTEXT *ctx, boo
void PROJ_DLL pj_context_set_user_writable_directory(PJ_CONTEXT* ctx, const std::string& path);
std::string PROJ_DLL pj_get_relative_share_proj(PJ_CONTEXT *ctx);
+std::vector<CoordOperation> pj_create_prepared_operations(PJ_CONTEXT *ctx,
+ const PJ *source_crs,
+ const PJ *target_crs,
+ PJ_OBJ_LIST* op_list);
+
+int pj_get_suggested_operation(PJ_CONTEXT *ctx,
+ const std::vector<CoordOperation>& opList,
+ const int iExcluded[2],
+ PJ_DIRECTION direction,
+ PJ_COORD coord);
+
+const PJ_UNITS *pj_list_linear_units();
+const PJ_UNITS *pj_list_angular_units();
+
/* classic public API */
#include "proj_api.h"
diff --git a/src/projections/stere.cpp b/src/projections/stere.cpp
index d95bb7fa..abc4aa13 100644
--- a/src/projections/stere.cpp
+++ b/src/projections/stere.cpp
@@ -5,7 +5,7 @@
#include <math.h>
PROJ_HEAD(stere, "Stereographic") "\n\tAzi, Sph&Ell\n\tlat_ts=";
-PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Sph&Ell\n\tsouth";
+PROJ_HEAD(ups, "Universal Polar Stereographic") "\n\tAzi, Ell\n\tsouth";
namespace { // anonymous namespace
@@ -320,8 +320,7 @@ PJ *PROJECTION(ups) {
/* International Ellipsoid */
P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI;
if (P->es == 0.0) {
- proj_errno_set(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
- return pj_default_destructor (P, ENOMEM);
+ return pj_default_destructor (P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
}
P->k0 = .994;
P->x0 = 2000000.;
diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp
index 1414e542..75a77ccd 100644
--- a/src/projections/tmerc.cpp
+++ b/src/projections/tmerc.cpp
@@ -22,7 +22,7 @@
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";
+PROJ_HEAD(utm, "Universal Transverse Mercator (UTM)") "\n\tCyl, Ell\n\tzone= south approx";
namespace { // anonymous namespace
struct pj_opaque_approx {
@@ -345,42 +345,31 @@ static PJ *setup_approx(PJ *P) {
/*****************************************************************************/
/* 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;
+inline
+static double gatg(const double *p1, int len_p1, double B, double cos_2B, double sin_2B) {
+ double h = 0, h1, h2 = 0;
- cos_2B = 2*cos(2*B);
- p = p1 + len_p1;
+ const double two_cos_2B = 2*cos_2B;
+ const double* p = p1 + len_p1;
h1 = *--p;
while (p - p1) {
- h = -h2 + cos_2B*h1 + *--p;
+ h = -h2 + two_cos_2B*h1 + *--p;
h2 = h1;
h1 = h;
}
- return (B + h*sin(2*B));
+ return (B + h*sin_2B);
}
/* 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;
+inline
+static double clenS(const double *a, int size,
+ double sin_arg_r, double cos_arg_r,
+ double sinh_arg_i, double cosh_arg_i,
+ double *R, double *I) {
+ double r, i, hr, hr1, hr2, hi, hi1, hi2;
/* 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);
+ const double* p = a + size;
r = 2*cos_arg_r*cosh_arg_i;
i = -2*sin_arg_r*sinh_arg_i;
@@ -427,28 +416,70 @@ static double clens(double *a, int size, double arg_r) {
static PJ_XY exact_e_fwd (PJ_LP lp, PJ *P) {
PJ_XY xy = {0.0,0.0};
struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
- double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
- double Cn = lp.phi, Ce = lp.lam;
/* ell. LAT, LNG -> Gaussian LAT, LNG */
- Cn = gatg (Q->cbg, PROJ_ETMERC_ORDER, Cn);
+ double Cn = gatg (Q->cbg, PROJ_ETMERC_ORDER, lp.phi, cos(2*lp.phi), sin(2*lp.phi));
/* 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);
+ const double sin_Cn = sin (Cn);
+ const double cos_Cn = cos (Cn);
+ const double sin_Ce = sin (lp.lam);
+ const double cos_Ce = cos (lp.lam);
+
+ const double cos_Cn_cos_Ce = cos_Cn*cos_Ce;
+ Cn = atan2 (sin_Cn, cos_Cn_cos_Ce);
+
+ const double inv_denom_tan_Ce = 1. / hypot (sin_Cn, cos_Cn_cos_Ce);
+ const double tan_Ce = sin_Ce*cos_Cn * inv_denom_tan_Ce;
+#if 0
+ // Variant of the above: found not to be measurably faster
+ const double sin_Ce_cos_Cn = sin_Ce*cos_Cn;
+ const double denom = sqrt(1 - sin_Ce_cos_Cn * sin_Ce_cos_Cn);
+ const double tan_Ce = sin_Ce_cos_Cn / denom;
#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);
+ double Ce = asinh ( tan_Ce ); /* Replaces: Ce = log(tan(FORTPI + Ce*0.5)); */
+
+/*
+ * Non-optimized version:
+ * const double sin_arg_r = sin(2*Cn);
+ * const double cos_arg_r = cos(2*Cn);
+ *
+ * Given:
+ * sin(2 * Cn) = 2 sin(Cn) cos(Cn)
+ * sin(atan(y)) = y / sqrt(1 + y^2)
+ * cos(atan(y)) = 1 / sqrt(1 + y^2)
+ * ==> sin(2 * Cn) = 2 tan_Cn / (1 + tan_Cn^2)
+ *
+ * cos(2 * Cn) = 2cos^2(Cn) - 1
+ * = 2 / (1 + tan_Cn^2) - 1
+ */
+ const double tmp_r = 2 * cos_Cn_cos_Ce * inv_denom_tan_Ce * inv_denom_tan_Ce;
+ const double sin_arg_r = sin_Cn * tmp_r;
+ const double cos_arg_r = cos_Cn_cos_Ce * tmp_r - 1;
+
+/*
+ * Non-optimized version:
+ * const double sinh_arg_i = sinh(2*Ce);
+ * const double cosh_arg_i = cosh(2*Ce);
+ *
+ * Given
+ * sinh(2 * Ce) = 2 sinh(Ce) cosh(Ce)
+ * sinh(asinh(y)) = y
+ * cosh(asinh(y)) = sqrt(1 + y^2)
+ * ==> sinh(2 * Ce) = 2 tan_Ce sqrt(1 + tan_Ce^2)
+ *
+ * cosh(2 * Ce) = 2cosh^2(Ce) - 1
+ * = 2 * (1 + tan_Ce^2) - 1
+ */
+ const double tmp_i = 1 + tan_Ce * tan_Ce;
+ const double sinh_arg_i = 2 * tan_Ce * sqrt(tmp_i);
+ const double cosh_arg_i = 2 * tmp_i - 1;
+
+ double dCn, dCe;
+ Cn += clenS (Q->gtu, PROJ_ETMERC_ORDER,
+ sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i,
+ &dCn, &dCe);
Ce += dCe;
if (fabs (Ce) <= 2.623395162778) {
xy.y = Q->Qn * Cn + Q->Zb; /* Northing */
@@ -463,32 +494,68 @@ static PJ_XY exact_e_fwd (PJ_LP lp, PJ *P) {
static PJ_LP exact_e_inv (PJ_XY xy, PJ *P) {
PJ_LP lp = {0.0,0.0};
struct pj_opaque_exact *Q = static_cast<struct pj_opaque_exact*>(P->opaque);
- double sin_Cn, cos_Cn, cos_Ce, sin_Ce, dCn, dCe;
- double Cn = xy.y, Ce = xy.x;
/* normalize N, E */
- Cn = (Cn - Q->Zb)/Q->Qn;
- Ce = Ce/Q->Qn;
+ double Cn = (xy.y - Q->Zb)/Q->Qn;
+ double Ce = xy.x/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);
+ const double sin_arg_r = sin(2*Cn);
+ const double cos_arg_r = cos(2*Cn);
+
+ //const double sinh_arg_i = sinh(2*Ce);
+ //const double cosh_arg_i = cosh(2*Ce);
+ const double exp_2_Ce = exp(2*Ce);
+ const double half_inv_exp_2_Ce = 0.5 / exp_2_Ce;
+ const double sinh_arg_i = 0.5 * exp_2_Ce - half_inv_exp_2_Ce;
+ const double cosh_arg_i = 0.5 * exp_2_Ce + half_inv_exp_2_Ce;
+
+ double dCn_ignored, dCe;
+ Cn += clenS(Q->utg, PROJ_ETMERC_ORDER,
+ sin_arg_r, cos_arg_r, sinh_arg_i, cosh_arg_i,
+ &dCn_ignored, &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);
+ const double sin_Cn = sin (Cn);
+ const double cos_Cn = cos (Cn);
+
+#if 0
+ // Non-optimized version:
+ double sin_Ce, cos_Ce;
+ Ce = atan (sinh (Ce)); // Replaces: Ce = 2*(atan(exp(Ce)) - FORTPI);
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));
+#else
+/*
+ * One can divide both member of Ce = atan2(...) by cos_Ce, which gives:
+ * Ce = atan2 (tan_Ce, cos_Cn) = atan2(sinh(Ce), cos_Cn)
+ *
+ * and the same for Cn = atan2(...)
+ * Cn = atan2 (sin_Cn, hypot (sin_Ce, cos_Ce*cos_Cn)/cos_Ce)
+ * = atan2 (sin_Cn, hypot (sin_Ce/cos_Ce, cos_Cn))
+ * = atan2 (sin_Cn, hypot (tan_Ce, cos_Cn))
+ * = atan2 (sin_Cn, hypot (sinhCe, cos_Cn))
+ */
+ const double sinhCe = sinh (Ce);
+ Ce = atan2 (sinhCe, cos_Cn);
+ const double modulus_Ce = hypot (sinhCe, cos_Cn);
+ Cn = atan2 (sin_Cn, modulus_Ce);
+#endif
+
/* Gaussian LAT, LNG -> ell. LAT, LNG */
- lp.phi = gatg (Q->cgb, PROJ_ETMERC_ORDER, Cn);
+
+ // Optimization of the computation of cos(2*Cn) and sin(2*Cn)
+ const double tmp = 2 * modulus_Ce / (sinhCe * sinhCe + 1);
+ const double sin_2_Cn = sin_Cn * tmp;
+ const double cos_2_Cn = tmp * modulus_Ce - 1.;
+ //const double cos_2_Cn = cos(2 * Cn);
+ //const double sin_2_Cn = sin(2 * Cn);
+
+ lp.phi = gatg (Q->cgb, PROJ_ETMERC_ORDER, Cn, cos_2_Cn, sin_2_Cn);
lp.lam = Ce;
}
else
@@ -573,7 +640,7 @@ static PJ *setup_exact(PJ *P) {
Q->gtu[5] = np*(212378941/319334400.0);
/* Gaussian latitude value of the origin latitude */
- Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0);
+ Z = gatg (Q->cbg, PROJ_ETMERC_ORDER, P->phi0, cos(2*P->phi0), sin(2*P->phi0));
/* Origin northing minus true northing at the origin latitude */
/* i.e. true northing = N - P->Zb */
@@ -627,8 +694,7 @@ PJ *PROJECTION(etmerc) {
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);
+ return pj_default_destructor(P, PJD_ERR_ELLIPSOID_USE_REQUIRED);
}
if (P->lam0 < -1000.0 || P->lam0 > 1000.0) {
return pj_default_destructor(P, PJD_ERR_INVALID_UTM_ZONE);
diff --git a/src/units.cpp b/src/units.cpp
index 34a71db1..53f43163 100644
--- a/src/units.cpp
+++ b/src/units.cpp
@@ -36,6 +36,12 @@ pj_units[] = {
{nullptr, nullptr, nullptr, 0.0}
};
+// For internal use
+const PJ_UNITS *pj_list_linear_units()
+{
+ return pj_units;
+}
+
const PJ_UNITS *proj_list_units()
{
return pj_units;
@@ -52,6 +58,12 @@ pj_angular_units[] = {
{nullptr, nullptr, nullptr, 0.0}
};
+// For internal use
+const PJ_UNITS *pj_list_angular_units()
+{
+ return pj_angular_units;
+}
+
const PJ_UNITS *proj_list_angular_units()
{
return pj_angular_units;
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-15 03:54:35 +0000