21 #include <grass/gis.h> 22 #include <grass/gprojects.h> 23 #include <grass/glocale.h> 27 #ifdef PROJ_VERSION_NUM 28 #undef PROJ_VERSION_NUM 30 #define PROJ_VERSION_NUM ((PROJ_VERSION_MAJOR)*1000000+(PROJ_VERSION_MINOR)*10000+(PROJ_VERSION_PATCH)*100) 34 #define MULTIPLY_LOOP(x,y,c,m) \ 36 for (i = 0; i < c; ++i) {\ 42 #define DIVIDE_LOOP(x,y,c,m) \ 44 for (i = 0; i < c; ++i) {\ 50 static double METERS_in = 1.0, METERS_out = 1.0;
53 #if PROJ_VERSION_MAJOR >= 6 54 int get_pj_area(
double *xmin,
double *xmax,
double *ymin,
double *ymax)
56 struct Cell_head window;
65 if (window.proj != PROJECTION_LL) {
70 const char *projstr =
NULL;
73 struct Key_Value *in_proj_info, *in_unit_info;
74 struct pj_info iproj, oproj, tproj;
79 G_warning(_(
"Can't get projection info of current location"));
85 G_warning(_(
"Can't get projection units of current location"));
90 if (
pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0) {
91 G_fatal_error(_(
"Can't get projection key values of current location"));
102 source_crs = proj_get_source_crs(
NULL, iproj.pj);
104 projstr = proj_as_proj_string(
NULL, source_crs, PJ_PROJ_5,
NULL);
108 proj_destroy(iproj.pj);
109 iproj.pj = source_crs;
116 G_asprintf(&tproj.def,
"+proj=pipeline +step +inv %s",
118 tproj.pj = proj_create(PJ_DEFAULT_CTX, tproj.def);
119 if (tproj.pj ==
NULL) {
120 G_warning(_(
"proj_create() failed for '%s'"), tproj.def);
123 proj_destroy(tproj.pj);
127 projstr = proj_as_proj_string(
NULL, tproj.pj, PJ_PROJ_5,
NULL);
128 if (projstr ==
NULL) {
129 G_warning(_(
"proj_create() failed for '%s'"), tproj.def);
132 proj_destroy(tproj.pj);
138 estep = (window.west + window.east) / 21.;
139 nstep = (window.north + window.south) / 21.;
140 for (i = 0; i < 20; i++) {
141 x[i] = window.west + estep * (i + 1);
144 x[i + 20] = window.west + estep * (i + 1);
145 y[i + 20] = window.south;
147 x[i + 40] = window.west;
148 y[i + 40] = window.south + nstep * (i + 1);
150 x[i + 60] = window.east;
151 y[i + 60] = window.south + nstep * (i + 1);
154 y[80] = window.north;
156 y[81] = window.south;
158 y[82] = window.north;
160 y[83] = window.south;
161 x[84] = (window.west + window.east) / 2.;
162 y[84] = (window.north + window.south) / 2.;
166 proj_destroy(tproj.pj);
168 *xmin = *xmax = x[84];
169 *ymin = *ymax = y[84];
170 for (i = 0; i < 84; i++) {
181 G_debug(1,
"get_pj_area(): xmin %g, xmax %g, ymin %g, ymax %g",
182 *xmin, *xmax, *ymin, *ymax);
187 char *get_pj_type_string(PJ *pj)
189 char *pj_type =
NULL;
191 switch (proj_get_type(pj)) {
192 case PJ_TYPE_UNKNOWN:
195 case PJ_TYPE_ELLIPSOID:
198 case PJ_TYPE_PRIME_MERIDIAN:
201 case PJ_TYPE_GEODETIC_REFERENCE_FRAME:
202 G_asprintf(&pj_type,
"geodetic reference frame");
204 case PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME:
205 G_asprintf(&pj_type,
"dynamic geodetic reference frame");
207 case PJ_TYPE_VERTICAL_REFERENCE_FRAME:
208 G_asprintf(&pj_type,
"vertical reference frame");
210 case PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME:
211 G_asprintf(&pj_type,
"dynamic vertical reference frame");
213 case PJ_TYPE_DATUM_ENSEMBLE:
220 case PJ_TYPE_GEODETIC_CRS:
223 case PJ_TYPE_GEOCENTRIC_CRS:
228 case PJ_TYPE_GEOGRAPHIC_CRS:
231 case PJ_TYPE_GEOGRAPHIC_2D_CRS:
234 case PJ_TYPE_GEOGRAPHIC_3D_CRS:
237 case PJ_TYPE_VERTICAL_CRS:
240 case PJ_TYPE_PROJECTED_CRS:
243 case PJ_TYPE_COMPOUND_CRS:
246 case PJ_TYPE_TEMPORAL_CRS:
249 case PJ_TYPE_ENGINEERING_CRS:
252 case PJ_TYPE_BOUND_CRS:
255 case PJ_TYPE_OTHER_CRS:
258 case PJ_TYPE_CONVERSION:
261 case PJ_TYPE_TRANSFORMATION:
264 case PJ_TYPE_CONCATENATED_OPERATION:
265 G_asprintf(&pj_type,
"concatenated operation");
267 case PJ_TYPE_OTHER_COORDINATE_OPERATION:
268 G_asprintf(&pj_type,
"other coordinate operation");
316 const struct pj_info *info_out,
317 struct pj_info *info_trans)
319 if (info_in->pj ==
NULL)
322 if (info_in->def ==
NULL)
323 G_fatal_error(_(
"Input coordinate system definition is NULL"));
326 #if PROJ_VERSION_MAJOR >= 6 331 info_trans->pj =
NULL;
332 info_trans->meters = 1.;
333 info_trans->zone = 0;
334 sprintf(info_trans->proj,
"pipeline");
337 if (info_trans->def) {
341 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
342 if (info_trans->pj ==
NULL) {
343 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
347 projstr = proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
348 if (projstr ==
NULL) {
349 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
359 info_trans->def =
G_store(projstr);
361 if (strstr(info_trans->def,
"axisswap")) {
362 G_warning(_(
"The transformation pipeline contains an '%s' step. " 363 "Remove this step if easting and northing are swapped in the output."),
367 G_debug(1,
"proj_create() pipeline: %s", info_trans->def);
372 ((
struct pj_info *)info_in)->meters = 1;
373 ((
struct pj_info *)info_out)->meters = 1;
378 else if (info_out->pj ==
NULL) {
379 const char *projstr =
NULL;
394 source_crs = proj_get_source_crs(
NULL, info_in->pj);
396 projstr = proj_as_proj_string(
NULL, source_crs, PJ_PROJ_5,
NULL);
399 proj_destroy(source_crs);
407 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
409 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
410 if (info_trans->pj ==
NULL) {
411 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
416 projstr = proj_as_proj_string(
NULL, info_trans->pj, PJ_PROJ_5,
NULL);
417 if (projstr ==
NULL) {
418 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
426 else if (info_in->def && info_out->pj && info_out->def) {
428 char *indefcrs =
NULL, *outdefcrs =
NULL;
429 char *insrid =
NULL, *outsrid =
NULL;
430 int use_insrid = 0, use_outsrid = 0;
431 PJ *source_crs, *target_crs;
432 PJ_AREA *pj_area =
NULL;
433 double xmin, xmax, ymin, ymax;
438 G_debug(1,
"source proj string: %s", info_in->def);
439 G_debug(1,
"source type: %s", get_pj_type_string(info_in->pj));
440 indefcrs = info_in->def;
441 source_crs = proj_get_source_crs(
NULL, info_in->pj);
445 projstr = proj_as_proj_string(
NULL, source_crs, PJ_PROJ_5,
NULL);
448 G_debug(1,
"Input CRS definition converted from '%s' to '%s'",
449 info_in->def, indefcrs);
451 proj_destroy(source_crs);
455 G_debug(1,
"target proj string: %s", info_out->def);
456 G_debug(1,
"target type: %s", get_pj_type_string(info_out->pj));
457 outdefcrs = info_out->def;
458 target_crs = proj_get_source_crs(
NULL, info_out->pj);
462 projstr = proj_as_proj_string(
NULL, target_crs, PJ_PROJ_5,
NULL);
465 G_debug(1,
"Output CRS definition converted from '%s' to '%s'",
466 info_out->def, outdefcrs);
468 proj_destroy(target_crs);
474 if (strncmp(info_in->srid,
"epsg", 4) == 0)
477 insrid =
G_store(info_in->srid);
480 if (info_out->srid) {
481 if (strncmp(info_out->srid,
"epsg", 4) == 0)
484 outsrid =
G_store(info_out->srid);
494 G_debug(1,
"Input CRS definition: %s", indef);
503 G_debug(1,
"Output CRS definition: %s", outdef);
506 source_crs = proj_create(
NULL, indef);
507 target_crs = proj_create(
NULL, outdef);
510 if (get_pj_area(&xmin, &xmax, &ymin, &ymax)) {
511 pj_area = proj_area_create();
512 proj_area_set_bbox(pj_area, xmin, ymin, xmax, ymax);
515 if (source_crs && target_crs) {
516 PJ_OPERATION_FACTORY_CONTEXT *operation_ctx;
517 PJ_OBJ_LIST *op_list;
519 operation_ctx = proj_create_operation_factory_context(
NULL,
NULL);
525 op_list = proj_create_operations(
NULL,
532 op_count = proj_list_get_count(op_list);
536 G_warning(_(
"Found %d possible transformations"), op_count);
537 for (i = 0; i < op_count; i++) {
539 const char *area_of_use, *projstr;
541 PJ_PROJ_INFO pj_info;
544 op = proj_list_get(
NULL, op_list, i);
545 op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
548 G_warning(_(
"proj_normalize_for_visualization() failed for operation %d"),
556 projstr = proj_as_proj_string(
NULL, op,
558 pj_info = proj_pj_info(op);
559 proj_get_area_of_use(
NULL, op, &w, &s, &e, &n, &area_of_use);
564 pj_info.description);
568 if (pj_info.accuracy > 0) {
573 #if PROJ_VERSION_NUM >= 6020000 574 str = proj_get_remarks(op);
579 str = proj_get_scope(op);
602 proj_list_destroy(op_list);
603 proj_operation_factory_context_destroy(operation_ctx);
606 proj_destroy(source_crs);
608 proj_destroy(target_crs);
611 G_debug(1,
"trying %s to %s", indef, outdef);
613 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
618 if (info_trans->pj) {
619 const char *projstr =
NULL;
621 projstr = proj_as_proj_string(
NULL, info_trans->pj,
623 if (projstr ==
NULL) {
624 G_debug(1,
"proj_create_crs_to_crs() failed with PROJ%d for input \"%s\", output \"%s\"",
625 PROJ_VERSION_MAJOR, indef, outdef);
630 G_debug(1,
"trying %s to %s", indef, outdef);
633 proj_destroy(info_trans->pj);
634 info_trans->pj =
NULL;
635 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
645 ((
struct pj_info *)info_in)->meters = 1;
648 ((
struct pj_info *)info_out)->meters = 1;
653 if (info_trans->pj) {
657 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ%d",
660 projstr = proj_as_proj_string(
NULL, info_trans->pj,
663 info_trans->def =
G_store(projstr);
671 pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, info_trans->pj);
674 G_warning(_(
"proj_normalize_for_visualization() failed for '%s'"),
678 proj_destroy(info_trans->pj);
679 info_trans->pj = pj_norm;
680 projstr = proj_as_proj_string(
NULL, info_trans->pj,
682 info_trans->def =
G_store(projstr);
690 G_debug(1,
"proj_create_crs_to_crs() pipeline: %s", info_trans->def);
693 proj_destroy(info_trans->pj);
694 info_trans->pj =
NULL;
698 if (info_trans->pj ==
NULL) {
699 G_debug(1,
"proj_create_crs_to_crs() failed with PROJ%d for input \"%s\", output \"%s\"",
700 PROJ_VERSION_MAJOR, indef, outdef);
702 G_warning(
"GPJ_init_transform(): falling back to proj_create()");
713 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s +step %s",
716 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
720 proj_area_destroy(pj_area);
729 if (info_trans->pj ==
NULL) {
730 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
735 info_trans->pj =
NULL;
736 info_trans->meters = 1.;
737 info_trans->zone = 0;
738 sprintf(info_trans->proj,
"pipeline");
741 if (info_trans->def) {
743 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
744 if (info_trans->pj ==
NULL) {
745 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
752 else if (info_out->pj ==
NULL) {
753 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s",
755 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
756 if (info_trans->pj ==
NULL) {
757 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
762 else if (info_in->def && info_out->pj && info_out->def) {
764 char *insrid =
NULL, *outsrid =
NULL;
768 if (strncmp(info_in->srid,
"EPSG", 4) == 0)
771 insrid =
G_store(info_in->srid);
774 if (info_out->pj && info_out->srid) {
775 if (strncmp(info_out->srid,
"EPSG", 4) == 0)
778 outsrid =
G_store(info_out->srid);
781 info_trans->pj =
NULL;
783 if (insrid && outsrid) {
787 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv +init=%s +step +init=%s",
791 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
797 if (info_trans->pj) {
798 G_debug(1,
"proj_create_crs_to_crs() succeeded with PROJ5");
824 G_asprintf(&(info_trans->def),
"+proj=pipeline +step +inv %s +step %s",
827 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
836 if (info_trans->pj ==
NULL) {
837 G_warning(_(
"proj_create() failed for '%s'"), info_trans->def);
844 if (info_out->pj ==
NULL) {
846 G_warning(_(
"Unable to create latlong equivalent for '%s'"),
885 const struct pj_info *info_out,
886 const struct pj_info *info_trans,
int dir,
887 double *
x,
double *y,
double *z)
893 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
896 if (info_in->pj ==
NULL)
899 if (info_trans->pj ==
NULL)
902 in_deg2rad = out_rad2deg = 1;
905 METERS_in = info_in->meters;
906 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
907 #if PROJ_VERSION_MAJOR >= 6 911 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
916 METERS_out = info_out->meters;
917 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
918 #if PROJ_VERSION_MAJOR >= 6 922 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
934 METERS_out = info_in->meters;
935 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
936 #if PROJ_VERSION_MAJOR >= 6 940 if (out_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
945 METERS_in = info_out->meters;
946 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
947 #if PROJ_VERSION_MAJOR >= 6 951 if (in_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
966 c.lpzt.lam = (*x) / RAD_TO_DEG;
967 c.lpzt.phi = (*y) / RAD_TO_DEG;
980 c.xyzt.x = *x * METERS_in;
981 c.xyzt.y = *y * METERS_in;
989 c = proj_trans(info_trans->pj, dir, c);
990 ok = proj_errno(info_trans->pj);
993 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1002 *x = c.lp.lam * RAD_TO_DEG;
1003 *y = c.lp.phi * RAD_TO_DEG;
1014 *x = c.xyzt.x / METERS_out;
1015 *y = c.xyzt.y / METERS_out;
1023 const struct pj_info *p_in, *p_out;
1025 if (info_out ==
NULL)
1028 if (dir == PJ_FWD) {
1037 METERS_in = p_in->meters;
1038 METERS_out = p_out->meters;
1043 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1044 u = (*x) / RAD_TO_DEG;
1045 v = (*y) / RAD_TO_DEG;
1052 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1055 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1059 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1060 *x = u * RAD_TO_DEG;
1061 *y = v * RAD_TO_DEG;
1064 *x = u / METERS_out;
1065 *y = v / METERS_out;
1101 const struct pj_info *info_out,
1102 const struct pj_info *info_trans,
int dir,
1103 double *
x,
double *y,
double *z,
int n)
1111 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1114 if (info_trans->pj ==
NULL)
1117 in_deg2rad = out_rad2deg = 1;
1118 if (dir == PJ_FWD) {
1120 METERS_in = info_in->meters;
1121 in_is_ll = !strncmp(info_in->proj,
"ll", 2);
1122 #if PROJ_VERSION_MAJOR >= 6 1126 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1131 METERS_out = info_out->meters;
1132 out_is_ll = !strncmp(info_out->proj,
"ll", 2);
1133 #if PROJ_VERSION_MAJOR >= 6 1137 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1149 METERS_out = info_in->meters;
1150 out_is_ll = !strncmp(info_in->proj,
"ll", 2);
1151 #if PROJ_VERSION_MAJOR >= 6 1155 if (out_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1160 METERS_in = info_out->meters;
1161 in_is_ll = !strncmp(info_out->proj,
"ll", 2);
1162 #if PROJ_VERSION_MAJOR >= 6 1166 if (in_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1178 z = G_malloc(
sizeof(
double) * n);
1180 for (i = 0; i < n; i++)
1195 for (i = 0; i < n; i++) {
1198 c.lpzt.lam = (*x) / RAD_TO_DEG;
1199 c.lpzt.phi = (*y) / RAD_TO_DEG;
1206 c = proj_trans(info_trans->pj, dir, c);
1207 if ((ok = proj_errno(info_trans->pj)) < 0)
1211 x[i] = c.lp.lam * RAD_TO_DEG;
1212 y[i] = c.lp.phi * RAD_TO_DEG;
1221 for (i = 0; i < n; i++) {
1224 c.lpzt.lam = (*x) / RAD_TO_DEG;
1225 c.lpzt.phi = (*y) / RAD_TO_DEG;
1232 c = proj_trans(info_trans->pj, dir, c);
1233 if ((ok = proj_errno(info_trans->pj)) < 0)
1237 x[i] = c.xy.x / METERS_out;
1238 y[i] = c.xy.y / METERS_out;
1245 for (i = 0; i < n; i++) {
1247 c.xyzt.x = x[i] * METERS_in;
1248 c.xyzt.y = y[i] * METERS_in;
1250 c = proj_trans(info_trans->pj, dir, c);
1251 if ((ok = proj_errno(info_trans->pj)) < 0)
1255 x[i] = c.lp.lam * RAD_TO_DEG;
1256 y[i] = c.lp.phi * RAD_TO_DEG;
1265 for (i = 0; i < n; i++) {
1267 c.xyzt.x = x[i] * METERS_in;
1268 c.xyzt.y = y[i] * METERS_in;
1270 c = proj_trans(info_trans->pj, dir, c);
1271 if ((ok = proj_errno(info_trans->pj)) < 0)
1274 x[i] = c.xy.x / METERS_out;
1275 y[i] = c.xy.y / METERS_out;
1283 G_warning(_(
"proj_trans() failed: %s"), proj_errno_string(ok));
1287 const struct pj_info *p_in, *p_out;
1289 if (dir == PJ_FWD) {
1298 METERS_in = p_in->meters;
1299 METERS_out = p_out->meters;
1302 z = G_malloc(
sizeof(
double) * n);
1304 for (i = 0; i < n; ++i)
1308 if (strncmp(p_in->proj,
"ll", 2) == 0) {
1309 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1311 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1316 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1321 if (strncmp(p_out->proj,
"ll", 2) == 0) {
1323 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1328 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1336 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1364 const struct pj_info *info_in,
const struct pj_info *info_out)
1368 struct pj_info info_trans;
1375 METERS_in = info_in->meters;
1376 METERS_out = info_out->meters;
1378 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1380 c.lpzt.lam = (*x) / RAD_TO_DEG;
1381 c.lpzt.phi = (*y) / RAD_TO_DEG;
1384 c = proj_trans(info_trans.pj, PJ_FWD, c);
1385 ok = proj_errno(info_trans.pj);
1387 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1389 *x = c.lp.lam * RAD_TO_DEG;
1390 *y = c.lp.phi * RAD_TO_DEG;
1394 *x = c.xy.x / METERS_out;
1395 *y = c.xy.y / METERS_out;
1400 c.xyzt.x = *x * METERS_in;
1401 c.xyzt.y = *y * METERS_in;
1404 c = proj_trans(info_trans.pj, PJ_FWD, c);
1405 ok = proj_errno(info_trans.pj);
1407 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1409 *x = c.lp.lam * RAD_TO_DEG;
1410 *y = c.lp.phi * RAD_TO_DEG;
1414 *x = c.xy.x / METERS_out;
1415 *y = c.xy.y / METERS_out;
1418 proj_destroy(info_trans.pj);
1421 G_warning(_(
"proj_trans() failed: %d"), ok);
1427 METERS_in = info_in->meters;
1428 METERS_out = info_out->meters;
1430 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1431 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1432 u = (*x) / RAD_TO_DEG;
1433 v = (*y) / RAD_TO_DEG;
1434 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1435 *x = u * RAD_TO_DEG;
1436 *y = v * RAD_TO_DEG;
1439 u = (*x) / RAD_TO_DEG;
1440 v = (*y) / RAD_TO_DEG;
1441 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1442 *x = u / METERS_out;
1443 *y = v / METERS_out;
1447 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1450 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1451 *x = u * RAD_TO_DEG;
1452 *y = v * RAD_TO_DEG;
1457 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1458 *x = u / METERS_out;
1459 *y = v / METERS_out;
1463 G_warning(_(
"pj_transform() failed: %s"), pj_strerrno(ok));
1493 const struct pj_info *info_in,
const struct pj_info *info_out)
1499 struct pj_info info_trans;
1506 METERS_in = info_in->meters;
1507 METERS_out = info_out->meters;
1510 h = G_malloc(
sizeof *h * count);
1512 for (i = 0; i <
count; ++i)
1517 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1519 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1520 for (i = 0; i <
count; i++) {
1522 c.lpzt.lam = x[i] / RAD_TO_DEG;
1523 c.lpzt.phi = y[i] / RAD_TO_DEG;
1525 c = proj_trans(info_trans.pj, PJ_FWD, c);
1526 if ((ok = proj_errno(info_trans.pj)) < 0)
1529 x[i] = c.lp.lam * RAD_TO_DEG;
1530 y[i] = c.lp.phi * RAD_TO_DEG;
1534 for (i = 0; i <
count; i++) {
1536 c.lpzt.lam = x[i] / RAD_TO_DEG;
1537 c.lpzt.phi = y[i] / RAD_TO_DEG;
1539 c = proj_trans(info_trans.pj, PJ_FWD, c);
1540 if ((ok = proj_errno(info_trans.pj)) < 0)
1543 x[i] = c.xy.x / METERS_out;
1544 y[i] = c.xy.y / METERS_out;
1550 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1551 for (i = 0; i <
count; i++) {
1553 c.xyzt.x = x[i] * METERS_in;
1554 c.xyzt.y = y[i] * METERS_in;
1556 c = proj_trans(info_trans.pj, PJ_FWD, c);
1557 if ((ok = proj_errno(info_trans.pj)) < 0)
1560 x[i] = c.lp.lam * RAD_TO_DEG;
1561 y[i] = c.lp.phi * RAD_TO_DEG;
1565 for (i = 0; i <
count; i++) {
1567 c.xyzt.x = x[i] * METERS_in;
1568 c.xyzt.y = y[i] * METERS_in;
1570 c = proj_trans(info_trans.pj, PJ_FWD, c);
1571 if ((ok = proj_errno(info_trans.pj)) < 0)
1574 x[i] = c.xy.x / METERS_out;
1575 y[i] = c.xy.y / METERS_out;
1581 proj_destroy(info_trans.pj);
1584 G_warning(_(
"proj_trans() failed: %d"), ok);
1587 METERS_in = info_in->meters;
1588 METERS_out = info_out->meters;
1591 h = G_malloc(
sizeof *h * count);
1593 for (i = 0; i <
count; ++i)
1597 if (strncmp(info_in->proj,
"ll", 2) == 0) {
1598 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1600 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1605 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1610 if (strncmp(info_out->proj,
"ll", 2) == 0) {
1612 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1617 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1625 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.
void G_get_window(struct Cell_head *window)
Get the current region.
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
struct Key_Value * G_get_projunits(void)
Gets units information for location.
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_unset_window()
Unset current region.
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
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 ...
int G_debug(int level, const char *msg,...)
Print debugging message.
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
char * G_store_upper(const char *s)
Copy string to allocated memory and convert copied string to upper case.
int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys, const struct Key_Value *in_units_keys)
Create a pj_info struct Co-ordinate System definition from a set of PROJ_INFO / PROJ_UNITS-style key-...
#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.