21 #include <grass/gis.h> 22 #include <grass/gprojects.h> 23 #include <grass/glocale.h> 26 #define MULTIPLY_LOOP(x,y,c,m) \ 29 for (i = 0; i < c; ++i) {\ 35 #define DIVIDE_LOOP(x,y,c,m) \ 38 for (i = 0; i < c; ++i) {\ 44 static double METERS_in = 1.0, METERS_out = 1.0;
85 const struct pj_info *info_out,
86 struct pj_info *info_trans)
88 if (info_in->pj ==
NULL)
92 if (!info_trans->def) {
93 if (info_out->pj !=
NULL && info_out->def !=
NULL)
94 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s +step %s",
95 info_in->def, info_out->def);
98 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
102 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
103 if (info_trans->pj ==
NULL) {
104 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
108 info_trans->meters = 1.;
109 info_trans->zone = 0;
110 sprintf(info_trans->proj,
"pipeline");
112 if (info_out->pj ==
NULL) {
114 G_warning(_(
"Unable to create latlong equivalent for '%s'"),
152 const struct pj_info *info_out,
153 const struct pj_info *info_trans,
int dir,
154 double *
x,
double *y,
double *z)
160 int in_is_ll, out_is_ll;
163 if (info_in->pj ==
NULL)
166 if (info_trans->pj ==
NULL)
171 METERS_in = info_in->meters;
172 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
174 METERS_out = info_out->meters;
175 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
184 METERS_out = info_in->meters;
185 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
187 METERS_in = info_out->meters;
188 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
199 c.lpzt.lam = (*x) / RAD_TO_DEG;
200 c.lpzt.phi = (*y) / RAD_TO_DEG;
208 c.xyzt.x = *x * METERS_in;
209 c.xyzt.y = *y * METERS_in;
217 c = proj_trans(info_trans->pj, dir, c);
218 ok = proj_errno(info_trans->pj);
221 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
228 *x = c.lpzt.lam * RAD_TO_DEG;
229 *y = c.lpzt.phi * RAD_TO_DEG;
235 *x = c.xyzt.x / METERS_out;
236 *y = c.xyzt.y / METERS_out;
244 const struct pj_info *p_in, *p_out;
246 if (info_out ==
NULL)
258 METERS_in = p_in->meters;
259 METERS_out = p_out->meters;
264 if (strncmp(p_in->proj,
"ll", 2) == 0) {
265 u = (*x) / RAD_TO_DEG;
266 v = (*y) / RAD_TO_DEG;
273 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
276 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
280 if (strncmp(p_out->proj,
"ll", 2) == 0) {
322 const struct pj_info *info_out,
323 const struct pj_info *info_trans,
int dir,
324 double *
x,
double *y,
double *z,
int n)
332 int in_is_ll, out_is_ll;
335 if (info_trans->pj ==
NULL)
340 METERS_in = info_in->meters;
341 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
343 METERS_out = info_out->meters;
344 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
353 METERS_out = info_in->meters;
354 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
356 METERS_in = info_out->meters;
357 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
366 z = G_malloc(
sizeof(
double) * n);
368 for (i = 0; i < n; i++)
383 for (i = 0; i < n; i++) {
385 c.lpzt.lam = x[i] / RAD_TO_DEG;
386 c.lpzt.phi = y[i] / RAD_TO_DEG;
388 c = proj_trans(info_trans->pj, dir, c);
389 if ((ok = proj_errno(info_trans->pj)) < 0)
392 x[i] = c.lp.lam * RAD_TO_DEG;
393 y[i] = c.lp.phi * RAD_TO_DEG;
397 for (i = 0; i < n; i++) {
399 c.lpzt.lam = x[i] / RAD_TO_DEG;
400 c.lpzt.phi = y[i] / RAD_TO_DEG;
402 c = proj_trans(info_trans->pj, dir, c);
403 if ((ok = proj_errno(info_trans->pj)) < 0)
406 x[i] = c.xy.x / METERS_out;
407 y[i] = c.xy.y / METERS_out;
414 for (i = 0; i < n; i++) {
416 c.xyzt.x = x[i] * METERS_in;
417 c.xyzt.y = y[i] * METERS_in;
419 c = proj_trans(info_trans->pj, dir, c);
420 if ((ok = proj_errno(info_trans->pj)) < 0)
423 x[i] = c.lp.lam * RAD_TO_DEG;
424 y[i] = c.lp.phi * RAD_TO_DEG;
428 for (i = 0; i < n; i++) {
430 c.xyzt.x = x[i] * METERS_in;
431 c.xyzt.y = y[i] * METERS_in;
433 c = proj_trans(info_trans->pj, dir, c);
434 if ((ok = proj_errno(info_trans->pj)) < 0)
437 x[i] = c.xy.x / METERS_out;
438 y[i] = c.xy.y / METERS_out;
446 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
450 const struct pj_info *p_in, *p_out;
461 METERS_in = p_in->meters;
462 METERS_out = p_out->meters;
465 z = G_malloc(
sizeof(
double) * n);
467 for (i = 0; i < n; ++i)
471 if (strncmp(p_in->proj,
"ll", 2) == 0) {
472 if (strncmp(p_out->proj,
"ll", 2) == 0) {
474 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
479 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
484 if (strncmp(p_out->proj,
"ll", 2) == 0) {
486 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
491 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
499 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
527 const struct pj_info *info_in,
const struct pj_info *info_out)
531 struct pj_info info_trans;
538 METERS_in = info_in->meters;
539 METERS_out = info_out->meters;
541 if (strncmp(info_in->proj,
"ll", 2) == 0) {
543 c.lpzt.lam = (*x) / RAD_TO_DEG;
544 c.lpzt.phi = (*y) / RAD_TO_DEG;
547 c = proj_trans(info_trans.pj, PJ_FWD, c);
548 ok = proj_errno(info_trans.pj);
550 if (strncmp(info_out->proj,
"ll", 2) == 0) {
552 *x = c.lp.lam * RAD_TO_DEG;
553 *y = c.lp.phi * RAD_TO_DEG;
557 *x = c.xy.x / METERS_out;
558 *y = c.xy.y / METERS_out;
563 c.xyzt.x = *x * METERS_in;
564 c.xyzt.y = *y * METERS_in;
567 c = proj_trans(info_trans.pj, PJ_FWD, c);
568 ok = proj_errno(info_trans.pj);
570 if (strncmp(info_out->proj,
"ll", 2) == 0) {
572 *x = c.lp.lam * RAD_TO_DEG;
573 *y = c.lp.phi * RAD_TO_DEG;
577 *x = c.xy.x / METERS_out;
578 *y = c.xy.y / METERS_out;
581 proj_destroy(info_trans.pj);
584 G_warning(_(
"proj_trans() failed: %d"), ok);
590 METERS_in = info_in->meters;
591 METERS_out = info_out->meters;
593 if (strncmp(info_in->proj,
"ll", 2) == 0) {
594 if (strncmp(info_out->proj,
"ll", 2) == 0) {
595 u = (*x) / RAD_TO_DEG;
596 v = (*y) / RAD_TO_DEG;
597 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
602 u = (*x) / RAD_TO_DEG;
603 v = (*y) / RAD_TO_DEG;
604 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
610 if (strncmp(info_out->proj,
"ll", 2) == 0) {
613 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
620 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
626 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
656 const struct pj_info *info_in,
const struct pj_info *info_out)
662 struct pj_info info_trans;
669 METERS_in = info_in->meters;
670 METERS_out = info_out->meters;
673 h = G_malloc(
sizeof *h * count);
675 for (i = 0; i <
count; ++i)
680 if (strncmp(info_in->proj,
"ll", 2) == 0) {
682 if (strncmp(info_out->proj,
"ll", 2) == 0) {
683 for (i = 0; i <
count; i++) {
685 c.lpzt.lam = x[i] / RAD_TO_DEG;
686 c.lpzt.phi = y[i] / RAD_TO_DEG;
688 c = proj_trans(info_trans.pj, PJ_FWD, c);
689 if ((ok = proj_errno(info_trans.pj)) < 0)
692 x[i] = c.lp.lam * RAD_TO_DEG;
693 y[i] = c.lp.phi * RAD_TO_DEG;
697 for (i = 0; i <
count; i++) {
699 c.lpzt.lam = x[i] / RAD_TO_DEG;
700 c.lpzt.phi = y[i] / RAD_TO_DEG;
702 c = proj_trans(info_trans.pj, PJ_FWD, c);
703 if ((ok = proj_errno(info_trans.pj)) < 0)
706 x[i] = c.xy.x / METERS_out;
707 y[i] = c.xy.y / METERS_out;
713 if (strncmp(info_out->proj,
"ll", 2) == 0) {
714 for (i = 0; i <
count; i++) {
716 c.xyzt.x = x[i] * METERS_in;
717 c.xyzt.y = y[i] * METERS_in;
719 c = proj_trans(info_trans.pj, PJ_FWD, c);
720 if ((ok = proj_errno(info_trans.pj)) < 0)
723 x[i] = c.lp.lam * RAD_TO_DEG;
724 y[i] = c.lp.phi * RAD_TO_DEG;
728 for (i = 0; i <
count; i++) {
730 c.xyzt.x = x[i] * METERS_in;
731 c.xyzt.y = y[i] * METERS_in;
733 c = proj_trans(info_trans.pj, PJ_FWD, c);
734 if ((ok = proj_errno(info_trans.pj)) < 0)
737 x[i] = c.xy.x / METERS_out;
738 y[i] = c.xy.y / METERS_out;
744 proj_destroy(info_trans.pj);
747 G_warning(_(
"proj_trans() failed: %d"), ok);
750 METERS_in = info_in->meters;
751 METERS_out = info_out->meters;
756 h = G_malloc(
sizeof *h * count);
758 for (i = 0; i <
count; ++i)
762 if (strncmp(info_in->proj,
"ll", 2) == 0) {
763 if (strncmp(info_out->proj,
"ll", 2) == 0) {
765 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
770 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
775 if (strncmp(info_out->proj,
"ll", 2) == 0) {
777 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
782 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
790 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
int pj_do_transform(int count, double *x, double *y, double *h, const struct pj_info *info_in, const struct pj_info *info_out)
Re-project an array of points between two co-ordinate systems with optional ellipsoidal height conver...
int pj_do_proj(double *x, double *y, const struct pj_info *info_in, const struct pj_info *info_out)
Re-project a point between two co-ordinate systems.
int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
Define a latitude / longitude co-ordinate system with the same ellipsoid and datum parameters as an e...
int GPJ_init_transform(const struct pj_info *info_in, const struct pj_info *info_out, struct pj_info *info_trans)
Create a PROJ transformation object to transform coordinates from an input SRS to an output SRS...
#define DIVIDE_LOOP(x, y, c, m)
int G_asprintf(char **out, const char *fmt,...)
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
int GPJ_transform_array(const struct pj_info *info_in, const struct pj_info *info_out, const struct pj_info *info_trans, int dir, double *x, double *y, double *z, int n)
Re-project an array of points between two co-ordinate systems using a transformation object prepared ...
#define MULTIPLY_LOOP(x, y, c, m)
int GPJ_transform(const struct pj_info *info_in, const struct pj_info *info_out, const struct pj_info *info_trans, int dir, double *x, double *y, double *z)
Re-project a point between two co-ordinate systems using a transformation object prepared with GPJ_pr...
void G_free(void *buf)
Free allocated memory.
void G_warning(const char *msg,...)
Print a warning message to stderr.