GRASS GIS 8 Programmer's Manual 8.3.0(2023)-exported
do_proj.c
Go to the documentation of this file.
1/**
2 \file do_proj.c
3
4 \brief GProj library - Functions for re-projecting point data
5
6 \author Original Author unknown, probably Soil Conservation Service
7 Eric Miller, Paul Kelly, Markus Metz
8
9 (C) 2003-2008,2018 by the GRASS Development Team
10
11 This program is free software under the GNU General Public
12 License (>=v2). Read the file COPYING that comes with GRASS
13 for details.
14**/
15
16#include <stdio.h>
17#include <string.h>
18#include <ctype.h>
19#include <math.h>
20
21#include <grass/gis.h>
22#include <grass/gprojects.h>
23#include <grass/glocale.h>
24
25/* a couple defines to simplify reading the function */
26#define MULTIPLY_LOOP(x, y, c, m) \
27 do { \
28 for (i = 0; i < c; ++i) { \
29 x[i] *= m; \
30 y[i] *= m; \
31 } \
32 } while (0)
33
34#define DIVIDE_LOOP(x, y, c, m) \
35 do { \
36 for (i = 0; i < c; ++i) { \
37 x[i] /= m; \
38 y[i] /= m; \
39 } \
40 } while (0)
41
42static double METERS_in = 1.0, METERS_out = 1.0;
43
44#ifdef HAVE_PROJ_H
45#if PROJ_VERSION_MAJOR >= 6
46int get_pj_area(const struct pj_info *iproj, double *xmin, double *xmax,
47 double *ymin, double *ymax)
48{
49 struct Cell_head window;
50
51 /* modules must set the current window, do not unset this window here */
52 /* G_unset_window(); */
53 G_get_set_window(&window);
54 *xmin = window.west;
55 *xmax = window.east;
56 *ymin = window.south;
57 *ymax = window.north;
58
59 if (window.proj != PROJECTION_LL) {
60 /* transform to ll equivalent */
61 double estep, nstep;
62 double x[85], y[85];
63 int i;
64 const char *projstr = NULL;
65 char *indef = NULL;
66 struct pj_info oproj, tproj; /* proj parameters */
67
68 oproj.pj = NULL;
69 oproj.proj[0] = '\0';
70 tproj.def = NULL;
71
72 if (proj_get_type(iproj->pj) == PJ_TYPE_BOUND_CRS) {
73 PJ *source_crs;
74
75 source_crs = proj_get_source_crs(NULL, iproj->pj);
76 if (source_crs) {
77 projstr =
78 proj_as_proj_string(NULL, source_crs, PJ_PROJ_5, NULL);
79 if (projstr) {
80 indef = G_store(projstr);
81 }
82 proj_destroy(source_crs);
83 }
84 }
85 else {
86 projstr = proj_as_proj_string(NULL, iproj->pj, PJ_PROJ_5, NULL);
87 if (projstr) {
88 indef = G_store(projstr);
89 }
90 }
91
92 if (indef == NULL)
93 indef = G_store(iproj->def);
94
95 /* needs +over to properly cross the anti-meridian
96 * the +over switch can be used to disable the default wrapping
97 * of output longitudes to the range -180 to 180 */
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);
101
102 if (tproj.pj == NULL) {
103 G_warning(_("proj_create() failed for '%s'"), tproj.def);
104 G_free(indef);
105 G_free(tproj.def);
106 proj_destroy(tproj.pj);
107
108 return 0;
109 }
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);
113 G_free(indef);
114 G_free(tproj.def);
115 proj_destroy(tproj.pj);
116
117 return 0;
118 }
119 else {
120 G_debug(1, "proj_create() projstr '%s'", projstr);
121 }
122 G_free(indef);
123
124 /* inpired by gdal/ogr/ogrct.cpp OGRProjCT::ListCoordinateOperations()
125 */
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);
130 y[i] = window.north;
131
132 x[i + 20] = window.west + estep * (i + 1);
133 y[i + 20] = window.south;
134
135 x[i + 40] = window.west;
136 y[i + 40] = window.south + nstep * (i + 1);
137
138 x[i + 60] = window.east;
139 y[i + 60] = window.south + nstep * (i + 1);
140 }
141 x[80] = window.west;
142 y[80] = window.north;
143 x[81] = window.west;
144 y[81] = window.south;
145 x[82] = window.east;
146 y[82] = window.north;
147 x[83] = window.east;
148 y[83] = window.south;
149 x[84] = (window.west + window.east) / 2.;
150 y[84] = (window.north + window.south) / 2.;
151
152 GPJ_transform_array(iproj, &oproj, &tproj, PJ_FWD, x, y, NULL, 85);
153
154 proj_destroy(tproj.pj);
155 G_free(tproj.def);
156
157 *xmin = *xmax = x[84];
158 *ymin = *ymax = y[84];
159 for (i = 0; i < 84; i++) {
160 if (*xmin > x[i])
161 *xmin = x[i];
162 if (*xmax < x[i])
163 *xmax = x[i];
164 if (*ymin > y[i])
165 *ymin = y[i];
166 if (*ymax < y[i])
167 *ymax = y[i];
168 }
169
170 /* The west longitude is generally lower than the east longitude,
171 * except for areas of interest that go across the anti-meridian.
172 * do not reduce global coverage to a small north-south strip
173 */
174 if (*xmin < -180 && *xmax < 180 && *xmin + 360 > *xmax) {
175 /* must be crossing the anti-meridian at 180W */
176 *xmin += 360;
177 }
178 else if (*xmax > 180 && *xmin > -180 && *xmax - 360 < *xmin) {
179 /* must be crossing the anti-meridian at 180E */
180 *xmax -= 360;
181 }
182
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);
187
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);
192
193 /* test valid values, as in
194 * gdal/ogr/ogrct.cpp
195 * OGRCoordinateTransformationOptions::SetAreaOfInterest()
196 */
197 if (fabs(*xmin) > 180) {
198 G_warning(_("Invalid west longitude %g"), *xmin);
199 return 0;
200 }
201 if (fabs(*xmax) > 180) {
202 G_warning(_("Invalid east longitude %g"), *xmax);
203 return 0;
204 }
205 if (fabs(*ymin) > 90) {
206 G_warning(_("Invalid south latitude %g"), *ymin);
207 return 0;
208 }
209 if (fabs(*ymax) > 90) {
210 G_warning(_("Invalid north latitude %g"), *ymax);
211 return 0;
212 }
213 if (*ymin > *ymax) {
214 G_warning(_("South %g is larger than north %g"), *ymin, *ymax);
215 return 0;
216 }
217 }
218 G_debug(1, "get_pj_area(): xmin %g, xmax %g, ymin %g, ymax %g", *xmin,
219 *xmax, *ymin, *ymax);
220
221 return 1;
222}
223
224char *get_pj_type_string(PJ *pj)
225{
226 char *pj_type = NULL;
227
228 switch (proj_get_type(pj)) {
229 case PJ_TYPE_UNKNOWN:
230 G_asprintf(&pj_type, "unknown");
231 break;
232 case PJ_TYPE_ELLIPSOID:
233 G_asprintf(&pj_type, "ellipsoid");
234 break;
235 case PJ_TYPE_PRIME_MERIDIAN:
236 G_asprintf(&pj_type, "prime meridian");
237 break;
238 case PJ_TYPE_GEODETIC_REFERENCE_FRAME:
239 G_asprintf(&pj_type, "geodetic reference frame");
240 break;
241 case PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME:
242 G_asprintf(&pj_type, "dynamic geodetic reference frame");
243 break;
244 case PJ_TYPE_VERTICAL_REFERENCE_FRAME:
245 G_asprintf(&pj_type, "vertical reference frame");
246 break;
247 case PJ_TYPE_DYNAMIC_VERTICAL_REFERENCE_FRAME:
248 G_asprintf(&pj_type, "dynamic vertical reference frame");
249 break;
250 case PJ_TYPE_DATUM_ENSEMBLE:
251 G_asprintf(&pj_type, "datum ensemble");
252 break;
253
254 /** Abstract type, not returned by proj_get_type() */
255 case PJ_TYPE_CRS:
256 G_asprintf(&pj_type, "crs");
257 break;
258 case PJ_TYPE_GEODETIC_CRS:
259 G_asprintf(&pj_type, "geodetic crs");
260 break;
261 case PJ_TYPE_GEOCENTRIC_CRS:
262 G_asprintf(&pj_type, "geocentric crs");
263 break;
264
265 /** proj_get_type() will never return that type, but
266 * PJ_TYPE_GEOGRAPHIC_2D_CRS or PJ_TYPE_GEOGRAPHIC_3D_CRS. */
267 case PJ_TYPE_GEOGRAPHIC_CRS:
268 G_asprintf(&pj_type, "geographic crs");
269 break;
270 case PJ_TYPE_GEOGRAPHIC_2D_CRS:
271 G_asprintf(&pj_type, "geographic 2D crs");
272 break;
273 case PJ_TYPE_GEOGRAPHIC_3D_CRS:
274 G_asprintf(&pj_type, "geographic 3D crs");
275 break;
276 case PJ_TYPE_VERTICAL_CRS:
277 G_asprintf(&pj_type, "vertical crs");
278 break;
279 case PJ_TYPE_PROJECTED_CRS:
280 G_asprintf(&pj_type, "projected crs");
281 break;
282 case PJ_TYPE_COMPOUND_CRS:
283 G_asprintf(&pj_type, "compound crs");
284 break;
285 case PJ_TYPE_TEMPORAL_CRS:
286 G_asprintf(&pj_type, "temporal crs");
287 break;
288 case PJ_TYPE_ENGINEERING_CRS:
289 G_asprintf(&pj_type, "engineering crs");
290 break;
291 case PJ_TYPE_BOUND_CRS:
292 G_asprintf(&pj_type, "bound crs");
293 break;
294 case PJ_TYPE_OTHER_CRS:
295 G_asprintf(&pj_type, "other crs");
296 break;
297 case PJ_TYPE_CONVERSION:
298 G_asprintf(&pj_type, "conversion");
299 break;
300 case PJ_TYPE_TRANSFORMATION:
301 G_asprintf(&pj_type, "transformation");
302 break;
303 case PJ_TYPE_CONCATENATED_OPERATION:
304 G_asprintf(&pj_type, "concatenated operation");
305 break;
306 case PJ_TYPE_OTHER_COORDINATE_OPERATION:
307 G_asprintf(&pj_type, "other coordinate operation");
308 break;
309 default:
310 G_asprintf(&pj_type, "unknown");
311 break;
312 }
313
314 return pj_type;
315}
316
317PJ *get_pj_object(const struct pj_info *in_gpj, char **in_defstr)
318{
319 PJ *in_pj = NULL;
320
321 *in_defstr = NULL;
322
323 /* 1. SRID, 2. WKT, 3. standard pj from pj_get_kv */
324 if (in_gpj->srid) {
325 G_debug(1, "Trying SRID '%s' ...", in_gpj->srid);
326 in_pj = proj_create(PJ_DEFAULT_CTX, in_gpj->srid);
327 if (!in_pj) {
328 G_warning(_("Unrecognized SRID '%s'"), in_gpj->srid);
329 }
330 else {
331 *in_defstr = G_store(in_gpj->srid);
332 /* PROJ will do the unit conversion if set up from srid
333 * -> disable unit conversion for GPJ_transform */
334 /* ugly hack */
335 ((struct pj_info *)in_gpj)->meters = 1;
336 }
337 }
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);
341 if (!in_pj) {
342 G_warning(_("Unrecognized WKT '%s'"), in_gpj->wkt);
343 }
344 else {
345 *in_defstr = G_store(in_gpj->wkt);
346 /* PROJ will do the unit conversion if set up from wkt
347 * -> disable unit conversion for GPJ_transform */
348 /* ugly hack */
349 ((struct pj_info *)in_gpj)->meters = 1;
350 }
351 }
352 if (!in_pj && in_gpj->pj) {
353 in_pj = proj_clone(PJ_DEFAULT_CTX, in_gpj->pj);
354 *in_defstr = G_store(proj_as_wkt(NULL, in_pj, PJ_WKT2_LATEST, NULL));
355 if (*in_defstr && !**in_defstr)
356 *in_defstr = NULL;
357 }
358
359 if (!in_pj) {
360 G_warning(_("Unable to create PROJ object"));
361
362 return NULL;
363 }
364
365 /* Even Rouault:
366 * if info_in->def contains a +towgs84/+nadgrids clause,
367 * this pipeline would apply it, whereas you probably only want
368 * the reverse projection, and no datum shift.
369 * The easiest would probably to mess up with the PROJ string.
370 * Otherwise with the PROJ API, you could
371 * instantiate a PJ object from the string,
372 * check if it is a BoundCRS with proj_get_source_crs(),
373 * and in that case, take the source CRS with proj_get_source_crs(),
374 * and do the inverse transform on it */
375
376 if (proj_get_type(in_pj) == PJ_TYPE_BOUND_CRS) {
377 PJ *source_crs;
378
379 G_debug(1, "found bound crs");
380 source_crs = proj_get_source_crs(NULL, in_pj);
381 if (source_crs) {
382 *in_defstr =
383 G_store(proj_as_wkt(NULL, source_crs, PJ_WKT2_LATEST, NULL));
384 if (*in_defstr && !**in_defstr)
385 *in_defstr = NULL;
386 in_pj = source_crs;
387 }
388 }
389
390 return in_pj;
391}
392#endif
393#endif
394
395/**
396 * \brief Create a PROJ transformation object to transform coordinates
397 * from an input SRS to an output SRS
398 *
399 * After the transformation has been initialized with this function,
400 * coordinates can be transformed from input SRS to output SRS with
401 * GPJ_transform() and direction = PJ_FWD, and back from output SRS to
402 * input SRS with direction = OJ_INV.
403 * If coordinates should be transformed between the input SRS and its
404 * latlong equivalent, an uninitialized info_out with
405 * info_out->pj = NULL can be passed to the function. In this case,
406 * coordinates will be transformed between the input SRS and its
407 * latlong equivalent, and for PROJ 5+, the transformation object is
408 * created accordingly, while for PROJ 4, the output SRS is created as
409 * latlong equivalent of the input SRS
410 *
411 * PROJ 5+:
412 * info_in->pj must not be null
413 * if info_out->pj is null, assume info_out to be the ll equivalent
414 * of info_in
415 * create info_trans as conversion from info_in to its ll equivalent
416 * NOTE: this is the inverse of the logic of PROJ 5 which by default
417 * converts from ll to a given SRS, not from a given SRS to ll
418 * thus PROJ 5+ itself uses an inverse transformation in the
419 * first step of the pipeline for proj_create_crs_to_crs()
420 * if info_trans->def is not NULL, this pipeline definition will be
421 * used to create a transformation object
422 * PROJ 4:
423 * info_in->pj must not be null
424 * if info_out->pj is null, create info_out as ll equivalent
425 * else do nothing, info_trans is not used
426 *
427 * \param info_in pointer to pj_info struct for input co-ordinate system
428 * \param info_out pointer to pj_info struct for output co-ordinate system
429 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
430 *5+)
431 *
432 * \return 1 on success, -1 on failure
433 **/
434int GPJ_init_transform(const struct pj_info *info_in,
435 const struct pj_info *info_out,
436 struct pj_info *info_trans)
437{
438 if (info_in->pj == NULL)
439 G_fatal_error(_("Input coordinate system is NULL"));
440
441 if (info_in->def == NULL)
442 G_fatal_error(_("Input coordinate system definition is NULL"));
443
444#ifdef HAVE_PROJ_H
445#if PROJ_VERSION_MAJOR >= 6
446
447 /* PROJ6+: enforce axis order easting, northing
448 * +axis=enu (works with proj-4.8+) */
449
450 info_trans->pj = NULL;
451 info_trans->meters = 1.;
452 info_trans->zone = 0;
453 sprintf(info_trans->proj, "pipeline");
454
455 /* user-provided pipeline */
456 if (info_trans->def) {
457 const char *projstr;
458
459 /* info_in->pj, info_in->proj, info_out->pj, info_out->proj
460 * must be set */
461 if (!info_in->pj || !info_in->proj[0] || !info_out->pj ||
462 !info_out->proj[0]) {
463 G_warning(_(
464 "A custom pipeline requires input and output projection info"));
465
466 return -1;
467 }
468
469 /* create a pj from user-defined transformation pipeline */
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);
473
474 return -1;
475 }
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);
479
480 return -1;
481 }
482 else {
483 /* make sure axis order is easting, northing
484 * proj_normalize_for_visualization() does not work here
485 * because source and target CRS are unknown to PROJ
486 * remove any "+step +proj=axisswap +order=2,1" ?
487 * */
488 info_trans->def = G_store(projstr);
489
490 if (strstr(info_trans->def, "axisswap")) {
491 G_warning(
492 _("The transformation pipeline contains an '%s' step. "
493 "Remove this step if easting and northing are swapped in "
494 "the output."),
495 "axisswap");
496 }
497
498 G_debug(1, "proj_create() pipeline: %s", info_trans->def);
499
500 /* the user-provided PROJ pipeline is supposed to do
501 * all the needed unit conversions */
502 /* ugly hack */
503 ((struct pj_info *)info_in)->meters = 1;
504 ((struct pj_info *)info_out)->meters = 1;
505 }
506 }
507 /* if no output CRS is defined,
508 * assume info_out to be ll equivalent of info_in */
509 else if (info_out->pj == NULL) {
510 const char *projstr = NULL;
511 char *indef = NULL;
512
513 /* get PROJ-style definition */
514 indef = G_store(info_in->def);
515 G_debug(1, "ll equivalent definition: %s", indef);
516
517 /* what about axis order?
518 * is it always enu?
519 * probably yes, as long as there is no +proj=axisswap step */
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);
524 G_free(indef);
525
526 return -1;
527 }
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);
531 G_free(indef);
532
533 return -1;
534 }
535 G_free(indef);
536 }
537 /* input and output CRS are available */
538 else if (info_in->def && info_out->pj && info_out->def) {
539 char *indef = NULL, *outdef = NULL;
540 char *insrid = NULL, *outsrid = NULL;
541 PJ *in_pj, *out_pj;
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;
547
548 /* get pj_area */
549 /* do it here because get_pj_area() will use
550 * the PROJ definition for simple transformation to the
551 * ll equivalent and we need to do unit conversion */
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);
555 }
556 else {
557 G_warning(_("Unable to determine area of interest for '%s'"),
558 info_in->def);
559
560 return -1;
561 }
562
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));
565
566 /* PROJ6+: EPSG must be uppercase EPSG */
567 if (info_in->srid) {
568 if (strncmp(info_in->srid, "epsg", 4) == 0) {
569 insrid = G_store_upper(info_in->srid);
570 G_free(info_in->srid);
571 ((struct pj_info *)info_in)->srid = insrid;
572 insrid = NULL;
573 }
574 }
575
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);
579
580 return -1;
581 }
582 G_debug(1, "Input CRS definition: %s", indef);
583
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));
586
587 /* PROJ6+: EPSG must be uppercase EPSG */
588 if (info_out->srid) {
589 if (strncmp(info_out->srid, "epsg", 4) == 0) {
590 outsrid = G_store_upper(info_out->srid);
591 G_free(info_out->srid);
592 ((struct pj_info *)info_out)->srid = outsrid;
593 outsrid = NULL;
594 }
595 }
596
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);
600
601 return -1;
602 }
603 G_debug(1, "Output CRS definition: %s", outdef);
604
605 /* check number of operations */
606
607 operation_ctx =
608 proj_create_operation_factory_context(PJ_DEFAULT_CTX, NULL);
609 /* proj_create_operations() works only if both source_crs
610 * and target_crs are found in the proj db
611 * if any is not found, proj can not get a list of operations
612 * and we have to take care of datumshift manually */
613 /* list all operations irrespecitve of area and
614 * grid availability */
615 op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
616 operation_ctx);
617 proj_operation_factory_context_destroy(operation_ctx);
618
619 op_count = 0;
620 if (op_list)
621 op_count = proj_list_get_count(op_list);
622 if (op_count > 1) {
623 int i;
624
625 G_important_message(_("Found %d possible transformations"),
626 op_count);
627 for (i = 0; i < op_count; i++) {
628 const char *area_of_use, *projstr;
629 double e, w, s, n;
630 PJ_PROJ_INFO pj_info;
631 PJ *op, *op_norm;
632
633 op = proj_list_get(PJ_DEFAULT_CTX, op_list, i);
634 op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
635
636 if (!op_norm) {
637 G_warning(_("proj_normalize_for_visualization() failed for "
638 "operation %d"),
639 i + 1);
640 }
641 else {
642 proj_destroy(op);
643 op = op_norm;
644 }
645
646 pj_info = proj_pj_info(op);
647 proj_get_area_of_use(NULL, op, &w, &s, &e, &n, &area_of_use);
648 G_important_message("************************");
649 G_important_message(_("Operation %d:"), i + 1);
650 if (pj_info.description) {
651 G_important_message(_("Description: %s"),
652 pj_info.description);
653 }
654 if (area_of_use) {
656 G_important_message(_("Area of use: %s"), area_of_use);
657 }
658 if (pj_info.accuracy > 0) {
660 G_important_message(_("Accuracy within area of use: %g m"),
661 pj_info.accuracy);
662 }
663#if PROJ_VERSION_NUM >= 6020000
664 const char *str = proj_get_remarks(op);
665
666 if (str && *str) {
668 G_important_message(_("Remarks: %s"), str);
669 }
670 str = proj_get_scope(op);
671 if (str && *str) {
673 G_important_message(_("Scope: %s"), str);
674 }
675#endif
676
677 projstr = proj_as_proj_string(NULL, op, PJ_PROJ_5, NULL);
678 if (projstr) {
680 G_important_message(_("PROJ string:"));
681 G_important_message("%s", projstr);
682 }
683 proj_destroy(op);
684 }
685 G_important_message("************************");
686
687 G_important_message(_("See also output of:"));
688 G_important_message("projinfo -o PROJ -s \"%s\" -t \"%s\"", indef,
689 outdef);
690 G_important_message(_("Please provide the appropriate PROJ string "
691 "with the %s option"),
692 "pipeline");
693 G_important_message("************************");
694 }
695
696 if (op_list)
697 proj_list_destroy(op_list);
698
699 /* following code copied from proj_create_crs_to_crs_from_pj()
700 * in proj src/4D_api.cpp
701 * but using PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT */
702
703 /* now use the current region as area of interest */
704 operation_ctx =
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);
714 /* The operations are sorted with the most relevant ones first:
715 * by descending area (intersection of the transformation area
716 * with the area of interest, or intersection of the
717 * transformation with the area of use of the CRS),
718 * and by increasing accuracy.
719 * Operations with unknown accuracy are sorted last,
720 * whatever their area.
721 */
722 op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
723 operation_ctx);
724 proj_operation_factory_context_destroy(operation_ctx);
725 op_count_area = 0;
726 if (op_list)
727 op_count_area = proj_list_get_count(op_list);
728 if (op_count_area == 0) {
729 /* no operations */
730 info_trans->pj = NULL;
731 }
732 else if (op_count_area == 1) {
733 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
734 }
735 else { /* op_count_area > 1 */
736 /* can't use pj_create_prepared_operations()
737 * this is a PROJ-internal function
738 * trust the sorting of PROJ and use the first one */
739 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
740 }
741 if (op_list)
742 proj_list_destroy(op_list);
743
744 /* try proj_create_crs_to_crs() */
745 /*
746 G_debug(1, "trying %s to %s", indef, outdef);
747 */
748
749 /* proj_create_crs_to_crs() does not work because it calls
750 * proj_create_crs_to_crs_from_pj() which calls
751 * proj_operation_factory_context_set_spatial_criterion()
752 * with PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION
753 * instead of
754 * PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT
755 *
756 * fixed in PROJ master, probably available with PROJ 7.3.x */
757
758 /*
759 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
760 indef,
761 outdef,
762 pj_area);
763 */
764
765 if (in_pj)
766 proj_destroy(in_pj);
767 if (out_pj)
768 proj_destroy(out_pj);
769
770 if (info_trans->pj) {
771 const char *projstr;
772 PJ *pj_norm = NULL;
773
774 G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ%d",
775 PROJ_VERSION_MAJOR);
776
777 projstr =
778 proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);
779
780 info_trans->def = G_store(projstr);
781
782 if (projstr) {
783 /* make sure axis order is easting, northing
784 * proj_normalize_for_visualization() requires
785 * source and target CRS
786 * -> does not work with ll equivalent of input:
787 * no target CRS in +proj=pipeline +step +inv %s */
788 pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX,
789 info_trans->pj);
790
791 if (!pj_norm) {
792 G_warning(
793 _("proj_normalize_for_visualization() failed for '%s'"),
794 info_trans->def);
795 }
796 else {
797 projstr =
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);
803 }
804 else {
805 proj_destroy(pj_norm);
806 G_warning(_("No PROJ definition for normalized version "
807 "of '%s'"),
808 info_trans->def);
809 }
810 }
811 G_important_message(_("Selected PROJ pipeline:"));
812 G_important_message(_("%s"), info_trans->def);
813 G_important_message("************************");
814 }
815 else {
816 proj_destroy(info_trans->pj);
817 info_trans->pj = NULL;
818 }
819 }
820
821 if (pj_area)
822 proj_area_destroy(pj_area);
823
824 if (insrid)
825 G_free(insrid);
826 if (outsrid)
827 G_free(outsrid);
828 G_free(indef);
829 G_free(outdef);
830 }
831 if (info_trans->pj == NULL) {
832 G_warning(_("proj_create() failed for '%s'"), info_trans->def);
833
834 return -1;
835 }
836#else /* PROJ 5 */
837 info_trans->pj = NULL;
838 info_trans->meters = 1.;
839 info_trans->zone = 0;
840 sprintf(info_trans->proj, "pipeline");
841
842 /* user-provided pipeline */
843 if (info_trans->def) {
844 /* create a pj from user-defined transformation pipeline */
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);
848
849 return -1;
850 }
851 }
852 /* if no output CRS is defined,
853 * assume info_out to be ll equivalent of info_in */
854 else if (info_out->pj == NULL) {
855 G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",
856 info_in->def);
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);
860
861 return -1;
862 }
863 }
864 else if (info_in->def && info_out->pj && info_out->def) {
865 char *indef = NULL, *outdef = NULL;
866 char *insrid = NULL, *outsrid = NULL;
867
868 /* PROJ5: EPSG must be lowercase epsg */
869 if (info_in->srid) {
870 if (strncmp(info_in->srid, "EPSG", 4) == 0)
871 insrid = G_store_lower(info_in->srid);
872 else
873 insrid = G_store(info_in->srid);
874 }
875
876 if (info_out->pj && info_out->srid) {
877 if (strncmp(info_out->srid, "EPSG", 4) == 0)
878 outsrid = G_store_lower(info_out->srid);
879 else
880 outsrid = G_store(info_out->srid);
881 }
882
883 info_trans->pj = NULL;
884
885 if (insrid && outsrid) {
886 G_asprintf(&indef, "%s", insrid);
887 G_asprintf(&outdef, "%s", outsrid);
888
889 G_asprintf(&(info_trans->def),
890 "+proj=pipeline +step +inv +init=%s +step +init=%s",
891 indef, outdef);
892
893 /* try proj_create_crs_to_crs() */
894 info_trans->pj =
895 proj_create_crs_to_crs(PJ_DEFAULT_CTX, indef, outdef, NULL);
896 }
897
898 if (info_trans->pj) {
899 G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ5");
900 }
901 else {
902 if (indef) {
903 G_free(indef);
904 indef = NULL;
905 }
906 if (insrid) {
907 G_asprintf(&indef, "+init=%s", insrid);
908 }
909 else {
910 G_asprintf(&indef, "%s", info_in->def);
911 }
912
913 if (outdef) {
914 G_free(outdef);
915 outdef = NULL;
916 }
917 if (outsrid) {
918 G_asprintf(&outdef, "+init=%s", outsrid);
919 }
920 else {
921 G_asprintf(&outdef, "%s", info_out->def);
922 }
923
924 /* try proj_create() with +proj=pipeline +step +inv %s +step %s" */
925 G_asprintf(&(info_trans->def),
926 "+proj=pipeline +step +inv %s +step %s", indef, outdef);
927
928 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
929 }
930 if (insrid)
931 G_free(insrid);
932 if (outsrid)
933 G_free(outsrid);
934 G_free(indef);
935 G_free(outdef);
936 }
937 if (info_trans->pj == NULL) {
938 G_warning(_("proj_create() failed for '%s'"), info_trans->def);
939
940 return -1;
941 }
942
943#endif
944#else /* PROJ 4 */
945 if (info_out->pj == NULL) {
946 if (GPJ_get_equivalent_latlong(info_out, info_in) < 0) {
947 G_warning(_("Unable to create latlong equivalent for '%s'"),
948 info_in->def);
949
950 return -1;
951 }
952 }
953#endif
954
955 return 1;
956}
957
958/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
959
960/**
961 * \brief Re-project a point between two co-ordinate systems using a
962 * transformation object prepared with GPJ_prepare_pj()
963 *
964 * This function takes pointers to three pj_info structures as arguments,
965 * and projects a point between the input and output co-ordinate system.
966 * The pj_info structure info_trans must have been initialized with
967 * GPJ_init_transform().
968 * The direction determines if a point is projected from input CRS to
969 * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
970 * The easting, northing, and height of the point are contained in the
971 * pointers passed to the function; these will be overwritten by the
972 * coordinates of the transformed point.
973 *
974 * \param info_in pointer to pj_info struct for input co-ordinate system
975 * \param info_out pointer to pj_info struct for output co-ordinate system
976 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
977 *5+) \param dir direction of the transformation (PJ_FWD or PJ_INV) \param x
978 *Pointer to a double containing easting or longitude \param y Pointer to a
979 *double containing northing or latitude \param z Pointer to a double containing
980 *height, or NULL
981 *
982 * \return Return value from PROJ proj_trans() function
983 **/
984
985int 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)
988{
989 int ok = 0;
990
991#ifdef HAVE_PROJ_H
992 /* PROJ 5+ variant */
993 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
994 PJ_COORD c;
995
996 if (info_in->pj == NULL)
997 G_fatal_error(_("No input projection"));
998
999 if (info_trans->pj == NULL)
1000 G_fatal_error(_("No transformation object"));
1001
1002 in_deg2rad = out_rad2deg = 1;
1003 if (dir == PJ_FWD) {
1004 /* info_in -> info_out */
1005 METERS_in = info_in->meters;
1006 in_is_ll = !strncmp(info_in->proj, "ll", 2);
1007#if PROJ_VERSION_MAJOR >= 6
1008 /* PROJ 6+: conversion to radians is not always needed:
1009 * if proj_angular_input(info_trans->pj, dir) == 1
1010 * -> convert from degrees to radians */
1011 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1012 in_deg2rad = 0;
1013 }
1014#endif
1015 if (info_out->pj) {
1016 METERS_out = info_out->meters;
1017 out_is_ll = !strncmp(info_out->proj, "ll", 2);
1018#if PROJ_VERSION_MAJOR >= 6
1019 /* PROJ 6+: conversion to radians is not always needed:
1020 * if proj_angular_input(info_trans->pj, dir) == 1
1021 * -> convert from degrees to radians */
1022 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1023 out_rad2deg = 0;
1024 }
1025#endif
1026 }
1027 else {
1028 METERS_out = 1.0;
1029 out_is_ll = 1;
1030 }
1031 }
1032 else {
1033 /* info_out -> info_in */
1034 METERS_out = info_in->meters;
1035 out_is_ll = !strncmp(info_in->proj, "ll", 2);
1036#if PROJ_VERSION_MAJOR >= 6
1037 /* PROJ 6+: conversion to radians is not always needed:
1038 * if proj_angular_input(info_trans->pj, dir) == 1
1039 * -> convert from degrees to radians */
1040 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1041 out_rad2deg = 0;
1042 }
1043#endif
1044 if (info_out->pj) {
1045 METERS_in = info_out->meters;
1046 in_is_ll = !strncmp(info_out->proj, "ll", 2);
1047#if PROJ_VERSION_MAJOR >= 6
1048 /* PROJ 6+: conversion to radians is not always needed:
1049 * if proj_angular_input(info_trans->pj, dir) == 1
1050 * -> convert from degrees to radians */
1051 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1052 in_deg2rad = 0;
1053 }
1054#endif
1055 }
1056 else {
1057 METERS_in = 1.0;
1058 in_is_ll = 1;
1059 }
1060 }
1061
1062 /* prepare */
1063 if (in_is_ll) {
1064 if (in_deg2rad) {
1065 /* convert degrees to radians */
1066 c.lpzt.lam = (*x) / RAD_TO_DEG;
1067 c.lpzt.phi = (*y) / RAD_TO_DEG;
1068 }
1069 else {
1070 c.lpzt.lam = (*x);
1071 c.lpzt.phi = (*y);
1072 }
1073 c.lpzt.z = 0;
1074 if (z)
1075 c.lpzt.z = *z;
1076 c.lpzt.t = 0;
1077 }
1078 else {
1079 /* convert to meters */
1080 c.xyzt.x = *x * METERS_in;
1081 c.xyzt.y = *y * METERS_in;
1082 c.xyzt.z = 0;
1083 if (z)
1084 c.xyzt.z = *z;
1085 c.xyzt.t = 0;
1086 }
1087
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);
1091
1092 /* transform */
1093 c = proj_trans(info_trans->pj, dir, c);
1094 ok = proj_errno(info_trans->pj);
1095
1096 if (ok < 0) {
1097 G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
1098 return ok;
1099 }
1100
1101 /* output */
1102 if (out_is_ll) {
1103 /* convert to degrees */
1104 if (out_rad2deg) {
1105 /* convert radians to degrees */
1106 *x = c.lp.lam * RAD_TO_DEG;
1107 *y = c.lp.phi * RAD_TO_DEG;
1108 }
1109 else {
1110 *x = c.lp.lam;
1111 *y = c.lp.phi;
1112 }
1113 if (z)
1114 *z = c.lpzt.z;
1115 }
1116 else {
1117 /* convert to map units */
1118 *x = c.xyzt.x / METERS_out;
1119 *y = c.xyzt.y / METERS_out;
1120 if (z)
1121 *z = c.xyzt.z;
1122 }
1123#else
1124 /* PROJ 4 variant */
1125 double u, v;
1126 double h = 0.0;
1127 const struct pj_info *p_in, *p_out;
1128
1129 if (info_out == NULL)
1130 G_fatal_error(_("No output projection"));
1131
1132 if (dir == PJ_FWD) {
1133 p_in = info_in;
1134 p_out = info_out;
1135 }
1136 else {
1137 p_in = info_out;
1138 p_out = info_in;
1139 }
1140
1141 METERS_in = p_in->meters;
1142 METERS_out = p_out->meters;
1143
1144 if (z)
1145 h = *z;
1146
1147 if (strncmp(p_in->proj, "ll", 2) == 0) {
1148 u = (*x) / RAD_TO_DEG;
1149 v = (*y) / RAD_TO_DEG;
1150 }
1151 else {
1152 u = *x * METERS_in;
1153 v = *y * METERS_in;
1154 }
1155
1156 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1157
1158 if (ok < 0) {
1159 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1160 return ok;
1161 }
1162
1163 if (strncmp(p_out->proj, "ll", 2) == 0) {
1164 *x = u * RAD_TO_DEG;
1165 *y = v * RAD_TO_DEG;
1166 }
1167 else {
1168 *x = u / METERS_out;
1169 *y = v / METERS_out;
1170 }
1171 if (z)
1172 *z = h;
1173#endif
1174
1175 return ok;
1176}
1177
1178/**
1179 * \brief Re-project an array of points between two co-ordinate systems
1180 * using a transformation object prepared with GPJ_prepare_pj()
1181 *
1182 * This function takes pointers to three pj_info structures as arguments,
1183 * and projects an array of pointd between the input and output
1184 * co-ordinate system. The pj_info structure info_trans must have been
1185 * initialized with GPJ_init_transform().
1186 * The direction determines if a point is projected from input CRS to
1187 * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
1188 * The easting, northing, and height of the point are contained in the
1189 * pointers passed to the function; these will be overwritten by the
1190 * coordinates of the transformed point.
1191 *
1192 * \param info_in pointer to pj_info struct for input co-ordinate system
1193 * \param info_out pointer to pj_info struct for output co-ordinate system
1194 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
1195 *5+) \param dir direction of the transformation (PJ_FWD or PJ_INV) \param x
1196 *pointer to an array of type double containing easting or longitude \param y
1197 *pointer to an array of type double containing northing or latitude \param z
1198 *pointer to an array of type double containing height, or NULL \param n number
1199 *of points in the arrays to be transformed
1200 *
1201 * \return Return value from PROJ proj_trans() function
1202 **/
1203
1204int GPJ_transform_array(const struct pj_info *info_in,
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)
1208{
1209 int ok;
1210 int i;
1211 int has_z = 1;
1212
1213#ifdef HAVE_PROJ_H
1214 /* PROJ 5+ variant */
1215 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1216 PJ_COORD c;
1217
1218 if (info_trans->pj == NULL)
1219 G_fatal_error(_("No transformation object"));
1220
1221 in_deg2rad = out_rad2deg = 1;
1222 if (dir == PJ_FWD) {
1223 /* info_in -> info_out */
1224 METERS_in = info_in->meters;
1225 in_is_ll = !strncmp(info_in->proj, "ll", 2);
1226#if PROJ_VERSION_MAJOR >= 6
1227 /* PROJ 6+: conversion to radians is not always needed:
1228 * if proj_angular_input(info_trans->pj, dir) == 1
1229 * -> convert from degrees to radians */
1230 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1231 in_deg2rad = 0;
1232 }
1233#endif
1234 if (info_out->pj) {
1235 METERS_out = info_out->meters;
1236 out_is_ll = !strncmp(info_out->proj, "ll", 2);
1237#if PROJ_VERSION_MAJOR >= 6
1238 /* PROJ 6+: conversion to radians is not always needed:
1239 * if proj_angular_input(info_trans->pj, dir) == 1
1240 * -> convert from degrees to radians */
1241 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1242 out_rad2deg = 0;
1243 }
1244#endif
1245 }
1246 else {
1247 METERS_out = 1.0;
1248 out_is_ll = 1;
1249 }
1250 }
1251 else {
1252 /* info_out -> info_in */
1253 METERS_out = info_in->meters;
1254 out_is_ll = !strncmp(info_in->proj, "ll", 2);
1255#if PROJ_VERSION_MAJOR >= 6
1256 /* PROJ 6+: conversion to radians is not always needed:
1257 * if proj_angular_input(info_trans->pj, dir) == 1
1258 * -> convert from degrees to radians */
1259 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1260 out_rad2deg = 0;
1261 }
1262#endif
1263 if (info_out->pj) {
1264 METERS_in = info_out->meters;
1265 in_is_ll = !strncmp(info_out->proj, "ll", 2);
1266#if PROJ_VERSION_MAJOR >= 6
1267 /* PROJ 6+: conversion to degrees is not always needed:
1268 * if proj_angular_output(info_trans->pj, dir) == 1
1269 * -> convert from degrees to radians */
1270 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1271 in_deg2rad = 0;
1272 }
1273#endif
1274 }
1275 else {
1276 METERS_in = 1.0;
1277 in_is_ll = 1;
1278 }
1279 }
1280
1281 if (z == NULL) {
1282 z = G_malloc(sizeof(double) * n);
1283 /* they say memset is only guaranteed for chars ;-( */
1284 for (i = 0; i < n; i++)
1285 z[i] = 0.0;
1286 has_z = 0;
1287 }
1288 ok = 0;
1289 if (in_is_ll) {
1290 c.lpzt.t = 0;
1291 if (out_is_ll) {
1292 /* what is more costly ?
1293 * calling proj_trans for each point
1294 * or having three loops over all points ?
1295 * proj_trans_array() itself calls proj_trans() in a loop
1296 * -> one loop over all points is better than
1297 * three loops over all points
1298 */
1299 for (i = 0; i < n; i++) {
1300 if (in_deg2rad) {
1301 /* convert degrees to radians */
1302 c.lpzt.lam = (*x) / RAD_TO_DEG;
1303 c.lpzt.phi = (*y) / RAD_TO_DEG;
1304 }
1305 else {
1306 c.lpzt.lam = (*x);
1307 c.lpzt.phi = (*y);
1308 }
1309 c.lpzt.z = z[i];
1310 c = proj_trans(info_trans->pj, dir, c);
1311 if ((ok = proj_errno(info_trans->pj)) < 0)
1312 break;
1313 if (out_rad2deg) {
1314 /* convert radians to degrees */
1315 x[i] = c.lp.lam * RAD_TO_DEG;
1316 y[i] = c.lp.phi * RAD_TO_DEG;
1317 }
1318 else {
1319 x[i] = c.lp.lam;
1320 y[i] = c.lp.phi;
1321 }
1322 }
1323 }
1324 else {
1325 for (i = 0; i < n; i++) {
1326 if (in_deg2rad) {
1327 /* convert degrees to radians */
1328 c.lpzt.lam = (*x) / RAD_TO_DEG;
1329 c.lpzt.phi = (*y) / RAD_TO_DEG;
1330 }
1331 else {
1332 c.lpzt.lam = (*x);
1333 c.lpzt.phi = (*y);
1334 }
1335 c.lpzt.z = z[i];
1336 c = proj_trans(info_trans->pj, dir, c);
1337 if ((ok = proj_errno(info_trans->pj)) < 0)
1338 break;
1339
1340 /* convert to map units */
1341 x[i] = c.xy.x / METERS_out;
1342 y[i] = c.xy.y / METERS_out;
1343 }
1344 }
1345 }
1346 else {
1347 c.xyzt.t = 0;
1348 if (out_is_ll) {
1349 for (i = 0; i < n; i++) {
1350 /* convert to meters */
1351 c.xyzt.x = x[i] * METERS_in;
1352 c.xyzt.y = y[i] * METERS_in;
1353 c.xyzt.z = z[i];
1354 c = proj_trans(info_trans->pj, dir, c);
1355 if ((ok = proj_errno(info_trans->pj)) < 0)
1356 break;
1357 if (out_rad2deg) {
1358 /* convert radians to degrees */
1359 x[i] = c.lp.lam * RAD_TO_DEG;
1360 y[i] = c.lp.phi * RAD_TO_DEG;
1361 }
1362 else {
1363 x[i] = c.lp.lam;
1364 y[i] = c.lp.phi;
1365 }
1366 }
1367 }
1368 else {
1369 for (i = 0; i < n; i++) {
1370 /* convert to meters */
1371 c.xyzt.x = x[i] * METERS_in;
1372 c.xyzt.y = y[i] * METERS_in;
1373 c.xyzt.z = z[i];
1374 c = proj_trans(info_trans->pj, dir, c);
1375 if ((ok = proj_errno(info_trans->pj)) < 0)
1376 break;
1377 /* convert to map units */
1378 x[i] = c.xy.x / METERS_out;
1379 y[i] = c.xy.y / METERS_out;
1380 }
1381 }
1382 }
1383 if (!has_z)
1384 G_free(z);
1385
1386 if (ok < 0) {
1387 G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
1388 }
1389#else
1390 /* PROJ 4 variant */
1391 const struct pj_info *p_in, *p_out;
1392
1393 if (dir == PJ_FWD) {
1394 p_in = info_in;
1395 p_out = info_out;
1396 }
1397 else {
1398 p_in = info_out;
1399 p_out = info_in;
1400 }
1401
1402 METERS_in = p_in->meters;
1403 METERS_out = p_out->meters;
1404
1405 if (z == NULL) {
1406 z = G_malloc(sizeof(double) * n);
1407 /* they say memset is only guaranteed for chars ;-( */
1408 for (i = 0; i < n; ++i)
1409 z[i] = 0.0;
1410 has_z = 0;
1411 }
1412 if (strncmp(p_in->proj, "ll", 2) == 0) {
1413 if (strncmp(p_out->proj, "ll", 2) == 0) {
1414 DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
1415 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1416 MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
1417 }
1418 else {
1419 DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
1420 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1421 DIVIDE_LOOP(x, y, n, METERS_out);
1422 }
1423 }
1424 else {
1425 if (strncmp(p_out->proj, "ll", 2) == 0) {
1426 MULTIPLY_LOOP(x, y, n, METERS_in);
1427 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1428 MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
1429 }
1430 else {
1431 MULTIPLY_LOOP(x, y, n, METERS_in);
1432 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1433 DIVIDE_LOOP(x, y, n, METERS_out);
1434 }
1435 }
1436 if (!has_z)
1437 G_free(z);
1438
1439 if (ok < 0)
1440 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1441#endif
1442
1443 return ok;
1444}
1445
1446/*
1447 * old API, to be deleted
1448 */
1449
1450/**
1451 * \brief Re-project a point between two co-ordinate systems
1452 *
1453 * This function takes pointers to two pj_info structures as arguments,
1454 * and projects a point between the co-ordinate systems represented by them.
1455 * The easting and northing of the point are contained in two pointers passed
1456 * to the function; these will be overwritten by the co-ordinates of the
1457 * re-projected point.
1458 *
1459 * \param x Pointer to a double containing easting or longitude
1460 * \param y Pointer to a double containing northing or latitude
1461 * \param info_in pointer to pj_info struct for input co-ordinate system
1462 * \param info_out pointer to pj_info struct for output co-ordinate system
1463 *
1464 * \return Return value from PROJ proj_trans() function
1465 **/
1466
1467int pj_do_proj(double *x, double *y, const struct pj_info *info_in,
1468 const struct pj_info *info_out)
1469{
1470 int ok;
1471
1472#ifdef HAVE_PROJ_H
1473 struct pj_info info_trans;
1474 PJ_COORD c;
1475
1476 if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
1477 return -1;
1478 }
1479
1480 METERS_in = info_in->meters;
1481 METERS_out = info_out->meters;
1482
1483 if (strncmp(info_in->proj, "ll", 2) == 0) {
1484 /* convert to radians */
1485 c.lpzt.lam = (*x) / RAD_TO_DEG;
1486 c.lpzt.phi = (*y) / RAD_TO_DEG;
1487 c.lpzt.z = 0;
1488 c.lpzt.t = 0;
1489 c = proj_trans(info_trans.pj, PJ_FWD, c);
1490 ok = proj_errno(info_trans.pj);
1491
1492 if (strncmp(info_out->proj, "ll", 2) == 0) {
1493 /* convert to degrees */
1494 *x = c.lp.lam * RAD_TO_DEG;
1495 *y = c.lp.phi * RAD_TO_DEG;
1496 }
1497 else {
1498 /* convert to map units */
1499 *x = c.xy.x / METERS_out;
1500 *y = c.xy.y / METERS_out;
1501 }
1502 }
1503 else {
1504 /* convert to meters */
1505 c.xyzt.x = *x * METERS_in;
1506 c.xyzt.y = *y * METERS_in;
1507 c.xyzt.z = 0;
1508 c.xyzt.t = 0;
1509 c = proj_trans(info_trans.pj, PJ_FWD, c);
1510 ok = proj_errno(info_trans.pj);
1511
1512 if (strncmp(info_out->proj, "ll", 2) == 0) {
1513 /* convert to degrees */
1514 *x = c.lp.lam * RAD_TO_DEG;
1515 *y = c.lp.phi * RAD_TO_DEG;
1516 }
1517 else {
1518 /* convert to map units */
1519 *x = c.xy.x / METERS_out;
1520 *y = c.xy.y / METERS_out;
1521 }
1522 }
1523 proj_destroy(info_trans.pj);
1524
1525 if (ok < 0) {
1526 G_warning(_("proj_trans() failed: %d"), ok);
1527 }
1528#else
1529 double u, v;
1530 double h = 0.0;
1531
1532 METERS_in = info_in->meters;
1533 METERS_out = info_out->meters;
1534
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;
1542 }
1543 else {
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;
1549 }
1550 }
1551 else {
1552 if (strncmp(info_out->proj, "ll", 2) == 0) {
1553 u = *x * METERS_in;
1554 v = *y * METERS_in;
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;
1558 }
1559 else {
1560 u = *x * METERS_in;
1561 v = *y * METERS_in;
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;
1565 }
1566 }
1567 if (ok < 0) {
1568 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1569 }
1570#endif
1571 return ok;
1572}
1573
1574/**
1575 * \brief Re-project an array of points between two co-ordinate systems with
1576 * optional ellipsoidal height conversion
1577 *
1578 * This function takes pointers to two pj_info structures as arguments,
1579 * and projects an array of points between the co-ordinate systems
1580 * represented by them. Pointers to the three arrays of easting, northing,
1581 * and ellipsoidal height of the point (this one may be NULL) are passed
1582 * to the function; these will be overwritten by the co-ordinates of the
1583 * re-projected points.
1584 *
1585 * \param count Number of points in the arrays to be transformed
1586 * \param x Pointer to an array of type double containing easting or longitude
1587 * \param y Pointer to an array of type double containing northing or latitude
1588 * \param h Pointer to an array of type double containing ellipsoidal height.
1589 * May be null in which case a two-dimensional re-projection will be
1590 * done
1591 * \param info_in pointer to pj_info struct for input co-ordinate system
1592 * \param info_out pointer to pj_info struct for output co-ordinate system
1593 *
1594 * \return Return value from PROJ proj_trans() function
1595 **/
1596
1597int pj_do_transform(int count, double *x, double *y, double *h,
1598 const struct pj_info *info_in,
1599 const struct pj_info *info_out)
1600{
1601 int ok;
1602 int i;
1603 int has_h = 1;
1604
1605#ifdef HAVE_PROJ_H
1606 struct pj_info info_trans;
1607 PJ_COORD c;
1608
1609 if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
1610 return -1;
1611 }
1612
1613 METERS_in = info_in->meters;
1614 METERS_out = info_out->meters;
1615
1616 if (h == NULL) {
1617 h = G_malloc(sizeof *h * count);
1618 /* they say memset is only guaranteed for chars ;-( */
1619 for (i = 0; i < count; ++i)
1620 h[i] = 0.0;
1621 has_h = 0;
1622 }
1623 ok = 0;
1624 if (strncmp(info_in->proj, "ll", 2) == 0) {
1625 c.lpzt.t = 0;
1626 if (strncmp(info_out->proj, "ll", 2) == 0) {
1627 for (i = 0; i < count; i++) {
1628 /* convert to radians */
1629 c.lpzt.lam = x[i] / RAD_TO_DEG;
1630 c.lpzt.phi = y[i] / RAD_TO_DEG;
1631 c.lpzt.z = h[i];
1632 c = proj_trans(info_trans.pj, PJ_FWD, c);
1633 if ((ok = proj_errno(info_trans.pj)) < 0)
1634 break;
1635 /* convert to degrees */
1636 x[i] = c.lp.lam * RAD_TO_DEG;
1637 y[i] = c.lp.phi * RAD_TO_DEG;
1638 }
1639 }
1640 else {
1641 for (i = 0; i < count; i++) {
1642 /* convert to radians */
1643 c.lpzt.lam = x[i] / RAD_TO_DEG;
1644 c.lpzt.phi = y[i] / RAD_TO_DEG;
1645 c.lpzt.z = h[i];
1646 c = proj_trans(info_trans.pj, PJ_FWD, c);
1647 if ((ok = proj_errno(info_trans.pj)) < 0)
1648 break;
1649 /* convert to map units */
1650 x[i] = c.xy.x / METERS_out;
1651 y[i] = c.xy.y / METERS_out;
1652 }
1653 }
1654 }
1655 else {
1656 c.xyzt.t = 0;
1657 if (strncmp(info_out->proj, "ll", 2) == 0) {
1658 for (i = 0; i < count; i++) {
1659 /* convert to meters */
1660 c.xyzt.x = x[i] * METERS_in;
1661 c.xyzt.y = y[i] * METERS_in;
1662 c.xyzt.z = h[i];
1663 c = proj_trans(info_trans.pj, PJ_FWD, c);
1664 if ((ok = proj_errno(info_trans.pj)) < 0)
1665 break;
1666 /* convert to degrees */
1667 x[i] = c.lp.lam * RAD_TO_DEG;
1668 y[i] = c.lp.phi * RAD_TO_DEG;
1669 }
1670 }
1671 else {
1672 for (i = 0; i < count; i++) {
1673 /* convert to meters */
1674 c.xyzt.x = x[i] * METERS_in;
1675 c.xyzt.y = y[i] * METERS_in;
1676 c.xyzt.z = h[i];
1677 c = proj_trans(info_trans.pj, PJ_FWD, c);
1678 if ((ok = proj_errno(info_trans.pj)) < 0)
1679 break;
1680 /* convert to map units */
1681 x[i] = c.xy.x / METERS_out;
1682 y[i] = c.xy.y / METERS_out;
1683 }
1684 }
1685 }
1686 if (!has_h)
1687 G_free(h);
1688 proj_destroy(info_trans.pj);
1689
1690 if (ok < 0) {
1691 G_warning(_("proj_trans() failed: %d"), ok);
1692 }
1693#else
1694 METERS_in = info_in->meters;
1695 METERS_out = info_out->meters;
1696
1697 if (h == NULL) {
1698 h = G_malloc(sizeof *h * count);
1699 /* they say memset is only guaranteed for chars ;-( */
1700 for (i = 0; i < count; ++i)
1701 h[i] = 0.0;
1702 has_h = 0;
1703 }
1704 if (strncmp(info_in->proj, "ll", 2) == 0) {
1705 if (strncmp(info_out->proj, "ll", 2) == 0) {
1706 DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
1707 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1708 MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
1709 }
1710 else {
1711 DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
1712 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1713 DIVIDE_LOOP(x, y, count, METERS_out);
1714 }
1715 }
1716 else {
1717 if (strncmp(info_out->proj, "ll", 2) == 0) {
1718 MULTIPLY_LOOP(x, y, count, METERS_in);
1719 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1720 MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
1721 }
1722 else {
1723 MULTIPLY_LOOP(x, y, count, METERS_in);
1724 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1725 DIVIDE_LOOP(x, y, count, METERS_out);
1726 }
1727 }
1728 if (!has_h)
1729 G_free(h);
1730
1731 if (ok < 0) {
1732 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1733 }
1734#endif
1735 return ok;
1736}
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:150
int G_asprintf(char **out, const char *fmt,...)
Definition: asprintf.c:69
#define NULL
Definition: ccmath.h:32
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:66
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...
Definition: do_proj.c:985
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.
Definition: do_proj.c:1467
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...
Definition: do_proj.c:1597
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 ...
Definition: do_proj.c:1204
#define MULTIPLY_LOOP(x, y, c, m)
Definition: do_proj.c:26
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.
Definition: do_proj.c:434
#define DIVIDE_LOOP(x, y, c, m)
Definition: do_proj.c:34
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...
Definition: get_proj.c:493
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
Definition: gis/error.c:130
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
void G_get_set_window(struct Cell_head *window)
Get the current working window (region)
int count
char * G_store_lower(const char *s)
Copy string to allocated memory and convert copied string to lower case.
Definition: strings.c:141
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:87
char * G_store_upper(const char *s)
Copy string to allocated memory and convert copied string to upper case.
Definition: strings.c:117
#define x