aboutsummaryrefslogtreecommitdiff
path: root/src/cs2cs.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/cs2cs.cpp')
-rw-r--r--src/cs2cs.cpp327
1 files changed, 250 insertions, 77 deletions
diff --git a/src/cs2cs.cpp b/src/cs2cs.cpp
index 6f302ae3..6f4c4a55 100644
--- a/src/cs2cs.cpp
+++ b/src/cs2cs.cpp
@@ -26,6 +26,8 @@
* DEALINGS IN THE SOFTWARE.
*****************************************************************************/
+#define FROM_PROJ_CPP
+
#include <ctype.h>
#include <locale.h>
#include <math.h>
@@ -33,6 +35,11 @@
#include <stdlib.h>
#include <string.h>
+#include <cassert>
+#include <string>
+
+#include <proj/internal/internal.hpp>
+
// PROJ include order is sensitive
// clang-format off
#include "proj.h"
@@ -41,9 +48,15 @@
// clang-format on
#define MAX_LINE 1000
-#define MAX_PARGS 100
-static projPJ fromProj, toProj;
+static PJ *transformation = nullptr;
+
+static bool srcIsGeog = false;
+static double srcToRadians = 0.0;
+
+static bool destIsGeog = false;
+static double destToRadians = 0.0;
+static bool destIsLatLong = false;
static int reversein = 0, /* != 0 reverse input arguments */
reverseout = 0, /* != 0 reverse output arguments */
@@ -116,19 +129,45 @@ static void process(FILE *fid)
}
if (data.u != HUGE_VAL) {
- if (pj_transform(fromProj, toProj, 1, 0, &(data.u), &(data.v),
- &z) != 0) {
- data.u = HUGE_VAL;
- data.v = HUGE_VAL;
- emess(-3, "pj_transform(): %s", pj_strerrno(pj_errno));
+
+ if (srcIsGeog) {
+ /* dmstor gives values to radians. Convert now to the SRS unit
+ */
+ data.u /= srcToRadians;
+ data.v /= srcToRadians;
}
+
+ PJ_COORD coord;
+ coord.xyzt.x = data.u;
+ coord.xyzt.y = data.v;
+ coord.xyzt.z = z;
+ coord.xyzt.t = HUGE_VAL;
+ coord = proj_trans(transformation, PJ_FWD, coord);
+ data.u = coord.xyz.x;
+ data.v = coord.xyz.y;
+ z = coord.xyz.z;
}
if (data.u == HUGE_VAL) /* error output */
fputs(oterr, stdout);
- else if (pj_is_latlong(toProj) && !oform) { /*ascii DMS output */
- if (reverseout) {
+ else if (destIsGeog && !oform) { /*ascii DMS output */
+
+ // rtodms() expect radians: convert from the output SRS unit
+ data.u *= destToRadians;
+ data.v *= destToRadians;
+
+ if (destIsLatLong) {
+ if (reverseout) {
+ fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+ putchar('\t');
+ fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+ } else {
+ fputs(rtodms(pline, data.u, 'N', 'S'), stdout);
+ putchar('\t');
+ fputs(rtodms(pline, data.v, 'E', 'W'), stdout);
+ }
+ } else if (reverseout) {
fputs(rtodms(pline, data.v, 'N', 'S'), stdout);
putchar('\t');
fputs(rtodms(pline, data.u, 'E', 'W'), stdout);
@@ -139,9 +178,9 @@ static void process(FILE *fid)
}
} else { /* x-y or decimal degree ascii output */
- if (proj_angular_output(toProj, PJ_FWD)) {
- data.v *= RAD_TO_DEG;
- data.u *= RAD_TO_DEG;
+ if (destIsGeog) {
+ data.v *= destToRadians * RAD_TO_DEG;
+ data.u *= destToRadians * RAD_TO_DEG;
}
if (reverseout) {
printf(oform, data.v);
@@ -167,17 +206,116 @@ static void process(FILE *fid)
}
/************************************************************************/
+/* instanciate_crs() */
+/************************************************************************/
+
+static PJ_OBJ *instanciate_crs(const std::string &definition,
+ const char *const *optionsImportCRS,
+ bool &isGeog, double &toRadians,
+ bool &isLatFirst) {
+ PJ_OBJ *crs = proj_obj_create_from_user_input(nullptr, definition.c_str(),
+ optionsImportCRS);
+ if (!crs) {
+ return nullptr;
+ }
+
+ isGeog = false;
+ toRadians = 0.0;
+ isLatFirst = false;
+
+ auto type = proj_obj_get_type(crs);
+ if (type == PJ_OBJ_TYPE_BOUND_CRS) {
+ auto base = proj_obj_get_source_crs(crs);
+ proj_obj_unref(crs);
+ crs = base;
+ type = proj_obj_get_type(crs);
+ }
+ if (type == PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS ||
+ type == PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+ auto cs = proj_obj_crs_get_coordinate_system(crs);
+ assert(cs);
+
+ isGeog = true;
+ const char *axisName = "";
+ proj_obj_cs_get_axis_info(cs, 0,
+ &axisName, // name,
+ nullptr, // abbrev
+ nullptr, // direction
+ &toRadians,
+ nullptr // unit name
+ );
+ isLatFirst =
+ NS_PROJ::internal::ci_find(std::string(axisName), "latitude") !=
+ std::string::npos;
+
+ proj_obj_unref(cs);
+ }
+
+ return crs;
+}
+
+/************************************************************************/
+/* get_geog_crs_proj_string_from_proj_crs() */
+/************************************************************************/
+
+static std::string get_geog_crs_proj_string_from_proj_crs(PJ_OBJ *src,
+ double &toRadians,
+ bool &isLatFirst) {
+ auto srcType = proj_obj_get_type(src);
+ if (srcType == PJ_OBJ_TYPE_BOUND_CRS) {
+ auto base = proj_obj_get_source_crs(src);
+ assert(base);
+ proj_obj_unref(src);
+ src = base;
+ srcType = proj_obj_get_type(src);
+ }
+ if (srcType != PJ_OBJ_TYPE_PROJECTED_CRS) {
+ return std::string();
+ }
+
+ auto base = proj_obj_get_source_crs(src);
+ assert(base);
+ auto baseType = proj_obj_get_type(base);
+ if (baseType != PJ_OBJ_TYPE_GEOGRAPHIC_2D_CRS &&
+ baseType != PJ_OBJ_TYPE_GEOGRAPHIC_3D_CRS) {
+ proj_obj_unref(base);
+ return std::string();
+ }
+
+ auto cs = proj_obj_crs_get_coordinate_system(base);
+ assert(cs);
+
+ const char *axisName = "";
+ proj_obj_cs_get_axis_info(cs, 0,
+ &axisName, // name,
+ nullptr, // abbrev
+ nullptr, // direction
+ &toRadians,
+ nullptr // unit name
+ );
+ isLatFirst = NS_PROJ::internal::ci_find(std::string(axisName),
+ "latitude") != std::string::npos;
+
+ proj_obj_unref(cs);
+
+ auto retCStr = proj_obj_as_proj_string(base, PJ_PROJ_5, nullptr);
+ std::string ret(retCStr ? retCStr : "");
+ proj_obj_unref(base);
+ return ret;
+}
+
+/************************************************************************/
/* main() */
/************************************************************************/
int main(int argc, char **argv) {
char *arg;
char **eargv = argv;
- char *from_argv[MAX_PARGS];
- char *to_argv[MAX_PARGS];
+ std::string fromStr;
+ std::string toStr;
FILE *fid;
- int from_argc = 0, to_argc = 0, eargc = 0, mon = 0;
- int have_to_flag = 0, inverse = 0, i;
+ int eargc = 0, mon = 0;
+ int have_to_flag = 0, inverse = 0;
int use_env_locale = 0;
/* This is just to check that pj_init() is locale-safe */
@@ -185,6 +323,11 @@ int main(int argc, char **argv) {
if (getenv("PROJ_USE_ENV_LOCALE") != nullptr)
use_env_locale = 1;
+ /* Enable compatibility mode for init=epsg:XXXX by default */
+ if (getenv("PROJ_USE_PROJ4_INIT_RULES") == nullptr) {
+ proj_context_use_proj4_init_rules(nullptr, true);
+ }
+
if ((emess_dat.Prog_name = strrchr(*argv, DIR_CHAR)) != nullptr)
++emess_dat.Prog_name;
else
@@ -194,9 +337,27 @@ int main(int argc, char **argv) {
(void)fprintf(stderr, usage, pj_get_release(), emess_dat.Prog_name);
exit(0);
}
+
+ // First pass to check if we have "cs2cs [-bla]* <SRC> <DEST>" syntax
+ int countNonOptionArg = 0;
+ for (int i = 1; i < argc; i++) {
+ if (argv[i][0] == '-') {
+ if (argv[i][1] == '\0') {
+ countNonOptionArg++;
+ }
+ } else {
+ if (strcmp(argv[i], "+to") == 0) {
+ countNonOptionArg = -1;
+ break;
+ }
+ countNonOptionArg++;
+ }
+ }
+ const bool isSrcDestSyntax = (countNonOptionArg == 2);
+
/* process run line arguments */
while (--argc > 0) { /* collect run line arguments */
- if (**++argv == '-')
+ if (**++argv == '-') {
for (arg = *argv;;) {
switch (*++arg) {
case '\0': /* position of "stdin" */
@@ -325,21 +486,28 @@ int main(int argc, char **argv) {
}
break;
}
- else if (strcmp(*argv, "+to") == 0) {
+ } else if (isSrcDestSyntax) {
+ if (fromStr.empty())
+ fromStr = *argv;
+ else
+ toStr = *argv;
+ } else if (strcmp(*argv, "+to") == 0) {
have_to_flag = 1;
} else if (**argv == '+') { /* + argument */
if (have_to_flag) {
- if (to_argc < MAX_PARGS)
- to_argv[to_argc++] = *argv + 1;
- else
- emess(1, "overflowed + argument table");
+ if (!toStr.empty())
+ toStr += ' ';
+ toStr += *argv;
} else {
- if (from_argc < MAX_PARGS)
- from_argv[from_argc++] = *argv + 1;
- else
- emess(1, "overflowed + argument table");
+ if (!fromStr.empty())
+ fromStr += ' ';
+ fromStr += *argv;
}
+ } else if (!have_to_flag) {
+ fromStr = *argv;
+ } else if (toStr.empty()) {
+ toStr = *argv;
} else /* assumed to be input file name(s) */
eargv[eargc++] = *argv;
}
@@ -351,17 +519,7 @@ int main(int argc, char **argv) {
* coordinate systems.
*/
if (inverse) {
- int argcount;
-
- for (i = 0; i < MAX_PARGS; i++) {
- arg = from_argv[i];
- from_argv[i] = to_argv[i];
- to_argv[i] = arg;
- }
-
- argcount = from_argc;
- from_argc = to_argc;
- to_argc = argcount;
+ std::swap(fromStr, toStr);
}
if (use_env_locale) {
@@ -369,48 +527,64 @@ int main(int argc, char **argv) {
setlocale(LC_ALL, "");
}
- if (from_argc == 0 && to_argc != 0) {
- /* we will generate the from proj as the latlong of the +to in a bit */
- } else if (!(fromProj = pj_init(from_argc, from_argv))) {
- printf("Using from definition: ");
- for (i = 0; i < from_argc; i++)
- printf("%s ", from_argv[i]);
- printf("\n");
-
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
+ if (fromStr.empty() && toStr.empty()) {
+ emess(3, "missing source and target coordinate systems");
}
- if (to_argc == 0) {
- if (!(toProj = pj_latlong_from_proj(fromProj))) {
- printf("Using to definition: ");
- for (i = 0; i < to_argc; i++)
- printf("%s ", to_argv[i]);
- printf("\n");
+ const char *const optionsProj4Mode[] = {"USE_PROJ4_INIT_RULES=YES",
+ nullptr};
+ const char *const *optionsImportCRS =
+ proj_context_get_use_proj4_init_rules(nullptr) ? optionsProj4Mode
+ : nullptr;
+
+ PJ_OBJ *src = nullptr;
+ if (!fromStr.empty()) {
+ bool ignored;
+ src = instanciate_crs(fromStr, optionsImportCRS, srcIsGeog,
+ srcToRadians, ignored);
+ if (!src) {
+ emess(3, "cannot instanciate source coordinate system");
+ }
+ }
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
+ PJ_OBJ *dst = nullptr;
+ if (!toStr.empty()) {
+ dst = instanciate_crs(toStr, optionsImportCRS, destIsGeog,
+ destToRadians, destIsLatLong);
+ if (!dst) {
+ emess(3, "cannot instanciate target coordinate system");
}
- } else if (!(toProj = pj_init(to_argc, to_argv))) {
- printf("Using to definition: ");
- for (i = 0; i < to_argc; i++)
- printf("%s ", to_argv[i]);
- printf("\n");
+ }
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
+ if (toStr.empty()) {
+ assert(src);
+ toStr = get_geog_crs_proj_string_from_proj_crs(src, destToRadians,
+ destIsLatLong);
+ if (toStr.empty()) {
+ emess(3,
+ "missing target CRS and source CRS is not a projected CRS");
+ }
+ destIsGeog = true;
+ } else if (fromStr.empty()) {
+ assert(dst);
+ bool ignored;
+ fromStr =
+ get_geog_crs_proj_string_from_proj_crs(dst, srcToRadians, ignored);
+ if (fromStr.empty()) {
+ emess(3,
+ "missing source CRS and target CRS is not a projected CRS");
+ }
+ srcIsGeog = true;
}
- if (from_argc == 0 && toProj != nullptr) {
- if (!(fromProj = pj_latlong_from_proj(toProj))) {
- printf("Using to definition: ");
- for (i = 0; i < to_argc; i++)
- printf("%s ", to_argv[i]);
- printf("\n");
+ proj_obj_unref(src);
+ proj_obj_unref(dst);
- emess(3, "projection initialization failure\ncause: %s",
- pj_strerrno(pj_errno));
- }
+ transformation = proj_create_crs_to_crs(nullptr, fromStr.c_str(),
+ toStr.c_str(), nullptr);
+ if (!transformation) {
+ emess(3, "cannot initialize transformation\ncause: %s",
+ pj_strerrno(pj_errno));
}
if (use_env_locale) {
@@ -420,19 +594,19 @@ int main(int argc, char **argv) {
if (mon) {
printf("%c ---- From Coordinate System ----\n", tag);
- pj_pr_list(fromProj);
+ printf("%s\n", fromStr.c_str());
printf("%c ---- To Coordinate System ----\n", tag);
- pj_pr_list(toProj);
+ printf("%s\n", toStr.c_str());
}
/* set input formatting control */
- if (!fromProj->is_latlong)
+ if (!srcIsGeog)
informat = strtod;
else {
informat = dmstor;
}
- if (!toProj->is_latlong && !oform)
+ if (!destIsGeog && !oform)
oform = "%.2f";
/* process input file list */
@@ -454,8 +628,7 @@ int main(int argc, char **argv) {
emess_dat.File_name = nullptr;
}
- pj_free(fromProj);
- pj_free(toProj);
+ proj_destroy(transformation);
pj_deallocate_grids();
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-27 12:48:07 +0000