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);
78 proj_as_proj_string(
NULL, source_crs, PJ_PROJ_5,
NULL);
82 proj_destroy(source_crs);
86 projstr = proj_as_proj_string(
NULL, iproj->pj, PJ_PROJ_5,
NULL);
98 G_asprintf(&tproj.def,
"+proj=pipeline +step +inv %s +over", indef);
99 G_debug(1,
"get_pj_area() tproj.def: %s", tproj.def);
100 tproj.pj = proj_create(PJ_DEFAULT_CTX, tproj.def);
102 if (tproj.pj ==
NULL) {
103 G_warning(_(
"proj_create() failed for '%s'"), tproj.def);
106 proj_destroy(tproj.pj);
110 projstr = proj_as_proj_string(
NULL, tproj.pj, PJ_PROJ_5,
NULL);
111 if (projstr ==
NULL) {
112 G_warning(_(
"proj_create() failed for '%s'"), tproj.def);
115 proj_destroy(tproj.pj);
120 G_debug(1,
"proj_create() projstr '%s'", projstr);
126 estep = (window.east - window.west) / 21.;
127 nstep = (window.north - window.south) / 21.;
128 for (i = 0; i < 20; i++) {
129 x[i] = window.west + estep * (i + 1);
132 x[i + 20] = window.west + estep * (i + 1);
133 y[i + 20] = window.south;
135 x[i + 40] = window.west;
136 y[i + 40] = window.south + nstep * (i + 1);
138 x[i + 60] = window.east;
139 y[i + 60] = window.south + nstep * (i + 1);
142 y[80] = window.north;
144 y[81] = window.south;
146 y[82] = window.north;
148 y[83] = window.south;
149 x[84] = (window.west + window.east) / 2.;
150 y[84] = (window.north + window.south) / 2.;
154 proj_destroy(tproj.pj);
157 *xmin = *xmax =
x[84];
158 *ymin = *ymax = y[84];
159 for (i = 0; i < 84; i++) {
174 if (*xmin < -180 && *xmax < 180 && *xmin + 360 > *xmax) {
178 else if (*xmax > 180 && *xmin > -180 && *xmax - 360 < *xmin) {
183 G_debug(1,
"input window north: %.8f", window.north);
184 G_debug(1,
"input window south: %.8f", window.south);
185 G_debug(1,
"input window east: %.8f", window.east);
186 G_debug(1,
"input window west: %.8f", window.west);
188 G_debug(1,
"transformed xmin: %.8f", *xmin);
189 G_debug(1,
"transformed xmax: %.8f", *xmax);
190 G_debug(1,
"transformed ymin: %.8f", *ymin);
191 G_debug(1,
"transformed ymax: %.8f", *ymax);
197 if (fabs(*xmin) > 180) {
198 G_warning(_(
"Invalid west longitude %g"), *xmin);
201 if (fabs(*xmax) > 180) {
202 G_warning(_(
"Invalid east longitude %g"), *xmax);
205 if (fabs(*ymin) > 90) {
206 G_warning(_(
"Invalid south latitude %g"), *ymin);
209 if (fabs(*ymax) > 90) {
210 G_warning(_(
"Invalid north latitude %g"), *ymax);
214 G_warning(_(
"South %g is larger than north %g"), *ymin, *ymax);
218 G_debug(1,
"get_pj_area(): xmin %g, xmax %g, ymin %g, ymax %g", *xmin,
219 *xmax, *ymin, *ymax);
224 char *get_pj_type_string(PJ *pj)
226 char *pj_type =
NULL;
228 switch (proj_get_type(pj)) {
229 case PJ_TYPE_UNKNOWN:
232 case PJ_TYPE_ELLIPSOID:
235 case PJ_TYPE_PRIME_MERIDIAN:
238 case PJ_TYPE_GEODETIC_REFERENCE_FRAME:
239 G_asprintf(&pj_type,
"geodetic reference frame");
241 case PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME:
242 G_asprintf(&pj_type,
"dynamic geodetic reference frame");
244 case PJ_TYPE_VERTICAL_REFERENCE_FRAME:
245 G_asprintf(&pj_type,
"vertical reference frame");
247 case PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME:
248 G_asprintf(&pj_type,
"dynamic vertical reference frame");
250 case PJ_TYPE_DATUM_ENSEMBLE:
258 case PJ_TYPE_GEODETIC_CRS:
261 case PJ_TYPE_GEOCENTRIC_CRS:
267 case PJ_TYPE_GEOGRAPHIC_CRS:
270 case PJ_TYPE_GEOGRAPHIC_2D_CRS:
273 case PJ_TYPE_GEOGRAPHIC_3D_CRS:
276 case PJ_TYPE_VERTICAL_CRS:
279 case PJ_TYPE_PROJECTED_CRS:
282 case PJ_TYPE_COMPOUND_CRS:
285 case PJ_TYPE_TEMPORAL_CRS:
288 case PJ_TYPE_ENGINEERING_CRS:
291 case PJ_TYPE_BOUND_CRS:
294 case PJ_TYPE_OTHER_CRS:
297 case PJ_TYPE_CONVERSION:
300 case PJ_TYPE_TRANSFORMATION:
303 case PJ_TYPE_CONCATENATED_OPERATION:
304 G_asprintf(&pj_type,
"concatenated operation");
306 case PJ_TYPE_OTHER_COORDINATE_OPERATION:
307 G_asprintf(&pj_type,
"other coordinate operation");
317 PJ *get_pj_object(
const struct pj_info *in_gpj,
char **in_defstr)
325 G_debug(1,
"Trying SRID '%s' ...", in_gpj->srid);
326 in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->srid);
328 G_warning(_(
"Unrecognized SRID '%s'"), in_gpj->srid);
331 *in_defstr =
G_store(in_gpj->srid);
335 ((
struct pj_info *)in_gpj)->meters = 1;
338 if (!in_pj && in_gpj->wkt) {
339 G_debug(1,
"Trying WKT '%s' ...", in_gpj->wkt);
340 in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->wkt);
342 G_warning(_(
"Unrecognized WKT '%s'"), in_gpj->wkt);
345 *in_defstr =
G_store(in_gpj->wkt);
349 ((
struct pj_info *)in_gpj)->meters = 1;
352 if (!in_pj && in_gpj->pj) {
353 in_pj = proj_clone(PJ_DEFAULT_CTX, in_gpj->pj);
355 if (*in_defstr && !**in_defstr)
360 G_warning(_(
"Unable to create PROJ object"));
376 if (proj_get_type(in_pj) == PJ_TYPE_BOUND_CRS) {
380 source_crs = proj_get_source_crs(
NULL, in_pj);
384 if (*in_defstr && !**in_defstr)
435 const struct pj_info *info_out,
436 struct pj_info *info_trans)
438 if (info_in->pj ==
NULL)
441 if (info_in->def ==
NULL)
442 G_fatal_error(_(
"Input coordinate system definition is NULL"));
445 #if PROJ_VERSION_MAJOR >= 6
450 info_trans->pj =
NULL;
451 info_trans->meters = 1.;
452 info_trans->zone = 0;
453 sprintf(info_trans->proj,
"pipeline");
456 if (info_trans->def) {
461 if (!info_in->pj || !info_in->proj[0] || !info_out->pj ||
462 !info_out->proj[0]) {
464 "A custom pipeline requires input and output projection info"));
470 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
471 if (info_trans->pj ==
NULL) {
472 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
476 projstr = proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
477 if (projstr ==
NULL) {
478 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
488 info_trans->def =
G_store(projstr);
490 if (strstr(info_trans->def,
"axisswap")) {
492 _(
"The transformation pipeline contains an '%s' step. "
493 "Remove this step if easting and northing are swapped in "
498 G_debug(1,
"proj_create() pipeline: %s", info_trans->def);
503 ((
struct pj_info *)info_in)->meters = 1;
504 ((
struct pj_info *)info_out)->meters = 1;
509 else if (info_out->pj ==
NULL) {
510 const char *projstr =
NULL;
515 G_debug(1,
"ll equivalent definition: %s", indef);
520 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s", indef);
521 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
522 if (info_trans->pj ==
NULL) {
523 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
528 projstr = proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
529 if (projstr ==
NULL) {
530 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
538 else if (info_in->def && info_out->pj && info_out->def) {
540 char *insrid =
NULL, *outsrid =
NULL;
542 PJ_OBJ_LIST *op_list;
543 PJ_OPERATION_FACTORY_CONTEXT *operation_ctx;
544 PJ_AREA *pj_area =
NULL;
545 double xmin, xmax, ymin, ymax;
546 int op_count = 0, op_count_area = 0;
552 if (get_pj_area(info_in, &xmin, &xmax, &ymin, &ymax)) {
553 pj_area = proj_area_create();
554 proj_area_set_bbox(pj_area, xmin, ymin, xmax, ymax);
557 G_warning(_(
"Unable to determine area of interest for '%s'"),
563 G_debug(1,
"source proj string: %s", info_in->def);
564 G_debug(1,
"source type: %s", get_pj_type_string(info_in->pj));
568 if (strncmp(info_in->srid,
"epsg", 4) == 0) {
571 ((
struct pj_info *)info_in)->srid = insrid;
576 in_pj = get_pj_object(info_in, &indef);
577 if (in_pj ==
NULL || indef ==
NULL) {
578 G_warning(_(
"Input CRS not available for '%s'"), info_in->def);
582 G_debug(1,
"Input CRS definition: %s", indef);
584 G_debug(1,
"target proj string: %s", info_out->def);
585 G_debug(1,
"target type: %s", get_pj_type_string(info_out->pj));
588 if (info_out->srid) {
589 if (strncmp(info_out->srid,
"epsg", 4) == 0) {
592 ((
struct pj_info *)info_out)->srid = outsrid;
597 out_pj = get_pj_object(info_out, &outdef);
598 if (out_pj ==
NULL || outdef ==
NULL) {
599 G_warning(_(
"Output CRS not available for '%s'"), info_out->def);
603 G_debug(1,
"Output CRS definition: %s", outdef);
608 proj_create_operation_factory_context(PJ_DEFAULT_CTX,
NULL);
615 op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
617 proj_operation_factory_context_destroy(operation_ctx);
621 op_count = proj_list_get_count(op_list);
627 for (i = 0; i < op_count; i++) {
628 const char *area_of_use, *projstr;
630 PJ_PROJ_INFO pj_info;
633 op = proj_list_get(PJ_DEFAULT_CTX, op_list, i);
634 op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
637 G_warning(_(
"proj_normalize_for_visualization() failed for "
646 pj_info = proj_pj_info(op);
647 proj_get_area_of_use(
NULL, op, &w, &s, &e, &n, &area_of_use);
650 if (pj_info.description) {
652 pj_info.description);
658 if (pj_info.accuracy > 0) {
663 #if PROJ_VERSION_NUM >= 6020000
664 const char *str = proj_get_remarks(op);
670 str = proj_get_scope(op);
677 projstr = proj_as_proj_string(
NULL, op, PJ_PROJ_5,
NULL);
691 "with the %s option"),
697 proj_list_destroy(op_list);
705 proj_create_operation_factory_context(PJ_DEFAULT_CTX,
NULL);
706 proj_operation_factory_context_set_area_of_interest(
707 PJ_DEFAULT_CTX, operation_ctx, xmin, ymin, xmax, ymax);
708 proj_operation_factory_context_set_spatial_criterion(
709 PJ_DEFAULT_CTX, operation_ctx,
710 PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT);
711 proj_operation_factory_context_set_grid_availability_use(
712 PJ_DEFAULT_CTX, operation_ctx,
713 PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
722 op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
724 proj_operation_factory_context_destroy(operation_ctx);
727 op_count_area = proj_list_get_count(op_list);
728 if (op_count_area == 0) {
730 info_trans->pj =
NULL;
732 else if (op_count_area == 1) {
733 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
739 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
742 proj_list_destroy(op_list);
768 proj_destroy(out_pj);
770 if (info_trans->pj) {
774 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ%d",
778 proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
780 info_trans->def =
G_store(projstr);
788 pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX,
793 _(
"proj_normalize_for_visualization() failed for '%s'"),
798 proj_as_proj_string(
NULL, pj_norm, PJ_PROJ_5,
NULL);
799 if (projstr && *projstr) {
800 proj_destroy(info_trans->pj);
801 info_trans->pj = pj_norm;
802 info_trans->def =
G_store(projstr);
805 proj_destroy(pj_norm);
806 G_warning(_(
"No PROJ definition for normalized version "
816 proj_destroy(info_trans->pj);
817 info_trans->pj =
NULL;
822 proj_area_destroy(pj_area);
831 if (info_trans->pj ==
NULL) {
832 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
837 info_trans->pj =
NULL;
838 info_trans->meters = 1.;
839 info_trans->zone = 0;
840 sprintf(info_trans->proj,
"pipeline");
843 if (info_trans->def) {
845 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
846 if (info_trans->pj ==
NULL) {
847 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
854 else if (info_out->pj ==
NULL) {
855 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
857 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
858 if (info_trans->pj ==
NULL) {
859 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
864 else if (info_in->def && info_out->pj && info_out->def) {
866 char *insrid =
NULL, *outsrid =
NULL;
870 if (strncmp(info_in->srid,
"EPSG", 4) == 0)
873 insrid =
G_store(info_in->srid);
876 if (info_out->pj && info_out->srid) {
877 if (strncmp(info_out->srid,
"EPSG", 4) == 0)
880 outsrid =
G_store(info_out->srid);
883 info_trans->pj =
NULL;
885 if (insrid && outsrid) {
890 "+proj=pipeline +step +inv +init=%s +step +init=%s",
895 proj_create_crs_to_crs(PJ_DEFAULT_CTX, indef, outdef,
NULL);
898 if (info_trans->pj) {
899 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ5");
926 "+proj=pipeline +step +inv %s +step %s", indef, outdef);
928 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
937 if (info_trans->pj ==
NULL) {
938 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
945 if (info_out->pj ==
NULL) {
947 G_warning(_(
"Unable to create latlong equivalent for '%s'"),
985 int GPJ_transform(
const struct pj_info *info_in,
const struct pj_info *info_out,
986 const struct pj_info *info_trans,
int dir,
double *
x,
987 double *y,
double *z)
993 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
996 if (info_in->pj ==
NULL)
999 if (info_trans->pj ==
NULL)
1002 in_deg2rad = out_rad2deg = 1;
1003 if (dir == PJ_FWD) {
1005 METERS_in = info_in->meters;
1006 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
1007 #if PROJ_VERSION_MAJOR >= 6
1011 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1016 METERS_out = info_out->meters;
1017 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
1018 #if PROJ_VERSION_MAJOR >= 6
1022 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1034 METERS_out = info_in->meters;
1035 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
1036 #if PROJ_VERSION_MAJOR >= 6
1040 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1045 METERS_in = info_out->meters;
1046 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
1047 #if PROJ_VERSION_MAJOR >= 6
1051 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1066 c.lpzt.lam = (*x) / RAD_TO_DEG;
1067 c.lpzt.phi = (*y) / RAD_TO_DEG;
1080 c.xyzt.x = *
x * METERS_in;
1081 c.xyzt.y = *y * METERS_in;
1088 G_debug(1,
"c.xyzt.x: %g", c.xyzt.x);
1089 G_debug(1,
"c.xyzt.y: %g", c.xyzt.y);
1090 G_debug(1,
"c.xyzt.z: %g", c.xyzt.z);
1093 c = proj_trans(info_trans->pj, dir, c);
1094 ok = proj_errno(info_trans->pj);
1097 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1106 *
x = c.lp.lam * RAD_TO_DEG;
1107 *y = c.lp.phi * RAD_TO_DEG;
1118 *
x = c.xyzt.x / METERS_out;
1119 *y = c.xyzt.y / METERS_out;
1127 const struct pj_info *p_in, *p_out;
1129 if (info_out ==
NULL)
1132 if (dir == PJ_FWD) {
1141 METERS_in = p_in->meters;
1142 METERS_out = p_out->meters;
1147 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1148 u = (*x) / RAD_TO_DEG;
1149 v = (*y) / RAD_TO_DEG;
1156 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1159 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1163 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1164 *
x = u * RAD_TO_DEG;
1165 *y = v * RAD_TO_DEG;
1168 *
x = u / METERS_out;
1169 *y = v / METERS_out;
1205 const struct pj_info *info_out,
1206 const struct pj_info *info_trans,
int dir,
double *
x,
1207 double *y,
double *z,
int n)
1215 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1218 if (info_trans->pj ==
NULL)
1221 in_deg2rad = out_rad2deg = 1;
1222 if (dir == PJ_FWD) {
1224 METERS_in = info_in->meters;
1225 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
1226 #if PROJ_VERSION_MAJOR >= 6
1230 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1235 METERS_out = info_out->meters;
1236 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
1237 #if PROJ_VERSION_MAJOR >= 6
1241 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1253 METERS_out = info_in->meters;
1254 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
1255 #if PROJ_VERSION_MAJOR >= 6
1259 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1264 METERS_in = info_out->meters;
1265 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
1266 #if PROJ_VERSION_MAJOR >= 6
1270 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1282 z = G_malloc(
sizeof(
double) * n);
1284 for (i = 0; i < n; i++)
1299 for (i = 0; i < n; i++) {
1302 c.lpzt.lam = (*x) / RAD_TO_DEG;
1303 c.lpzt.phi = (*y) / RAD_TO_DEG;
1310 c = proj_trans(info_trans->pj, dir, c);
1311 if ((ok = proj_errno(info_trans->pj)) < 0)
1315 x[i] = c.lp.lam * RAD_TO_DEG;
1316 y[i] = c.lp.phi * RAD_TO_DEG;
1325 for (i = 0; i < n; i++) {
1328 c.lpzt.lam = (*x) / RAD_TO_DEG;
1329 c.lpzt.phi = (*y) / RAD_TO_DEG;
1336 c = proj_trans(info_trans->pj, dir, c);
1337 if ((ok = proj_errno(info_trans->pj)) < 0)
1341 x[i] = c.xy.x / METERS_out;
1342 y[i] = c.xy.y / METERS_out;
1349 for (i = 0; i < n; i++) {
1351 c.xyzt.x =
x[i] * METERS_in;
1352 c.xyzt.y = y[i] * METERS_in;
1354 c = proj_trans(info_trans->pj, dir, c);
1355 if ((ok = proj_errno(info_trans->pj)) < 0)
1359 x[i] = c.lp.lam * RAD_TO_DEG;
1360 y[i] = c.lp.phi * RAD_TO_DEG;
1369 for (i = 0; i < n; i++) {
1371 c.xyzt.x =
x[i] * METERS_in;
1372 c.xyzt.y = y[i] * METERS_in;
1374 c = proj_trans(info_trans->pj, dir, c);
1375 if ((ok = proj_errno(info_trans->pj)) < 0)
1378 x[i] = c.xy.x / METERS_out;
1379 y[i] = c.xy.y / METERS_out;
1387 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1391 const struct pj_info *p_in, *p_out;
1393 if (dir == PJ_FWD) {
1402 METERS_in = p_in->meters;
1403 METERS_out = p_out->meters;
1406 z = G_malloc(
sizeof(
double) * n);
1408 for (i = 0; i < n; ++i)
1412 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1413 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1415 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1420 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1425 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1427 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1432 ok = pj_transform(info_in->pj, info_out->pj, n, 1,
x, y, z);
1440 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1468 const struct pj_info *info_out)
1473 struct pj_info info_trans;
1480 METERS_in = info_in->meters;
1481 METERS_out = info_out->meters;
1483 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1485 c.lpzt.lam = (*x) / RAD_TO_DEG;
1486 c.lpzt.phi = (*y) / RAD_TO_DEG;
1489 c = proj_trans(info_trans.pj, PJ_FWD, c);
1490 ok = proj_errno(info_trans.pj);
1492 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1494 *
x = c.lp.lam * RAD_TO_DEG;
1495 *y = c.lp.phi * RAD_TO_DEG;
1499 *
x = c.xy.x / METERS_out;
1500 *y = c.xy.y / METERS_out;
1505 c.xyzt.x = *
x * METERS_in;
1506 c.xyzt.y = *y * METERS_in;
1509 c = proj_trans(info_trans.pj, PJ_FWD, c);
1510 ok = proj_errno(info_trans.pj);
1512 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1514 *
x = c.lp.lam * RAD_TO_DEG;
1515 *y = c.lp.phi * RAD_TO_DEG;
1519 *
x = c.xy.x / METERS_out;
1520 *y = c.xy.y / METERS_out;
1523 proj_destroy(info_trans.pj);
1526 G_warning(_(
"proj_trans() failed: %d"), ok);
1532 METERS_in = info_in->meters;
1533 METERS_out = info_out->meters;
1535 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1536 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1537 u = (*x) / RAD_TO_DEG;
1538 v = (*y) / RAD_TO_DEG;
1539 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1540 *
x = u * RAD_TO_DEG;
1541 *y = v * RAD_TO_DEG;
1544 u = (*x) / RAD_TO_DEG;
1545 v = (*y) / RAD_TO_DEG;
1546 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1547 *
x = u / METERS_out;
1548 *y = v / METERS_out;
1552 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1555 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1556 *
x = u * RAD_TO_DEG;
1557 *y = v * RAD_TO_DEG;
1562 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1563 *
x = u / METERS_out;
1564 *y = v / METERS_out;
1568 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1598 const struct pj_info *info_in,
1599 const struct pj_info *info_out)
1606 struct pj_info info_trans;
1613 METERS_in = info_in->meters;
1614 METERS_out = info_out->meters;
1617 h = G_malloc(
sizeof *h *
count);
1619 for (i = 0; i <
count; ++i)
1624 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1626 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1627 for (i = 0; i <
count; i++) {
1629 c.lpzt.lam =
x[i] / RAD_TO_DEG;
1630 c.lpzt.phi = y[i] / RAD_TO_DEG;
1632 c = proj_trans(info_trans.pj, PJ_FWD, c);
1633 if ((ok = proj_errno(info_trans.pj)) < 0)
1636 x[i] = c.lp.lam * RAD_TO_DEG;
1637 y[i] = c.lp.phi * RAD_TO_DEG;
1641 for (i = 0; i <
count; i++) {
1643 c.lpzt.lam =
x[i] / RAD_TO_DEG;
1644 c.lpzt.phi = y[i] / RAD_TO_DEG;
1646 c = proj_trans(info_trans.pj, PJ_FWD, c);
1647 if ((ok = proj_errno(info_trans.pj)) < 0)
1650 x[i] = c.xy.x / METERS_out;
1651 y[i] = c.xy.y / METERS_out;
1657 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1658 for (i = 0; i <
count; i++) {
1660 c.xyzt.x =
x[i] * METERS_in;
1661 c.xyzt.y = y[i] * METERS_in;
1663 c = proj_trans(info_trans.pj, PJ_FWD, c);
1664 if ((ok = proj_errno(info_trans.pj)) < 0)
1667 x[i] = c.lp.lam * RAD_TO_DEG;
1668 y[i] = c.lp.phi * RAD_TO_DEG;
1672 for (i = 0; i <
count; i++) {
1674 c.xyzt.x =
x[i] * METERS_in;
1675 c.xyzt.y = y[i] * METERS_in;
1677 c = proj_trans(info_trans.pj, PJ_FWD, c);
1678 if ((ok = proj_errno(info_trans.pj)) < 0)
1681 x[i] = c.xy.x / METERS_out;
1682 y[i] = c.xy.y / METERS_out;
1688 proj_destroy(info_trans.pj);
1691 G_warning(_(
"proj_trans() failed: %d"), ok);
1694 METERS_in = info_in->meters;
1695 METERS_out = info_out->meters;
1698 h = G_malloc(
sizeof *h *
count);
1700 for (i = 0; i <
count; ++i)
1704 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1705 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1707 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1712 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1717 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1719 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1724 ok = pj_transform(info_in->pj, info_out->pj,
count, 1,
x, y, h);
1732 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.