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) {\
42static double METERS_in = 1.0, METERS_out = 1.0;
45#if PROJ_VERSION_MAJOR >= 6
46int get_pj_area(
const struct pj_info *iproj,
double *xmin,
double *xmax,
47 double *ymin,
double *ymax)
49 struct Cell_head window;
59 if (window.proj != PROJECTION_LL) {
64 const char *projstr =
NULL;
66 struct pj_info oproj, tproj;
72 if (proj_get_type(iproj->pj) == PJ_TYPE_BOUND_CRS) {
75 source_crs = proj_get_source_crs(
NULL, iproj->pj);
77 projstr = proj_as_proj_string(
NULL, source_crs, PJ_PROJ_5,
NULL);
81 proj_destroy(source_crs);
85 projstr = proj_as_proj_string(
NULL, iproj->pj, PJ_PROJ_5,
NULL);
94 G_asprintf(&tproj.def,
"+proj=pipeline +step +inv %s",
96 G_debug(1,
"get_pj_area() tproj.def: %s", tproj.def);
97 tproj.pj = proj_create(PJ_DEFAULT_CTX, tproj.def);
99 if (tproj.pj ==
NULL) {
100 G_warning(_(
"proj_create() failed for '%s'"), tproj.def);
103 proj_destroy(tproj.pj);
107 projstr = proj_as_proj_string(
NULL, tproj.pj, PJ_PROJ_5,
NULL);
108 if (projstr ==
NULL) {
109 G_warning(_(
"proj_create() failed for '%s'"), tproj.def);
112 proj_destroy(tproj.pj);
117 G_debug(1,
"proj_create() projstr '%s'", projstr);
121 estep = (window.west + window.east) / 21.;
122 nstep = (window.north + window.south) / 21.;
123 for (i = 0; i < 20; i++) {
124 x[i] = window.west + estep * (i + 1);
127 x[i + 20] = window.west + estep * (i + 1);
128 y[i + 20] = window.south;
130 x[i + 40] = window.west;
131 y[i + 40] = window.south + nstep * (i + 1);
133 x[i + 60] = window.east;
134 y[i + 60] = window.south + nstep * (i + 1);
137 y[80] = window.north;
139 y[81] = window.south;
141 y[82] = window.north;
143 y[83] = window.south;
144 x[84] = (window.west + window.east) / 2.;
145 y[84] = (window.north + window.south) / 2.;
149 proj_destroy(tproj.pj);
152 *xmin = *xmax =
x[84];
153 *ymin = *ymax = y[84];
154 for (i = 0; i < 84; i++) {
165 G_debug(1,
"input window north: %.8f", window.north);
166 G_debug(1,
"input window south: %.8f", window.south);
167 G_debug(1,
"input window east: %.8f", window.east);
168 G_debug(1,
"input window west: %.8f", window.west);
170 G_debug(1,
"transformed xmin: %.8f", *xmin);
171 G_debug(1,
"transformed xmax: %.8f", *xmax);
172 G_debug(1,
"transformed ymin: %.8f", *ymin);
173 G_debug(1,
"transformed ymax: %.8f", *ymax);
175 G_debug(1,
"get_pj_area(): xmin %g, xmax %g, ymin %g, ymax %g",
176 *xmin, *xmax, *ymin, *ymax);
181char *get_pj_type_string(PJ *pj)
183 char *pj_type =
NULL;
185 switch (proj_get_type(pj)) {
186 case PJ_TYPE_UNKNOWN:
189 case PJ_TYPE_ELLIPSOID:
192 case PJ_TYPE_PRIME_MERIDIAN:
195 case PJ_TYPE_GEODETIC_REFERENCE_FRAME:
196 G_asprintf(&pj_type,
"geodetic reference frame");
198 case PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME:
199 G_asprintf(&pj_type,
"dynamic geodetic reference frame");
201 case PJ_TYPE_VERTICAL_REFERENCE_FRAME:
202 G_asprintf(&pj_type,
"vertical reference frame");
204 case PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME:
205 G_asprintf(&pj_type,
"dynamic vertical reference frame");
207 case PJ_TYPE_DATUM_ENSEMBLE:
214 case PJ_TYPE_GEODETIC_CRS:
217 case PJ_TYPE_GEOCENTRIC_CRS:
222 case PJ_TYPE_GEOGRAPHIC_CRS:
225 case PJ_TYPE_GEOGRAPHIC_2D_CRS:
228 case PJ_TYPE_GEOGRAPHIC_3D_CRS:
231 case PJ_TYPE_VERTICAL_CRS:
234 case PJ_TYPE_PROJECTED_CRS:
237 case PJ_TYPE_COMPOUND_CRS:
240 case PJ_TYPE_TEMPORAL_CRS:
243 case PJ_TYPE_ENGINEERING_CRS:
246 case PJ_TYPE_BOUND_CRS:
249 case PJ_TYPE_OTHER_CRS:
252 case PJ_TYPE_CONVERSION:
255 case PJ_TYPE_TRANSFORMATION:
258 case PJ_TYPE_CONCATENATED_OPERATION:
259 G_asprintf(&pj_type,
"concatenated operation");
261 case PJ_TYPE_OTHER_COORDINATE_OPERATION:
262 G_asprintf(&pj_type,
"other coordinate operation");
273PJ *get_pj_object(
const struct pj_info *in_gpj,
char **in_defstr)
281 G_debug(1,
"Trying SRID '%s' ...", in_gpj->srid);
282 in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->srid);
284 G_warning(_(
"Unrecognized SRID '%s'"), in_gpj->srid);
287 *in_defstr =
G_store(in_gpj->srid);
291 ((
struct pj_info *)in_gpj)->meters = 1;
294 if (!in_pj && in_gpj->wkt) {
295 G_debug(1,
"Trying WKT '%s' ...", in_gpj->wkt);
296 in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->wkt);
298 G_warning(_(
"Unrecognized WKT '%s'"), in_gpj->wkt);
301 *in_defstr =
G_store(in_gpj->wkt);
305 ((
struct pj_info *)in_gpj)->meters = 1;
308 if (!in_pj && in_gpj->pj) {
309 in_pj = proj_clone(PJ_DEFAULT_CTX, in_gpj->pj);
311 if (*in_defstr && !**in_defstr)
316 G_warning(_(
"Unable to create PROJ object"));
332 if (proj_get_type(in_pj) == PJ_TYPE_BOUND_CRS) {
336 source_crs = proj_get_source_crs(
NULL, in_pj);
338 *in_defstr =
G_store(proj_as_wkt(
NULL, source_crs, PJ_WKT2_LATEST,
NULL));
339 if (*in_defstr && !**in_defstr)
389 const struct pj_info *info_out,
390 struct pj_info *info_trans)
392 if (info_in->pj ==
NULL)
395 if (info_in->def ==
NULL)
396 G_fatal_error(_(
"Input coordinate system definition is NULL"));
399#if PROJ_VERSION_MAJOR >= 6
404 info_trans->pj =
NULL;
405 info_trans->meters = 1.;
406 info_trans->zone = 0;
407 sprintf(info_trans->proj,
"pipeline");
410 if (info_trans->def) {
415 if (!info_in->pj || !info_in->proj[0] ||
416 !info_out->pj ||info_out->proj[0]) {
417 G_warning(_(
"A custom pipeline requires input and output projection info"));
424 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
425 if (info_trans->pj ==
NULL) {
426 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
430 projstr = proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
431 if (projstr ==
NULL) {
432 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
442 info_trans->def =
G_store(projstr);
444 if (strstr(info_trans->def,
"axisswap")) {
445 G_warning(_(
"The transformation pipeline contains an '%s' step. "
446 "Remove this step if easting and northing are swapped in the output."),
450 G_debug(1,
"proj_create() pipeline: %s", info_trans->def);
455 ((
struct pj_info *)info_in)->meters = 1;
456 ((
struct pj_info *)info_out)->meters = 1;
461 else if (info_out->pj ==
NULL) {
462 const char *projstr =
NULL;
467 G_debug(1,
"ll equivalent definition: %s", indef);
472 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
474 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
475 if (info_trans->pj ==
NULL) {
476 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
481 projstr = proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
482 if (projstr ==
NULL) {
483 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
491 else if (info_in->def && info_out->pj && info_out->def) {
493 char *insrid =
NULL, *outsrid =
NULL;
495 PJ_OBJ_LIST *op_list;
496 PJ_OPERATION_FACTORY_CONTEXT *operation_ctx;
497 PJ_AREA *pj_area =
NULL;
498 double xmin, xmax, ymin, ymax;
499 int op_count = 0, op_count_area = 0;
505 if (get_pj_area(info_in, &xmin, &xmax, &ymin, &ymax)) {
506 pj_area = proj_area_create();
507 proj_area_set_bbox(pj_area, xmin, ymin, xmax, ymax);
510 G_debug(1,
"source proj string: %s", info_in->def);
511 G_debug(1,
"source type: %s", get_pj_type_string(info_in->pj));
515 if (strncmp(info_in->srid,
"epsg", 4) == 0) {
518 ((
struct pj_info *)info_in)->srid = insrid;
523 in_pj = get_pj_object(info_in, &indef);
524 if (in_pj ==
NULL || indef ==
NULL) {
525 G_warning(_(
"Input CRS not available for '%s'"), info_in->def);
529 G_debug(1,
"Input CRS definition: %s", indef);
531 G_debug(1,
"target proj string: %s", info_out->def);
532 G_debug(1,
"target type: %s", get_pj_type_string(info_out->pj));
535 if (info_out->srid) {
536 if (strncmp(info_out->srid,
"epsg", 4) == 0) {
539 ((
struct pj_info *)info_out)->srid = outsrid;
544 out_pj = get_pj_object(info_out, &outdef);
545 if (out_pj ==
NULL || outdef ==
NULL) {
546 G_warning(_(
"Output CRS not available for '%s'"), info_out->def);
550 G_debug(1,
"Output CRS definition: %s", outdef);
554 operation_ctx = proj_create_operation_factory_context(PJ_DEFAULT_CTX,
NULL);
561 op_list = proj_create_operations(PJ_DEFAULT_CTX,
565 proj_operation_factory_context_destroy(operation_ctx);
569 op_count = proj_list_get_count(op_list);
574 for (i = 0; i < op_count; i++) {
575 const char *area_of_use, *projstr;
577 PJ_PROJ_INFO pj_info;
580 op = proj_list_get(PJ_DEFAULT_CTX, op_list, i);
581 op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
584 G_warning(_(
"proj_normalize_for_visualization() failed for operation %d"),
592 pj_info = proj_pj_info(op);
593 proj_get_area_of_use(
NULL, op, &w, &s, &e, &n, &area_of_use);
596 if (pj_info.description) {
598 pj_info.description);
605 if (pj_info.accuracy > 0) {
610#if PROJ_VERSION_NUM >= 6020000
611 const char *str = proj_get_remarks(op);
616 str = proj_get_scope(op);
623 projstr = proj_as_proj_string(
NULL, op,
643 proj_list_destroy(op_list);
651 operation_ctx = proj_create_operation_factory_context(PJ_DEFAULT_CTX,
NULL);
652 proj_operation_factory_context_set_area_of_interest(PJ_DEFAULT_CTX,
656 proj_operation_factory_context_set_spatial_criterion(PJ_DEFAULT_CTX,
658 PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT);
659 proj_operation_factory_context_set_grid_availability_use(PJ_DEFAULT_CTX,
661 PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
670 op_list = proj_create_operations(PJ_DEFAULT_CTX,
674 proj_operation_factory_context_destroy(operation_ctx);
677 op_count_area = proj_list_get_count(op_list);
678 if (op_count_area == 0) {
680 info_trans->pj =
NULL;
682 else if (op_count_area == 1) {
683 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
689 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
692 proj_list_destroy(op_list);
695 G_debug(1,
"trying %s to %s", indef, outdef);
716 proj_destroy(out_pj);
718 if (info_trans->pj) {
722 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ%d",
725 projstr = proj_as_proj_string(
NULL, info_trans->pj,
728 info_trans->def =
G_store(projstr);
736 pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, info_trans->pj);
739 G_warning(_(
"proj_normalize_for_visualization() failed for '%s'"),
743 projstr = proj_as_proj_string(
NULL, pj_norm,
745 if (projstr && *projstr) {
746 proj_destroy(info_trans->pj);
747 info_trans->pj = pj_norm;
748 info_trans->def =
G_store(projstr);
751 proj_destroy(pj_norm);
752 G_warning(_(
"No PROJ definition for normalized version of '%s'"),
761 proj_destroy(info_trans->pj);
762 info_trans->pj =
NULL;
767 proj_area_destroy(pj_area);
776 if (info_trans->pj ==
NULL) {
777 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
782 info_trans->pj =
NULL;
783 info_trans->meters = 1.;
784 info_trans->zone = 0;
785 sprintf(info_trans->proj,
"pipeline");
788 if (info_trans->def) {
790 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
791 if (info_trans->pj ==
NULL) {
792 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
799 else if (info_out->pj ==
NULL) {
800 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
802 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
803 if (info_trans->pj ==
NULL) {
804 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
809 else if (info_in->def && info_out->pj && info_out->def) {
811 char *insrid =
NULL, *outsrid =
NULL;
815 if (strncmp(info_in->srid,
"EPSG", 4) == 0)
818 insrid =
G_store(info_in->srid);
821 if (info_out->pj && info_out->srid) {
822 if (strncmp(info_out->srid,
"EPSG", 4) == 0)
825 outsrid =
G_store(info_out->srid);
828 info_trans->pj =
NULL;
830 if (insrid && outsrid) {
834 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv +init=%s +step +init=%s",
838 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
844 if (info_trans->pj) {
845 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ5");
871 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s +step %s",
874 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
883 if (info_trans->pj ==
NULL) {
884 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
891 if (info_out->pj ==
NULL) {
893 G_warning(_(
"Unable to create latlong equivalent for '%s'"),
932 const struct pj_info *info_out,
933 const struct pj_info *info_trans,
int dir,
934 double *
x,
double *y,
double *z)
940 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
943 if (info_in->pj ==
NULL)
946 if (info_trans->pj ==
NULL)
949 in_deg2rad = out_rad2deg = 1;
952 METERS_in = info_in->meters;
953 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
954#if PROJ_VERSION_MAJOR >= 6
958 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
963 METERS_out = info_out->meters;
964 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
965#if PROJ_VERSION_MAJOR >= 6
969 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
981 METERS_out = info_in->meters;
982 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
983#if PROJ_VERSION_MAJOR >= 6
987 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
992 METERS_in = info_out->meters;
993 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
994#if PROJ_VERSION_MAJOR >= 6
998 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1013 c.lpzt.lam = (*x) / RAD_TO_DEG;
1014 c.lpzt.phi = (*y) / RAD_TO_DEG;
1027 c.xyzt.x = *
x * METERS_in;
1028 c.xyzt.y = *y * METERS_in;
1035 G_debug(1,
"c.xyzt.x: %g", c.xyzt.x);
1036 G_debug(1,
"c.xyzt.y: %g", c.xyzt.y);
1037 G_debug(1,
"c.xyzt.z: %g", c.xyzt.z);
1040 c = proj_trans(info_trans->pj, dir, c);
1041 ok = proj_errno(info_trans->pj);
1044 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1053 *
x = c.lp.lam * RAD_TO_DEG;
1054 *y = c.lp.phi * RAD_TO_DEG;
1065 *
x = c.xyzt.x / METERS_out;
1066 *y = c.xyzt.y / METERS_out;
1074 const struct pj_info *p_in, *p_out;
1076 if (info_out ==
NULL)
1079 if (dir == PJ_FWD) {
1088 METERS_in = p_in->meters;
1089 METERS_out = p_out->meters;
1094 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1095 u = (*x) / RAD_TO_DEG;
1096 v = (*y) / RAD_TO_DEG;
1103 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1106 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1110 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1111 *
x = u * RAD_TO_DEG;
1112 *y = v * RAD_TO_DEG;
1115 *
x = u / METERS_out;
1116 *y = v / METERS_out;
1152 const struct pj_info *info_out,
1153 const struct pj_info *info_trans,
int dir,
1154 double *
x,
double *y,
double *z,
int n)
1162 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1165 if (info_trans->pj ==
NULL)
1168 in_deg2rad = out_rad2deg = 1;
1169 if (dir == PJ_FWD) {
1171 METERS_in = info_in->meters;
1172 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
1173#if PROJ_VERSION_MAJOR >= 6
1177 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1182 METERS_out = info_out->meters;
1183 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
1184#if PROJ_VERSION_MAJOR >= 6
1188 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1200 METERS_out = info_in->meters;
1201 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
1202#if PROJ_VERSION_MAJOR >= 6
1206 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1211 METERS_in = info_out->meters;
1212 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
1213#if PROJ_VERSION_MAJOR >= 6
1217 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1229 z = G_malloc(
sizeof(
double) * n);
1231 for (i = 0; i < n; i++)
1246 for (i = 0; i < n; i++) {
1249 c.lpzt.lam = (*x) / RAD_TO_DEG;
1250 c.lpzt.phi = (*y) / RAD_TO_DEG;
1257 c = proj_trans(info_trans->pj, dir, c);
1258 if ((ok = proj_errno(info_trans->pj)) < 0)
1262 x[i] = c.lp.lam * RAD_TO_DEG;
1263 y[i] = c.lp.phi * RAD_TO_DEG;
1272 for (i = 0; i < n; i++) {
1275 c.lpzt.lam = (*x) / RAD_TO_DEG;
1276 c.lpzt.phi = (*y) / RAD_TO_DEG;
1283 c = proj_trans(info_trans->pj, dir, c);
1284 if ((ok = proj_errno(info_trans->pj)) < 0)
1288 x[i] = c.xy.x / METERS_out;
1289 y[i] = c.xy.y / METERS_out;
1296 for (i = 0; i < n; i++) {
1298 c.xyzt.x =
x[i] * METERS_in;
1299 c.xyzt.y = y[i] * METERS_in;
1301 c = proj_trans(info_trans->pj, dir, c);
1302 if ((ok = proj_errno(info_trans->pj)) < 0)
1306 x[i] = c.lp.lam * RAD_TO_DEG;
1307 y[i] = c.lp.phi * RAD_TO_DEG;
1316 for (i = 0; i < n; i++) {
1318 c.xyzt.x =
x[i] * METERS_in;
1319 c.xyzt.y = y[i] * METERS_in;
1321 c = proj_trans(info_trans->pj, dir, c);
1322 if ((ok = proj_errno(info_trans->pj)) < 0)
1325 x[i] = c.xy.x / METERS_out;
1326 y[i] = c.xy.y / METERS_out;
1334 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1338 const struct pj_info *p_in, *p_out;
1340 if (dir == PJ_FWD) {
1349 METERS_in = p_in->meters;
1350 METERS_out = p_out->meters;
1353 z = G_malloc(
sizeof(
double) * n);
1355 for (i = 0; i < n; ++i)
1359 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1360 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1362 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1367 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1372 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1374 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1379 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1387 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1415 const struct pj_info *info_in,
const struct pj_info *info_out)
1419 struct pj_info info_trans;
1426 METERS_in = info_in->meters;
1427 METERS_out = info_out->meters;
1429 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1431 c.lpzt.lam = (*x) / RAD_TO_DEG;
1432 c.lpzt.phi = (*y) / RAD_TO_DEG;
1435 c = proj_trans(info_trans.pj, PJ_FWD, c);
1436 ok = proj_errno(info_trans.pj);
1438 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1440 *
x = c.lp.lam * RAD_TO_DEG;
1441 *y = c.lp.phi * RAD_TO_DEG;
1445 *
x = c.xy.x / METERS_out;
1446 *y = c.xy.y / METERS_out;
1451 c.xyzt.x = *
x * METERS_in;
1452 c.xyzt.y = *y * METERS_in;
1455 c = proj_trans(info_trans.pj, PJ_FWD, c);
1456 ok = proj_errno(info_trans.pj);
1458 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1460 *
x = c.lp.lam * RAD_TO_DEG;
1461 *y = c.lp.phi * RAD_TO_DEG;
1465 *
x = c.xy.x / METERS_out;
1466 *y = c.xy.y / METERS_out;
1469 proj_destroy(info_trans.pj);
1472 G_warning(_(
"proj_trans() failed: %d"), ok);
1478 METERS_in = info_in->meters;
1479 METERS_out = info_out->meters;
1481 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1482 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1483 u = (*x) / RAD_TO_DEG;
1484 v = (*y) / RAD_TO_DEG;
1485 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1486 *
x = u * RAD_TO_DEG;
1487 *y = v * RAD_TO_DEG;
1490 u = (*x) / RAD_TO_DEG;
1491 v = (*y) / RAD_TO_DEG;
1492 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1493 *
x = u / METERS_out;
1494 *y = v / METERS_out;
1498 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1501 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1502 *
x = u * RAD_TO_DEG;
1503 *y = v * RAD_TO_DEG;
1508 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1509 *
x = u / METERS_out;
1510 *y = v / METERS_out;
1514 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1544 const struct pj_info *info_in,
const struct pj_info *info_out)
1550 struct pj_info info_trans;
1557 METERS_in = info_in->meters;
1558 METERS_out = info_out->meters;
1561 h = G_malloc(
sizeof *h *
count);
1563 for (i = 0; i <
count; ++i)
1568 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1570 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1571 for (i = 0; i <
count; i++) {
1573 c.lpzt.lam =
x[i] / RAD_TO_DEG;
1574 c.lpzt.phi = y[i] / RAD_TO_DEG;
1576 c = proj_trans(info_trans.pj, PJ_FWD, c);
1577 if ((ok = proj_errno(info_trans.pj)) < 0)
1580 x[i] = c.lp.lam * RAD_TO_DEG;
1581 y[i] = c.lp.phi * RAD_TO_DEG;
1585 for (i = 0; i <
count; i++) {
1587 c.lpzt.lam =
x[i] / RAD_TO_DEG;
1588 c.lpzt.phi = y[i] / RAD_TO_DEG;
1590 c = proj_trans(info_trans.pj, PJ_FWD, c);
1591 if ((ok = proj_errno(info_trans.pj)) < 0)
1594 x[i] = c.xy.x / METERS_out;
1595 y[i] = c.xy.y / METERS_out;
1601 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1602 for (i = 0; i <
count; i++) {
1604 c.xyzt.x =
x[i] * METERS_in;
1605 c.xyzt.y = y[i] * METERS_in;
1607 c = proj_trans(info_trans.pj, PJ_FWD, c);
1608 if ((ok = proj_errno(info_trans.pj)) < 0)
1611 x[i] = c.lp.lam * RAD_TO_DEG;
1612 y[i] = c.lp.phi * RAD_TO_DEG;
1616 for (i = 0; i <
count; i++) {
1618 c.xyzt.x =
x[i] * METERS_in;
1619 c.xyzt.y = y[i] * METERS_in;
1621 c = proj_trans(info_trans.pj, PJ_FWD, c);
1622 if ((ok = proj_errno(info_trans.pj)) < 0)
1625 x[i] = c.xy.x / METERS_out;
1626 y[i] = c.xy.y / METERS_out;
1632 proj_destroy(info_trans.pj);
1635 G_warning(_(
"proj_trans() failed: %d"), ok);
1638 METERS_in = info_in->meters;
1639 METERS_out = info_out->meters;
1642 h = G_malloc(
sizeof *h *
count);
1644 for (i = 0; i <
count; ++i)
1648 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1649 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1651 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1656 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1661 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1663 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1668 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1676 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
void G_free(void *buf)
Free allocated memory.
int G_asprintf(char **out, const char *fmt,...)
int G_debug(int level, const char *msg,...)
Print debugging message.
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...
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 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 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_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 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...
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
void G_warning(const char *msg,...)
Print a warning message to stderr.
void G_get_set_window(struct Cell_head *window)
Get the current working window (region)
char * G_store_lower(const char *s)
Copy string to allocated memory and convert copied string to lower case.
char * G_store(const char *s)
Copy string to allocated memory.
char * G_store_upper(const char *s)
Copy string to allocated memory and convert copied string to upper case.