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;
45 #if PROJ_VERSION_MAJOR >= 6
46 int 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);
181 char *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");
273 PJ *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++) {
576 const char *area_of_use, *projstr;
578 PJ_PROJ_INFO pj_info;
581 op = proj_list_get(PJ_DEFAULT_CTX, op_list, i);
582 op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
585 G_warning(_(
"proj_normalize_for_visualization() failed for operation %d"),
593 pj_info = proj_pj_info(op);
594 proj_get_area_of_use(
NULL, op, &w, &s, &e, &n, &area_of_use);
597 if (pj_info.description) {
599 pj_info.description);
606 if (pj_info.accuracy > 0) {
611 #if PROJ_VERSION_NUM >= 6020000
612 str = proj_get_remarks(op);
617 str = proj_get_scope(op);
624 projstr = proj_as_proj_string(
NULL, op,
644 proj_list_destroy(op_list);
652 operation_ctx = proj_create_operation_factory_context(PJ_DEFAULT_CTX,
NULL);
653 proj_operation_factory_context_set_area_of_interest(PJ_DEFAULT_CTX,
657 proj_operation_factory_context_set_spatial_criterion(PJ_DEFAULT_CTX,
659 PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT);
660 proj_operation_factory_context_set_grid_availability_use(PJ_DEFAULT_CTX,
662 PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
671 op_list = proj_create_operations(PJ_DEFAULT_CTX,
675 proj_operation_factory_context_destroy(operation_ctx);
678 op_count_area = proj_list_get_count(op_list);
679 if (op_count_area == 0) {
681 info_trans->pj =
NULL;
683 else if (op_count_area == 1) {
684 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
690 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
693 proj_list_destroy(op_list);
696 G_debug(1,
"trying %s to %s", indef, outdef);
717 proj_destroy(out_pj);
719 if (info_trans->pj) {
723 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ%d",
726 projstr = proj_as_proj_string(
NULL, info_trans->pj,
729 info_trans->def =
G_store(projstr);
737 pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, info_trans->pj);
740 G_warning(_(
"proj_normalize_for_visualization() failed for '%s'"),
744 projstr = proj_as_proj_string(
NULL, pj_norm,
746 if (projstr && *projstr) {
747 proj_destroy(info_trans->pj);
748 info_trans->pj = pj_norm;
749 info_trans->def =
G_store(projstr);
752 proj_destroy(pj_norm);
753 G_warning(_(
"No PROJ definition for normalized version of '%s'"),
762 proj_destroy(info_trans->pj);
763 info_trans->pj =
NULL;
768 proj_area_destroy(pj_area);
777 if (info_trans->pj ==
NULL) {
778 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
783 info_trans->pj =
NULL;
784 info_trans->meters = 1.;
785 info_trans->zone = 0;
786 sprintf(info_trans->proj,
"pipeline");
789 if (info_trans->def) {
791 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
792 if (info_trans->pj ==
NULL) {
793 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
800 else if (info_out->pj ==
NULL) {
801 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
803 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
804 if (info_trans->pj ==
NULL) {
805 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
810 else if (info_in->def && info_out->pj && info_out->def) {
812 char *insrid =
NULL, *outsrid =
NULL;
816 if (strncmp(info_in->srid,
"EPSG", 4) == 0)
819 insrid =
G_store(info_in->srid);
822 if (info_out->pj && info_out->srid) {
823 if (strncmp(info_out->srid,
"EPSG", 4) == 0)
826 outsrid =
G_store(info_out->srid);
829 info_trans->pj =
NULL;
831 if (insrid && outsrid) {
835 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv +init=%s +step +init=%s",
839 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
845 if (info_trans->pj) {
846 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ5");
872 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s +step %s",
875 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
884 if (info_trans->pj ==
NULL) {
885 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
892 if (info_out->pj ==
NULL) {
894 G_warning(_(
"Unable to create latlong equivalent for '%s'"),
933 const struct pj_info *info_out,
934 const struct pj_info *info_trans,
int dir,
935 double *
x,
double *y,
double *z)
941 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
944 if (info_in->pj ==
NULL)
947 if (info_trans->pj ==
NULL)
950 in_deg2rad = out_rad2deg = 1;
953 METERS_in = info_in->meters;
954 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
955 #if PROJ_VERSION_MAJOR >= 6
959 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
964 METERS_out = info_out->meters;
965 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
966 #if PROJ_VERSION_MAJOR >= 6
970 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
982 METERS_out = info_in->meters;
983 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
984 #if PROJ_VERSION_MAJOR >= 6
988 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
993 METERS_in = info_out->meters;
994 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
995 #if PROJ_VERSION_MAJOR >= 6
999 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1014 c.lpzt.lam = (*x) / RAD_TO_DEG;
1015 c.lpzt.phi = (*y) / RAD_TO_DEG;
1028 c.xyzt.x = *
x * METERS_in;
1029 c.xyzt.y = *y * METERS_in;
1036 G_debug(1,
"c.xyzt.x: %g", c.xyzt.x);
1037 G_debug(1,
"c.xyzt.y: %g", c.xyzt.y);
1038 G_debug(1,
"c.xyzt.z: %g", c.xyzt.z);
1041 c = proj_trans(info_trans->pj, dir, c);
1042 ok = proj_errno(info_trans->pj);
1045 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1054 *
x = c.lp.lam * RAD_TO_DEG;
1055 *y = c.lp.phi * RAD_TO_DEG;
1066 *
x = c.xyzt.x / METERS_out;
1067 *y = c.xyzt.y / METERS_out;
1075 const struct pj_info *p_in, *p_out;
1077 if (info_out ==
NULL)
1080 if (dir == PJ_FWD) {
1089 METERS_in = p_in->meters;
1090 METERS_out = p_out->meters;
1095 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1096 u = (*x) / RAD_TO_DEG;
1097 v = (*y) / RAD_TO_DEG;
1104 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1107 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1111 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1112 *
x = u * RAD_TO_DEG;
1113 *y = v * RAD_TO_DEG;
1116 *
x = u / METERS_out;
1117 *y = v / METERS_out;
1153 const struct pj_info *info_out,
1154 const struct pj_info *info_trans,
int dir,
1155 double *
x,
double *y,
double *z,
int n)
1163 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1166 if (info_trans->pj ==
NULL)
1169 in_deg2rad = out_rad2deg = 1;
1170 if (dir == PJ_FWD) {
1172 METERS_in = info_in->meters;
1173 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
1174 #if PROJ_VERSION_MAJOR >= 6
1178 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1183 METERS_out = info_out->meters;
1184 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
1185 #if PROJ_VERSION_MAJOR >= 6
1189 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1201 METERS_out = info_in->meters;
1202 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
1203 #if PROJ_VERSION_MAJOR >= 6
1207 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1212 METERS_in = info_out->meters;
1213 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
1214 #if PROJ_VERSION_MAJOR >= 6
1218 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1230 z = G_malloc(
sizeof(
double) * n);
1232 for (i = 0; i < n; i++)
1247 for (i = 0; i < n; i++) {
1250 c.lpzt.lam = (*x) / RAD_TO_DEG;
1251 c.lpzt.phi = (*y) / RAD_TO_DEG;
1258 c = proj_trans(info_trans->pj, dir, c);
1259 if ((ok = proj_errno(info_trans->pj)) < 0)
1263 x[i] = c.lp.lam * RAD_TO_DEG;
1264 y[i] = c.lp.phi * RAD_TO_DEG;
1273 for (i = 0; i < n; i++) {
1276 c.lpzt.lam = (*x) / RAD_TO_DEG;
1277 c.lpzt.phi = (*y) / RAD_TO_DEG;
1284 c = proj_trans(info_trans->pj, dir, c);
1285 if ((ok = proj_errno(info_trans->pj)) < 0)
1289 x[i] = c.xy.x / METERS_out;
1290 y[i] = c.xy.y / METERS_out;
1297 for (i = 0; i < n; i++) {
1299 c.xyzt.x =
x[i] * METERS_in;
1300 c.xyzt.y = y[i] * METERS_in;
1302 c = proj_trans(info_trans->pj, dir, c);
1303 if ((ok = proj_errno(info_trans->pj)) < 0)
1307 x[i] = c.lp.lam * RAD_TO_DEG;
1308 y[i] = c.lp.phi * RAD_TO_DEG;
1317 for (i = 0; i < n; i++) {
1319 c.xyzt.x =
x[i] * METERS_in;
1320 c.xyzt.y = y[i] * METERS_in;
1322 c = proj_trans(info_trans->pj, dir, c);
1323 if ((ok = proj_errno(info_trans->pj)) < 0)
1326 x[i] = c.xy.x / METERS_out;
1327 y[i] = c.xy.y / METERS_out;
1335 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1339 const struct pj_info *p_in, *p_out;
1341 if (dir == PJ_FWD) {
1350 METERS_in = p_in->meters;
1351 METERS_out = p_out->meters;
1354 z = G_malloc(
sizeof(
double) * n);
1356 for (i = 0; i < n; ++i)
1360 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1361 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1363 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1368 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1373 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1375 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1380 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1388 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1416 const struct pj_info *info_in,
const struct pj_info *info_out)
1420 struct pj_info info_trans;
1427 METERS_in = info_in->meters;
1428 METERS_out = info_out->meters;
1430 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1432 c.lpzt.lam = (*x) / RAD_TO_DEG;
1433 c.lpzt.phi = (*y) / RAD_TO_DEG;
1436 c = proj_trans(info_trans.pj, PJ_FWD, c);
1437 ok = proj_errno(info_trans.pj);
1439 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1441 *
x = c.lp.lam * RAD_TO_DEG;
1442 *y = c.lp.phi * RAD_TO_DEG;
1446 *
x = c.xy.x / METERS_out;
1447 *y = c.xy.y / METERS_out;
1452 c.xyzt.x = *
x * METERS_in;
1453 c.xyzt.y = *y * METERS_in;
1456 c = proj_trans(info_trans.pj, PJ_FWD, c);
1457 ok = proj_errno(info_trans.pj);
1459 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1461 *
x = c.lp.lam * RAD_TO_DEG;
1462 *y = c.lp.phi * RAD_TO_DEG;
1466 *
x = c.xy.x / METERS_out;
1467 *y = c.xy.y / METERS_out;
1470 proj_destroy(info_trans.pj);
1473 G_warning(_(
"proj_trans() failed: %d"), ok);
1479 METERS_in = info_in->meters;
1480 METERS_out = info_out->meters;
1482 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1483 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1484 u = (*x) / RAD_TO_DEG;
1485 v = (*y) / RAD_TO_DEG;
1486 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1487 *
x = u * RAD_TO_DEG;
1488 *y = v * RAD_TO_DEG;
1491 u = (*x) / RAD_TO_DEG;
1492 v = (*y) / RAD_TO_DEG;
1493 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1494 *
x = u / METERS_out;
1495 *y = v / METERS_out;
1499 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1502 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1503 *
x = u * RAD_TO_DEG;
1504 *y = v * RAD_TO_DEG;
1509 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1510 *
x = u / METERS_out;
1511 *y = v / METERS_out;
1515 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1545 const struct pj_info *info_in,
const struct pj_info *info_out)
1551 struct pj_info info_trans;
1558 METERS_in = info_in->meters;
1559 METERS_out = info_out->meters;
1562 h = G_malloc(
sizeof *h *
count);
1564 for (i = 0; i <
count; ++i)
1569 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1571 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1572 for (i = 0; i <
count; i++) {
1574 c.lpzt.lam =
x[i] / RAD_TO_DEG;
1575 c.lpzt.phi = y[i] / RAD_TO_DEG;
1577 c = proj_trans(info_trans.pj, PJ_FWD, c);
1578 if ((ok = proj_errno(info_trans.pj)) < 0)
1581 x[i] = c.lp.lam * RAD_TO_DEG;
1582 y[i] = c.lp.phi * RAD_TO_DEG;
1586 for (i = 0; i <
count; i++) {
1588 c.lpzt.lam =
x[i] / RAD_TO_DEG;
1589 c.lpzt.phi = y[i] / RAD_TO_DEG;
1591 c = proj_trans(info_trans.pj, PJ_FWD, c);
1592 if ((ok = proj_errno(info_trans.pj)) < 0)
1595 x[i] = c.xy.x / METERS_out;
1596 y[i] = c.xy.y / METERS_out;
1602 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1603 for (i = 0; i <
count; i++) {
1605 c.xyzt.x =
x[i] * METERS_in;
1606 c.xyzt.y = y[i] * METERS_in;
1608 c = proj_trans(info_trans.pj, PJ_FWD, c);
1609 if ((ok = proj_errno(info_trans.pj)) < 0)
1612 x[i] = c.lp.lam * RAD_TO_DEG;
1613 y[i] = c.lp.phi * RAD_TO_DEG;
1617 for (i = 0; i <
count; i++) {
1619 c.xyzt.x =
x[i] * METERS_in;
1620 c.xyzt.y = y[i] * METERS_in;
1622 c = proj_trans(info_trans.pj, PJ_FWD, c);
1623 if ((ok = proj_errno(info_trans.pj)) < 0)
1626 x[i] = c.xy.x / METERS_out;
1627 y[i] = c.xy.y / METERS_out;
1633 proj_destroy(info_trans.pj);
1636 G_warning(_(
"proj_trans() failed: %d"), ok);
1639 METERS_in = info_in->meters;
1640 METERS_out = info_out->meters;
1643 h = G_malloc(
sizeof *h *
count);
1645 for (i = 0; i <
count; ++i)
1649 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1650 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1652 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1657 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1662 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1664 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1669 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1677 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(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.
char * G_store_lower(const char *s)
Copy string to allocated memory and convert copied string to lower case.