From 608c7c6f568cab02abcebcbc162db47dfdb6da8a Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 12 Jul 2017 10:15:16 +0200 Subject: Implemented proj_create_crs_to_crs() --- src/PJ_cart.c | 23 +++++++++++++++++++++ src/pj_obs_api.c | 31 +++++++++++++++++++++++++++- src/proj.def | 62 ++++++++++++++++++++++++++++---------------------------- src/proj.h | 1 + 4 files changed, 85 insertions(+), 32 deletions(-) (limited to 'src') diff --git a/src/PJ_cart.c b/src/PJ_cart.c index 26acac2d..9155dd03 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -422,6 +422,29 @@ int pj_cart_selftest (void) { /* Clean up */ proj_destroy (P); + + /* test proj_create_crs_to_crs() */ + P = proj_create_crs_to_crs(0, "epsg:25832", "epsg:25833"); + if (P==0) + return 30; + + a.coo.xy.x = 700000.0; + a.coo.xy.y = 6000000.0; + b.coo.xy.x = 307788.8761171057; + b.coo.xy.y = 5999669.3036037628; + + a = proj_trans_obs(P, PJ_FWD, a); + if (dist > 1e-7) + return 31; + proj_destroy(P); + + /* let's make sure that only entries in init-files results in a usable PJ */ + P = proj_create_crs_to_crs(0, "proj=utm +zone=32 +datum=WGS84", "proj=utm +zone=33 +datum=WGS84"); + if (P != 0) { + proj_destroy(P); + return 32; + } + return 0; } #endif diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c index d8395b0f..e52d9c1c 100644 --- a/src/pj_obs_api.c +++ b/src/pj_obs_api.c @@ -314,6 +314,36 @@ PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv) { return pj_init_ctx (ctx, argc, argv); } +/*****************************************************************************/ +PJ *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *srid_from, const char *srid_to) { +/****************************************************************************** + Create a transformation pipeline between two known coordinate reference + systems. + + srid_from and srid_to should be the value part of a +init=... parameter + set, i.e. "epsg:25833" or "IGNF:AMST63". Any projection definition that is + can be found in a init-file in PROJ_LIB is a valid input to this function. + + For now the function mimics the cs2cs app: An input and an output CRS is + given and coordinates are transformed via a hub datum (WGS84). This + transformation strategy is referred to as "early-binding" by the EPSG. The + function can be extended to support "late-binding" transformations in the + future without affecting users of the function. + + Example call: + + PJ *P = proj_create_crs_to_crs(0, "epsg:25832", "epsg:25833"); + +******************************************************************************/ + PJ *P; + char buffer[256]; + + sprintf(buffer, "+proj=pipeline +step +init=%s +inv +step +init=%s", srid_from, srid_to); + P = proj_create(ctx, buffer); + + return P; +} + PJ *proj_destroy (PJ *P) { pj_free (P); return 0; @@ -447,4 +477,3 @@ double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians) int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs); int proj_transform_coord (PJ *P, enum proj_direction direction, size_t n, PJ_COORD *coord); PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags); -PJ *proj_create_crs_to_crs (PJ_CONTEXT *ctx, const char *def_from, const char *def_to); diff --git a/src/proj.def b/src/proj.def index 314c4385..c8f8be10 100644 --- a/src/proj.def +++ b/src/proj.def @@ -92,42 +92,42 @@ EXPORTS proj_create @90 proj_create_argv @91 - proj_destroy @92 + proj_create_crs_to_crs @92 + proj_destroy @93 - proj_trans_obs @93 - proj_trans_coord @94 - proj_transform @95 - proj_roundtrip @96 + proj_trans_obs @94 + proj_trans_coord @95 + proj_transform @96 + proj_roundtrip @97 - proj_coord @97 - proj_coord_error @98 - proj_obs_error @99 + proj_coord @98 + proj_coord_error @99 + proj_obs_error @100 - proj_errno @100 - proj_errno_set @101 - proj_errno_reset @102 - proj_errno_restore @103 - proj_context_errno_set @104 + proj_errno @101 + proj_errno_set @102 + proj_errno_reset @103 + proj_errno_restore @104 + proj_context_errno_set @105 - proj_context_create @105 - proj_context_set @106 - proj_context_inherit @107 - proj_context_destroy @108 + proj_context_create @106 + proj_context_set @107 + proj_context_inherit @108 + proj_context_destroy @109 - proj_lp_dist @109 - proj_xy_dist @110 - proj_xyz_dist @111 + proj_lp_dist @110 + proj_xy_dist @111 + proj_xyz_dist @112 - proj_log_level @112 - proj_log_func @113 - proj_log_error @114 - proj_log_debug @115 - proj_log_trace @116 + proj_log_level @113 + proj_log_func @114 + proj_log_error @115 + proj_log_debug @116 + proj_log_trace @117 - proj_definition_retrieve @117 - proj_release @118 - proj_torad @119 - proj_todeg @120 + proj_definition_retrieve @118 + proj_release @119 + proj_torad @120 + proj_todeg @121 - - pj_find_file @121 + pj_find_file @122 diff --git a/src/proj.h b/src/proj.h index 1977acde..b5a659d4 100644 --- a/src/proj.h +++ b/src/proj.h @@ -286,6 +286,7 @@ void proj_context_destroy (PJ_CONTEXT *ctx); /* Manage the transformation definition object PJ */ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition); PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv); +PJ *proj_create_crs_to_crs(PJ_CONTEXT *ctx, const char *srid_from, const char *srid_to); PJ *proj_destroy (PJ *P); -- cgit v1.2.3 From 8731336c9101b1948b401e0b0d8149598d56fa90 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 12 Jul 2017 10:42:44 +0200 Subject: Implemented proj_obs() --- src/pj_obs_api.c | 20 +++++++++++++++++++- src/proj.h | 4 ++-- 2 files changed, 21 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c index e52d9c1c..1b353b46 100644 --- a/src/pj_obs_api.c +++ b/src/pj_obs_api.c @@ -42,6 +42,7 @@ #include +/* Initialize PJ_COORD struct */ PJ_COORD proj_coord (double x, double y, double z, double t) { PJ_COORD res; res.v[0] = x; @@ -51,6 +52,24 @@ PJ_COORD proj_coord (double x, double y, double z, double t) { return res; } +/* Initialize PJ_OBS struct */ +PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags) { + PJ_OBS res; + res.coo.v[0] = x; + res.coo.v[1] = y; + res.coo.v[2] = z; + res.coo.v[3] = t; + res.anc.v[0] = o; + res.anc.v[1] = p; + res.anc.v[2] = k; + res.id = id; + res.flags = flags; + + return res; +} + + + /* Geodesic distance between two points with angular 2D coordinates */ double proj_lp_dist (PJ *P, LP a, LP b) { double s12, azi1, azi2; @@ -476,4 +495,3 @@ double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians) int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs); int proj_transform_coord (PJ *P, enum proj_direction direction, size_t n, PJ_COORD *coord); -PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags); diff --git a/src/proj.h b/src/proj.h index b5a659d4..d15a798f 100644 --- a/src/proj.h +++ b/src/proj.h @@ -312,9 +312,9 @@ size_t proj_transform ( ); -/* not a constructor, but an initializer */ +/* Initializers */ PJ_COORD proj_coord (double x, double y, double z, double t); - +PJ_OBS proj_obs (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags); /* Measure internal consistency - in forward or inverse direction */ double proj_roundtrip (PJ *P, enum proj_direction direction, int n, PJ_OBS obs); -- cgit v1.2.3 From c3fd68227b408a0172afb449651d3dae6b38612e Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 12 Jul 2017 11:52:14 +0200 Subject: Implemented proj_transform_coord() and proj_transform_obs() --- src/PJ_cart.c | 39 ++++++++++++++++++++++++++++++++++----- src/pj_obs_api.c | 52 +++++++++++++++++++++++++++++++++++++++++----------- src/proj.def | 55 +++++++++++++++++++++++++++++-------------------------- src/proj.h | 2 ++ 4 files changed, 106 insertions(+), 42 deletions(-) (limited to 'src') diff --git a/src/PJ_cart.c b/src/PJ_cart.c index 9155dd03..c6d1dd74 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -235,6 +235,7 @@ int pj_cart_selftest (void) { PJ_CONTEXT *ctx; PJ *P; PJ_OBS a, b, obs[2]; + PJ_COORD coord[2]; int err; size_t n, sz; double dist, h, t; @@ -419,14 +420,42 @@ int pj_cart_selftest (void) { if (50 != obs[1].coo.lpz.z) return 27; /* NOTE: unchanged */ if (50==h) return 28; - /* Clean up */ - proj_destroy (P); + /* test proj_transform_obs() */ + + obs[0].coo = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0); + obs[1].coo = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0); + + if (proj_transform_obs(P, PJ_FWD, 2, obs)) + return 30; + if (a.coo.lpz.lam != obs[0].coo.lpz.lam) return 31; + if (a.coo.lpz.phi != obs[0].coo.lpz.phi) return 32; + if (a.coo.lpz.z != obs[0].coo.lpz.z) return 33; + if (b.coo.lpz.lam != obs[1].coo.lpz.lam) return 34; + if (b.coo.lpz.phi != obs[1].coo.lpz.phi) return 35; + if (b.coo.lpz.z != obs[1].coo.lpz.z) return 36; + + /* test proj_transform_coord() */ + + coord[0] = proj_coord (PJ_TORAD(12), PJ_TORAD(55), 45, 0); + coord[1] = proj_coord (PJ_TORAD(12), PJ_TORAD(56), 50, 0); + if (proj_transform_coord(P, PJ_FWD, 2, coord)) + return 40; + + if (a.coo.lpz.lam != coord[0].lpz.lam) return 41; + if (a.coo.lpz.phi != coord[0].lpz.phi) return 42; + if (a.coo.lpz.z != coord[0].lpz.z) return 43; + if (b.coo.lpz.lam != coord[1].lpz.lam) return 44; + if (b.coo.lpz.phi != coord[1].lpz.phi) return 45; + if (b.coo.lpz.z != coord[1].lpz.z) return 46; + + /* Clean up after transform* tests */ + proj_destroy (P); /* test proj_create_crs_to_crs() */ P = proj_create_crs_to_crs(0, "epsg:25832", "epsg:25833"); if (P==0) - return 30; + return 50; a.coo.xy.x = 700000.0; a.coo.xy.y = 6000000.0; @@ -435,14 +464,14 @@ int pj_cart_selftest (void) { a = proj_trans_obs(P, PJ_FWD, a); if (dist > 1e-7) - return 31; + return 51; proj_destroy(P); /* let's make sure that only entries in init-files results in a usable PJ */ P = proj_create_crs_to_crs(0, "proj=utm +zone=32 +datum=WGS84", "proj=utm +zone=33 +datum=WGS84"); if (P != 0) { proj_destroy(P); - return 32; + return 52; } return 0; diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c index 1b353b46..81a6e337 100644 --- a/src/pj_obs_api.c +++ b/src/pj_obs_api.c @@ -216,7 +216,7 @@ size_t proj_transform ( In most cases, the stride will be identical for x, y,z, and t, since they will typically be either individual arrays (stride = sizeof(double)), or strided views into an array of application specific data structures (stride = sizeof (...)). - + But in order to support cases where x, y, z, and t come from heterogeneous sources, individual strides, sx, sy, sz, st, are used. @@ -244,7 +244,7 @@ size_t proj_transform ( if (0==ny) y = &null_broadcast; if (0==nz) z = &null_broadcast; if (0==nt) t = &null_broadcast; - + /* nothing to do? */ if (0==nx+ny+nz+nt) return 0; @@ -315,10 +315,46 @@ size_t proj_transform ( *z = coord.xyzt.z; if (nt==1) *t = coord.xyzt.t; - + return i; } +/*****************************************************************************/ +int proj_transform_obs (PJ *P, enum proj_direction direction, size_t n, PJ_OBS *obs) { +/****************************************************************************** + Batch transform an array of PJ_OBS. + + Returns 0 if all observations are transformed without error, otherwise + returns error number. +******************************************************************************/ + size_t i; + for (i=0; i Date: Wed, 12 Jul 2017 14:07:31 +0200 Subject: Added proj_has_inverse(). Fixes #155. --- src/PJ_cart.c | 17 +++++++++++++++++ src/pj_obs_api.c | 4 ++++ src/proj.def | 3 ++- src/proj.h | 5 +++++ 4 files changed, 28 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_cart.c b/src/PJ_cart.c index c6d1dd74..1109b01b 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -473,6 +473,23 @@ int pj_cart_selftest (void) { proj_destroy(P); return 52; } + proj_destroy(P); + + /* Test proj_has_inverse() */ + P = proj_create(0, "+proj=august"); /* august has no inverse */ + if (proj_has_inverse(P)) { + proj_destroy(P); + return 60; + } + proj_destroy(P); + + P = proj_create(0, "+proj=merc"); /* merc has an inverse */ + if (!proj_has_inverse(P)) { + proj_destroy(P); + return 61; + } + proj_destroy(P); + return 0; } diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c index 81a6e337..4fac1dac 100644 --- a/src/pj_obs_api.c +++ b/src/pj_obs_api.c @@ -525,3 +525,7 @@ void *proj_release (void *buffer) { double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);} double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);} + +int proj_has_inverse(PJ *P) { + return (P->inv != 0 || P->inv3d != 0 || P->invobs != 0); +} diff --git a/src/proj.def b/src/proj.def index a61151ed..abb17765 100644 --- a/src/proj.def +++ b/src/proj.def @@ -132,5 +132,6 @@ EXPORTS proj_release @122 proj_torad @123 proj_todeg @124 + proj_has_inverse @125 - pj_find_file @125 + pj_find_file @126 diff --git a/src/proj.h b/src/proj.h index f94fb384..c0b8b7cc 100644 --- a/src/proj.h +++ b/src/proj.h @@ -337,6 +337,7 @@ void proj_errno_set (PJ *P, int err); int proj_errno_reset (PJ *P); void proj_errno_restore (PJ *P, int err); + /* Build a fully expanded proj_create() compatible representation of P */ char *proj_definition_retrieve (PJ *P); /* ...and get rid of it safely */ @@ -349,6 +350,10 @@ void *proj_release (void *buffer); double proj_torad (double angle_in_degrees); double proj_todeg (double angle_in_radians); +/* Check if a projection has an inverse mapping */ +int proj_has_inverse(PJ *P); + + #ifdef __cplusplus } #endif -- cgit v1.2.3 From 40bc7c0bb8b422774be37fa7af0f80877bb52696 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Thu, 13 Jul 2017 14:00:47 +0200 Subject: Add proj_dmstor() and proj_rtodms() to proj.h API. Resolves #172. --- src/PJ_cart.c | 15 +++++++++++++++ src/pj_obs_api.c | 8 ++++++++ src/proj.def | 4 +++- src/proj.h | 3 +++ 4 files changed, 29 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PJ_cart.c b/src/PJ_cart.c index 1109b01b..a40ebdaf 100644 --- a/src/PJ_cart.c +++ b/src/PJ_cart.c @@ -242,6 +242,7 @@ int pj_cart_selftest (void) { char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"}; char *arg = {" +proj=utm +zone=32 +ellps=GRS80"}; char *arg_def; + char buf[40]; /* An utm projection on the GRS80 ellipsoid */ P = proj_create (0, arg); @@ -490,6 +491,20 @@ int pj_cart_selftest (void) { } proj_destroy(P); + /* test proj_rtodms() and proj_dmstor() */ + if (strcmp("180dN", proj_rtodms(buf, M_PI, 'N', 'S'))) + return 70; + + if (proj_dmstor(&buf[0], NULL) != M_PI) + return 71; + + if (strcmp("114d35'29.612\"S", proj_rtodms(buf, -2.0, 'N', 'S'))) + return 72; + + /* we can't expect perfect numerical accuracy so testing with a tolerance */ + if (fabs(-2.0 - proj_dmstor(&buf[0], NULL)) > 1e-7) + return 73; + return 0; } diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c index 4fac1dac..05fe58f0 100644 --- a/src/pj_obs_api.c +++ b/src/pj_obs_api.c @@ -529,3 +529,11 @@ double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians) int proj_has_inverse(PJ *P) { return (P->inv != 0 || P->inv3d != 0 || P->invobs != 0); } + +double proj_dmstor(const char *is, char **rs) { + return dmstor(is, rs); +} + +char* proj_rtodms(char *s, double r, int pos, int neg) { + return rtodms(s, r, pos, neg); +} diff --git a/src/proj.def b/src/proj.def index abb17765..5a762291 100644 --- a/src/proj.def +++ b/src/proj.def @@ -133,5 +133,7 @@ EXPORTS proj_torad @123 proj_todeg @124 proj_has_inverse @125 + proj_rtodms @126 + proj_dmstor @127 - pj_find_file @126 + pj_find_file @128 diff --git a/src/proj.h b/src/proj.h index c0b8b7cc..398368c4 100644 --- a/src/proj.h +++ b/src/proj.h @@ -353,6 +353,9 @@ double proj_todeg (double angle_in_radians); /* Check if a projection has an inverse mapping */ int proj_has_inverse(PJ *P); +double proj_dmstor(const char *is, char **rs); +char* proj_rtodms(char *s, double r, int pos, int neg); + #ifdef __cplusplus } -- cgit v1.2.3