21 #include <grass/gis.h> 22 #include <grass/gprojects.h> 23 #include <grass/glocale.h> 26 #define MULTIPLY_LOOP(x,y,c,m) \ 28 for (i = 0; i < c; ++i) {\ 34 #define DIVIDE_LOOP(x,y,c,m) \ 36 for (i = 0; i < c; ++i) {\ 42 static double METERS_in = 1.0, METERS_out = 1.0;
83 const struct pj_info *info_out,
84 struct pj_info *info_trans)
86 if (info_in->pj ==
NULL)
90 info_trans->pj =
NULL;
91 if (!info_trans->def) {
92 if (info_in->srid && info_out->pj && info_out->srid) {
93 char *insrid, *outsrid;
95 #if PROJ_VERSION_MAJOR >= 6 97 if (strncmp(info_in->srid,
"epsg", 4) == 0)
100 insrid =
G_store(info_in->srid);
102 if (strncmp(info_out->srid,
"epsg", 4) == 0)
105 outsrid =
G_store(info_out->srid);
112 if (strncmp(info_in->srid,
"EPSG", 4) == 0)
115 insrid =
G_store(info_in->srid);
117 if (strncmp(info_out->srid,
"EPSG", 4) == 0)
120 outsrid =
G_store(info_out->srid);
124 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
129 if (info_trans->pj ==
NULL) {
130 G_warning(_(
"proj_create_crs_to_crs() failed for '%s' and '%s'"),
133 #if PROJ_VERSION_MAJOR >= 6 135 const char *str = proj_as_proj_string(
NULL, info_trans->pj,
139 info_trans->def =
G_store(str);
145 if (info_trans->pj ==
NULL) {
150 if (info_out->pj !=
NULL && info_out->def !=
NULL)
151 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s +step %s",
152 info_in->def, info_out->def);
155 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
160 if (info_trans->pj ==
NULL)
161 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
162 if (info_trans->pj ==
NULL) {
163 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
167 info_trans->meters = 1.;
168 info_trans->zone = 0;
169 sprintf(info_trans->proj,
"pipeline");
171 if (info_out->pj ==
NULL) {
173 G_warning(_(
"Unable to create latlong equivalent for '%s'"),
211 const struct pj_info *info_out,
212 const struct pj_info *info_trans,
int dir,
213 double *
x,
double *y,
double *z)
219 int in_is_ll, out_is_ll;
222 if (info_in->pj ==
NULL)
225 if (info_trans->pj ==
NULL)
230 METERS_in = info_in->meters;
231 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
233 METERS_out = info_out->meters;
234 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
243 METERS_out = info_in->meters;
244 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
246 METERS_in = info_out->meters;
247 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
261 c.lpzt.lam = (*x) / RAD_TO_DEG;
262 c.lpzt.phi = (*y) / RAD_TO_DEG;
270 c.xyzt.x = *x * METERS_in;
271 c.xyzt.y = *y * METERS_in;
279 c = proj_trans(info_trans->pj, dir, c);
280 ok = proj_errno(info_trans->pj);
283 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
293 *x = c.lpzt.lam * RAD_TO_DEG;
294 *y = c.lpzt.phi * RAD_TO_DEG;
300 *x = c.xyzt.x / METERS_out;
301 *y = c.xyzt.y / METERS_out;
309 const struct pj_info *p_in, *p_out;
311 if (info_out ==
NULL)
323 METERS_in = p_in->meters;
324 METERS_out = p_out->meters;
329 if (strncmp(p_in->proj,
"ll", 2) == 0) {
330 u = (*x) / RAD_TO_DEG;
331 v = (*y) / RAD_TO_DEG;
338 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
341 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
345 if (strncmp(p_out->proj,
"ll", 2) == 0) {
387 const struct pj_info *info_out,
388 const struct pj_info *info_trans,
int dir,
389 double *
x,
double *y,
double *z,
int n)
397 int in_is_ll, out_is_ll;
400 if (info_trans->pj ==
NULL)
405 METERS_in = info_in->meters;
406 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
408 METERS_out = info_out->meters;
409 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
418 METERS_out = info_in->meters;
419 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
421 METERS_in = info_out->meters;
422 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
431 z = G_malloc(
sizeof(
double) * n);
433 for (i = 0; i < n; i++)
448 for (i = 0; i < n; i++) {
450 c.lpzt.lam = x[i] / RAD_TO_DEG;
451 c.lpzt.phi = y[i] / RAD_TO_DEG;
453 c = proj_trans(info_trans->pj, dir, c);
454 if ((ok = proj_errno(info_trans->pj)) < 0)
457 x[i] = c.lp.lam * RAD_TO_DEG;
458 y[i] = c.lp.phi * RAD_TO_DEG;
462 for (i = 0; i < n; i++) {
464 c.lpzt.lam = x[i] / RAD_TO_DEG;
465 c.lpzt.phi = y[i] / RAD_TO_DEG;
467 c = proj_trans(info_trans->pj, dir, c);
468 if ((ok = proj_errno(info_trans->pj)) < 0)
471 x[i] = c.xy.x / METERS_out;
472 y[i] = c.xy.y / METERS_out;
479 for (i = 0; i < n; i++) {
481 c.xyzt.x = x[i] * METERS_in;
482 c.xyzt.y = y[i] * METERS_in;
484 c = proj_trans(info_trans->pj, dir, c);
485 if ((ok = proj_errno(info_trans->pj)) < 0)
488 x[i] = c.lp.lam * RAD_TO_DEG;
489 y[i] = c.lp.phi * RAD_TO_DEG;
493 for (i = 0; i < n; i++) {
495 c.xyzt.x = x[i] * METERS_in;
496 c.xyzt.y = y[i] * METERS_in;
498 c = proj_trans(info_trans->pj, dir, c);
499 if ((ok = proj_errno(info_trans->pj)) < 0)
502 x[i] = c.xy.x / METERS_out;
503 y[i] = c.xy.y / METERS_out;
511 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
515 const struct pj_info *p_in, *p_out;
526 METERS_in = p_in->meters;
527 METERS_out = p_out->meters;
530 z = G_malloc(
sizeof(
double) * n);
532 for (i = 0; i < n; ++i)
536 if (strncmp(p_in->proj,
"ll", 2) == 0) {
537 if (strncmp(p_out->proj,
"ll", 2) == 0) {
539 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
544 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
549 if (strncmp(p_out->proj,
"ll", 2) == 0) {
551 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
556 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
564 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
592 const struct pj_info *info_in,
const struct pj_info *info_out)
596 struct pj_info info_trans;
603 METERS_in = info_in->meters;
604 METERS_out = info_out->meters;
606 if (strncmp(info_in->proj,
"ll", 2) == 0) {
608 c.lpzt.lam = (*x) / RAD_TO_DEG;
609 c.lpzt.phi = (*y) / RAD_TO_DEG;
612 c = proj_trans(info_trans.pj, PJ_FWD, c);
613 ok = proj_errno(info_trans.pj);
615 if (strncmp(info_out->proj,
"ll", 2) == 0) {
617 *x = c.lp.lam * RAD_TO_DEG;
618 *y = c.lp.phi * RAD_TO_DEG;
622 *x = c.xy.x / METERS_out;
623 *y = c.xy.y / METERS_out;
628 c.xyzt.x = *x * METERS_in;
629 c.xyzt.y = *y * METERS_in;
632 c = proj_trans(info_trans.pj, PJ_FWD, c);
633 ok = proj_errno(info_trans.pj);
635 if (strncmp(info_out->proj,
"ll", 2) == 0) {
637 *x = c.lp.lam * RAD_TO_DEG;
638 *y = c.lp.phi * RAD_TO_DEG;
642 *x = c.xy.x / METERS_out;
643 *y = c.xy.y / METERS_out;
646 proj_destroy(info_trans.pj);
649 G_warning(_(
"proj_trans() failed: %d"), ok);
655 METERS_in = info_in->meters;
656 METERS_out = info_out->meters;
658 if (strncmp(info_in->proj,
"ll", 2) == 0) {
659 if (strncmp(info_out->proj,
"ll", 2) == 0) {
660 u = (*x) / RAD_TO_DEG;
661 v = (*y) / RAD_TO_DEG;
662 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
667 u = (*x) / RAD_TO_DEG;
668 v = (*y) / RAD_TO_DEG;
669 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
675 if (strncmp(info_out->proj,
"ll", 2) == 0) {
678 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
685 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
691 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
721 const struct pj_info *info_in,
const struct pj_info *info_out)
727 struct pj_info info_trans;
734 METERS_in = info_in->meters;
735 METERS_out = info_out->meters;
738 h = G_malloc(
sizeof *h * count);
740 for (i = 0; i <
count; ++i)
745 if (strncmp(info_in->proj,
"ll", 2) == 0) {
747 if (strncmp(info_out->proj,
"ll", 2) == 0) {
748 for (i = 0; i <
count; i++) {
750 c.lpzt.lam = x[i] / RAD_TO_DEG;
751 c.lpzt.phi = y[i] / RAD_TO_DEG;
753 c = proj_trans(info_trans.pj, PJ_FWD, c);
754 if ((ok = proj_errno(info_trans.pj)) < 0)
757 x[i] = c.lp.lam * RAD_TO_DEG;
758 y[i] = c.lp.phi * RAD_TO_DEG;
762 for (i = 0; i <
count; i++) {
764 c.lpzt.lam = x[i] / RAD_TO_DEG;
765 c.lpzt.phi = y[i] / RAD_TO_DEG;
767 c = proj_trans(info_trans.pj, PJ_FWD, c);
768 if ((ok = proj_errno(info_trans.pj)) < 0)
771 x[i] = c.xy.x / METERS_out;
772 y[i] = c.xy.y / METERS_out;
778 if (strncmp(info_out->proj,
"ll", 2) == 0) {
779 for (i = 0; i <
count; i++) {
781 c.xyzt.x = x[i] * METERS_in;
782 c.xyzt.y = y[i] * METERS_in;
784 c = proj_trans(info_trans.pj, PJ_FWD, c);
785 if ((ok = proj_errno(info_trans.pj)) < 0)
788 x[i] = c.lp.lam * RAD_TO_DEG;
789 y[i] = c.lp.phi * RAD_TO_DEG;
793 for (i = 0; i <
count; i++) {
795 c.xyzt.x = x[i] * METERS_in;
796 c.xyzt.y = y[i] * METERS_in;
798 c = proj_trans(info_trans.pj, PJ_FWD, c);
799 if ((ok = proj_errno(info_trans.pj)) < 0)
802 x[i] = c.xy.x / METERS_out;
803 y[i] = c.xy.y / METERS_out;
809 proj_destroy(info_trans.pj);
812 G_warning(_(
"proj_trans() failed: %d"), ok);
815 METERS_in = info_in->meters;
816 METERS_out = info_out->meters;
819 h = G_malloc(
sizeof *h * count);
821 for (i = 0; i <
count; ++i)
825 if (strncmp(info_in->proj,
"ll", 2) == 0) {
826 if (strncmp(info_out->proj,
"ll", 2) == 0) {
828 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
833 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
838 if (strncmp(info_out->proj,
"ll", 2) == 0) {
840 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
845 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
853 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...
char * G_store(const char *s)
Copy string to allocated memory.
#define DIVIDE_LOOP(x, y, c, m)
int G_asprintf(char **out, const char *fmt,...)
char * G_store_lower(const char *s)
Copy string to allocated memory and convert copied string to lower case.
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 ...
char * G_store_upper(const char *s)
Copy string to allocated memory and convert copied string to upper case.
#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.