GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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 * using PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION
702 * this can cause problems and artefacts
703 * switch to PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT
704 * in case of problems
705 * but results can be different from gdalwarp:
706 * shifted geolocation in some areas
707 * in these cases there is no right or wrong,
708 * different pipelines are all regarded as valid by PROJ
709 * depending on the area of interest
710 *
711 * see also:
712 * OGRProjCT::ListCoordinateOperations() in GDAL ogr/ogrct.cpp
713 * create_operation_to_geog_crs() in PROJ src/4D_api.cpp
714 * proj_create_crs_to_crs_from_pj() in PROJ src/4D_api.cpp
715 * proj_operation_factory_context_set_spatial_criterion() in PROJ
716 * src/iso19111/c_api.cpp
717 * */
718
719 /* now use the current region as area of interest */
720 operation_ctx =
721 proj_create_operation_factory_context(PJ_DEFAULT_CTX, NULL);
722 proj_operation_factory_context_set_area_of_interest(
723 PJ_DEFAULT_CTX, operation_ctx, xmin, ymin, xmax, ymax);
724 proj_operation_factory_context_set_spatial_criterion(
725 PJ_DEFAULT_CTX, operation_ctx,
726 PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
727 /* from GDAL OGRProjCT::ListCoordinateOperations() */
728 proj_operation_factory_context_set_grid_availability_use(
729 PJ_DEFAULT_CTX, operation_ctx,
730#if PROJ_VERSION_NUM >= 7000000
731 proj_context_is_network_enabled(PJ_DEFAULT_CTX)
732 ? PROJ_GRID_AVAILABILITY_KNOWN_AVAILABLE
733 :
734#endif
735 PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
736
737 /* The operations are sorted with the most relevant ones first:
738 * by descending area (intersection of the transformation area
739 * with the area of interest, or intersection of the
740 * transformation with the area of use of the CRS),
741 * and by increasing accuracy.
742 * Operations with unknown accuracy are sorted last,
743 * whatever their area.
744 */
745 op_list = proj_create_operations(PJ_DEFAULT_CTX, in_pj, out_pj,
746 operation_ctx);
747 proj_operation_factory_context_destroy(operation_ctx);
748 op_count_area = 0;
749 if (op_list)
750 op_count_area = proj_list_get_count(op_list);
751 if (op_count_area == 0) {
752 /* no operations */
753 info_trans->pj = NULL;
754 }
755 else if (op_count_area == 1) {
756 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
757 }
758 else { /* op_count_area > 1 */
759 /* can't use pj_create_prepared_operations()
760 * this is a PROJ-internal function
761 * trust the sorting of PROJ and use the first one */
762 info_trans->pj = proj_list_get(PJ_DEFAULT_CTX, op_list, 0);
763 }
764 if (op_list)
765 proj_list_destroy(op_list);
766
767 /* try proj_create_crs_to_crs() */
768 /*
769 G_debug(1, "trying %s to %s", indef, outdef);
770 */
771
772 /* proj_create_crs_to_crs() does not work because it calls
773 * proj_create_crs_to_crs_from_pj() which calls
774 * proj_operation_factory_context_set_spatial_criterion()
775 * with PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION
776 * instead of
777 * PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT
778 *
779 * fixed in PROJ master, probably available with PROJ 7.3.x */
780
781 /*
782 info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
783 indef,
784 outdef,
785 pj_area);
786 */
787
788 if (in_pj)
789 proj_destroy(in_pj);
790 if (out_pj)
791 proj_destroy(out_pj);
792
793 if (info_trans->pj) {
794 const char *projstr;
795 PJ *pj_norm = NULL;
796
797 G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ%d",
798 PROJ_VERSION_MAJOR);
799
800 projstr =
801 proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);
802
803 info_trans->def = G_store(projstr);
804
805 if (projstr) {
806 /* make sure axis order is easting, northing
807 * proj_normalize_for_visualization() requires
808 * source and target CRS
809 * -> does not work with ll equivalent of input:
810 * no target CRS in +proj=pipeline +step +inv %s */
811 pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX,
812 info_trans->pj);
813
814 if (!pj_norm) {
815 G_warning(
816 _("proj_normalize_for_visualization() failed for '%s'"),
817 info_trans->def);
818 }
819 else {
820 projstr =
821 proj_as_proj_string(NULL, pj_norm, PJ_PROJ_5, NULL);
822 if (projstr && *projstr) {
823 proj_destroy(info_trans->pj);
824 info_trans->pj = pj_norm;
825 info_trans->def = G_store(projstr);
826 }
827 else {
828 proj_destroy(pj_norm);
829 G_warning(_("No PROJ definition for normalized version "
830 "of '%s'"),
831 info_trans->def);
832 }
833 }
834 G_important_message(_("Selected PROJ pipeline:"));
835 G_important_message(_("%s"), info_trans->def);
836 G_important_message("************************");
837 }
838 else {
839 proj_destroy(info_trans->pj);
840 info_trans->pj = NULL;
841 }
842 }
843
844 if (pj_area)
845 proj_area_destroy(pj_area);
846
847 if (insrid)
848 G_free(insrid);
849 if (outsrid)
850 G_free(outsrid);
851 G_free(indef);
852 G_free(outdef);
853 }
854 if (info_trans->pj == NULL) {
855 G_warning(_("proj_create() failed for '%s'"), info_trans->def);
856
857 return -1;
858 }
859#else /* PROJ 5 */
860 info_trans->pj = NULL;
861 info_trans->meters = 1.;
862 info_trans->zone = 0;
863 sprintf(info_trans->proj, "pipeline");
864
865 /* user-provided pipeline */
866 if (info_trans->def) {
867 /* create a pj from user-defined transformation pipeline */
868 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
869 if (info_trans->pj == NULL) {
870 G_warning(_("proj_create() failed for '%s'"), info_trans->def);
871
872 return -1;
873 }
874 }
875 /* if no output CRS is defined,
876 * assume info_out to be ll equivalent of info_in */
877 else if (info_out->pj == NULL) {
878 G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",
879 info_in->def);
880 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
881 if (info_trans->pj == NULL) {
882 G_warning(_("proj_create() failed for '%s'"), info_trans->def);
883
884 return -1;
885 }
886 }
887 else if (info_in->def && info_out->pj && info_out->def) {
888 char *indef = NULL, *outdef = NULL;
889 char *insrid = NULL, *outsrid = NULL;
890
891 /* PROJ5: EPSG must be lowercase epsg */
892 if (info_in->srid) {
893 if (strncmp(info_in->srid, "EPSG", 4) == 0)
894 insrid = G_store_lower(info_in->srid);
895 else
896 insrid = G_store(info_in->srid);
897 }
898
899 if (info_out->pj && info_out->srid) {
900 if (strncmp(info_out->srid, "EPSG", 4) == 0)
901 outsrid = G_store_lower(info_out->srid);
902 else
903 outsrid = G_store(info_out->srid);
904 }
905
906 info_trans->pj = NULL;
907
908 if (insrid && outsrid) {
909 G_asprintf(&indef, "%s", insrid);
910 G_asprintf(&outdef, "%s", outsrid);
911
912 G_asprintf(&(info_trans->def),
913 "+proj=pipeline +step +inv +init=%s +step +init=%s",
914 indef, outdef);
915
916 /* try proj_create_crs_to_crs() */
917 info_trans->pj =
918 proj_create_crs_to_crs(PJ_DEFAULT_CTX, indef, outdef, NULL);
919 }
920
921 if (info_trans->pj) {
922 G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ5");
923 }
924 else {
925 if (indef) {
926 G_free(indef);
927 indef = NULL;
928 }
929 if (insrid) {
930 G_asprintf(&indef, "+init=%s", insrid);
931 }
932 else {
933 G_asprintf(&indef, "%s", info_in->def);
934 }
935
936 if (outdef) {
937 G_free(outdef);
938 outdef = NULL;
939 }
940 if (outsrid) {
941 G_asprintf(&outdef, "+init=%s", outsrid);
942 }
943 else {
944 G_asprintf(&outdef, "%s", info_out->def);
945 }
946
947 /* try proj_create() with +proj=pipeline +step +inv %s +step %s" */
948 G_asprintf(&(info_trans->def),
949 "+proj=pipeline +step +inv %s +step %s", indef, outdef);
950
951 info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
952 }
953 if (insrid)
954 G_free(insrid);
955 if (outsrid)
956 G_free(outsrid);
957 G_free(indef);
958 G_free(outdef);
959 }
960 if (info_trans->pj == NULL) {
961 G_warning(_("proj_create() failed for '%s'"), info_trans->def);
962
963 return -1;
964 }
965
966#endif
967#else /* PROJ 4 */
968 if (info_out->pj == NULL) {
969 if (GPJ_get_equivalent_latlong(info_out, info_in) < 0) {
970 G_warning(_("Unable to create latlong equivalent for '%s'"),
971 info_in->def);
972
973 return -1;
974 }
975 }
976#endif
977
978 return 1;
979}
980
981/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
982
983/**
984 * \brief Re-project a point between two co-ordinate systems using a
985 * transformation object prepared with GPJ_prepare_pj()
986 *
987 * This function takes pointers to three pj_info structures as arguments,
988 * and projects a point between the input and output co-ordinate system.
989 * The pj_info structure info_trans must have been initialized with
990 * GPJ_init_transform().
991 * The direction determines if a point is projected from input CRS to
992 * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
993 * The easting, northing, and height of the point are contained in the
994 * pointers passed to the function; these will be overwritten by the
995 * coordinates of the transformed point.
996 *
997 * \param info_in pointer to pj_info struct for input co-ordinate system
998 * \param info_out pointer to pj_info struct for output co-ordinate system
999 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
1000 *5+) \param dir direction of the transformation (PJ_FWD or PJ_INV) \param x
1001 *Pointer to a double containing easting or longitude \param y Pointer to a
1002 *double containing northing or latitude \param z Pointer to a double containing
1003 *height, or NULL
1004 *
1005 * \return Return value from PROJ proj_trans() function
1006 **/
1007
1008int GPJ_transform(const struct pj_info *info_in, const struct pj_info *info_out,
1009 const struct pj_info *info_trans, int dir, double *x,
1010 double *y, double *z)
1011{
1012 int ok = 0;
1013
1014#ifdef HAVE_PROJ_H
1015 /* PROJ 5+ variant */
1016 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1017 PJ_COORD c;
1018
1019 if (info_in->pj == NULL)
1020 G_fatal_error(_("No input projection"));
1021
1022 if (info_trans->pj == NULL)
1023 G_fatal_error(_("No transformation object"));
1024
1025 in_deg2rad = out_rad2deg = 1;
1026 if (dir == PJ_FWD) {
1027 /* info_in -> info_out */
1028 METERS_in = info_in->meters;
1029 in_is_ll = !strncmp(info_in->proj, "ll", 2);
1030#if PROJ_VERSION_MAJOR >= 6
1031 /* PROJ 6+: conversion to radians is not always needed:
1032 * if proj_angular_input(info_trans->pj, dir) == 1
1033 * -> convert from degrees to radians */
1034 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1035 in_deg2rad = 0;
1036 }
1037#endif
1038 if (info_out->pj) {
1039 METERS_out = info_out->meters;
1040 out_is_ll = !strncmp(info_out->proj, "ll", 2);
1041#if PROJ_VERSION_MAJOR >= 6
1042 /* PROJ 6+: conversion to radians is not always needed:
1043 * if proj_angular_input(info_trans->pj, dir) == 1
1044 * -> convert from degrees to radians */
1045 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1046 out_rad2deg = 0;
1047 }
1048#endif
1049 }
1050 else {
1051 METERS_out = 1.0;
1052 out_is_ll = 1;
1053 }
1054 }
1055 else {
1056 /* info_out -> info_in */
1057 METERS_out = info_in->meters;
1058 out_is_ll = !strncmp(info_in->proj, "ll", 2);
1059#if PROJ_VERSION_MAJOR >= 6
1060 /* PROJ 6+: conversion to radians is not always needed:
1061 * if proj_angular_input(info_trans->pj, dir) == 1
1062 * -> convert from degrees to radians */
1063 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1064 out_rad2deg = 0;
1065 }
1066#endif
1067 if (info_out->pj) {
1068 METERS_in = info_out->meters;
1069 in_is_ll = !strncmp(info_out->proj, "ll", 2);
1070#if PROJ_VERSION_MAJOR >= 6
1071 /* PROJ 6+: conversion to radians is not always needed:
1072 * if proj_angular_input(info_trans->pj, dir) == 1
1073 * -> convert from degrees to radians */
1074 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1075 in_deg2rad = 0;
1076 }
1077#endif
1078 }
1079 else {
1080 METERS_in = 1.0;
1081 in_is_ll = 1;
1082 }
1083 }
1084
1085 /* prepare */
1086 if (in_is_ll) {
1087 if (in_deg2rad) {
1088 /* convert degrees to radians */
1089 c.lpzt.lam = (*x) / RAD_TO_DEG;
1090 c.lpzt.phi = (*y) / RAD_TO_DEG;
1091 }
1092 else {
1093 c.lpzt.lam = (*x);
1094 c.lpzt.phi = (*y);
1095 }
1096 c.lpzt.z = 0;
1097 if (z)
1098 c.lpzt.z = *z;
1099 c.lpzt.t = 0;
1100 }
1101 else {
1102 /* convert to meters */
1103 c.xyzt.x = *x * METERS_in;
1104 c.xyzt.y = *y * METERS_in;
1105 c.xyzt.z = 0;
1106 if (z)
1107 c.xyzt.z = *z;
1108 c.xyzt.t = 0;
1109 }
1110
1111 G_debug(1, "c.xyzt.x: %g", c.xyzt.x);
1112 G_debug(1, "c.xyzt.y: %g", c.xyzt.y);
1113 G_debug(1, "c.xyzt.z: %g", c.xyzt.z);
1114
1115 /* transform */
1116 c = proj_trans(info_trans->pj, dir, c);
1117 ok = proj_errno(info_trans->pj);
1118
1119 if (ok < 0) {
1120 G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
1121 return ok;
1122 }
1123
1124 /* output */
1125 if (out_is_ll) {
1126 /* convert to degrees */
1127 if (out_rad2deg) {
1128 /* convert radians to degrees */
1129 *x = c.lp.lam * RAD_TO_DEG;
1130 *y = c.lp.phi * RAD_TO_DEG;
1131 }
1132 else {
1133 *x = c.lp.lam;
1134 *y = c.lp.phi;
1135 }
1136 if (z)
1137 *z = c.lpzt.z;
1138 }
1139 else {
1140 /* convert to map units */
1141 *x = c.xyzt.x / METERS_out;
1142 *y = c.xyzt.y / METERS_out;
1143 if (z)
1144 *z = c.xyzt.z;
1145 }
1146#else
1147 /* PROJ 4 variant */
1148 double u, v;
1149 double h = 0.0;
1150 const struct pj_info *p_in, *p_out;
1151
1152 if (info_out == NULL)
1153 G_fatal_error(_("No output projection"));
1154
1155 if (dir == PJ_FWD) {
1156 p_in = info_in;
1157 p_out = info_out;
1158 }
1159 else {
1160 p_in = info_out;
1161 p_out = info_in;
1162 }
1163
1164 METERS_in = p_in->meters;
1165 METERS_out = p_out->meters;
1166
1167 if (z)
1168 h = *z;
1169
1170 if (strncmp(p_in->proj, "ll", 2) == 0) {
1171 u = (*x) / RAD_TO_DEG;
1172 v = (*y) / RAD_TO_DEG;
1173 }
1174 else {
1175 u = *x * METERS_in;
1176 v = *y * METERS_in;
1177 }
1178
1179 ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1180
1181 if (ok < 0) {
1182 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1183 return ok;
1184 }
1185
1186 if (strncmp(p_out->proj, "ll", 2) == 0) {
1187 *x = u * RAD_TO_DEG;
1188 *y = v * RAD_TO_DEG;
1189 }
1190 else {
1191 *x = u / METERS_out;
1192 *y = v / METERS_out;
1193 }
1194 if (z)
1195 *z = h;
1196#endif
1197
1198 return ok;
1199}
1200
1201/**
1202 * \brief Re-project an array of points between two co-ordinate systems
1203 * using a transformation object prepared with GPJ_prepare_pj()
1204 *
1205 * This function takes pointers to three pj_info structures as arguments,
1206 * and projects an array of pointd between the input and output
1207 * co-ordinate system. The pj_info structure info_trans must have been
1208 * initialized with GPJ_init_transform().
1209 * The direction determines if a point is projected from input CRS to
1210 * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
1211 * The easting, northing, and height of the point are contained in the
1212 * pointers passed to the function; these will be overwritten by the
1213 * coordinates of the transformed point.
1214 *
1215 * \param info_in pointer to pj_info struct for input co-ordinate system
1216 * \param info_out pointer to pj_info struct for output co-ordinate system
1217 * \param info_trans pointer to pj_info struct for a transformation object (PROJ
1218 *5+) \param dir direction of the transformation (PJ_FWD or PJ_INV) \param x
1219 *pointer to an array of type double containing easting or longitude \param y
1220 *pointer to an array of type double containing northing or latitude \param z
1221 *pointer to an array of type double containing height, or NULL \param n number
1222 *of points in the arrays to be transformed
1223 *
1224 * \return Return value from PROJ proj_trans() function
1225 **/
1226
1227int GPJ_transform_array(const struct pj_info *info_in,
1228 const struct pj_info *info_out,
1229 const struct pj_info *info_trans, int dir, double *x,
1230 double *y, double *z, int n)
1231{
1232 int ok;
1233 int i;
1234 int has_z = 1;
1235
1236#ifdef HAVE_PROJ_H
1237 /* PROJ 5+ variant */
1238 int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1239 PJ_COORD c;
1240
1241 if (info_trans->pj == NULL)
1242 G_fatal_error(_("No transformation object"));
1243
1244 in_deg2rad = out_rad2deg = 1;
1245 if (dir == PJ_FWD) {
1246 /* info_in -> info_out */
1247 METERS_in = info_in->meters;
1248 in_is_ll = !strncmp(info_in->proj, "ll", 2);
1249#if PROJ_VERSION_MAJOR >= 6
1250 /* PROJ 6+: conversion to radians is not always needed:
1251 * if proj_angular_input(info_trans->pj, dir) == 1
1252 * -> convert from degrees to radians */
1253 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1254 in_deg2rad = 0;
1255 }
1256#endif
1257 if (info_out->pj) {
1258 METERS_out = info_out->meters;
1259 out_is_ll = !strncmp(info_out->proj, "ll", 2);
1260#if PROJ_VERSION_MAJOR >= 6
1261 /* PROJ 6+: conversion to radians is not always needed:
1262 * if proj_angular_input(info_trans->pj, dir) == 1
1263 * -> convert from degrees to radians */
1264 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1265 out_rad2deg = 0;
1266 }
1267#endif
1268 }
1269 else {
1270 METERS_out = 1.0;
1271 out_is_ll = 1;
1272 }
1273 }
1274 else {
1275 /* info_out -> info_in */
1276 METERS_out = info_in->meters;
1277 out_is_ll = !strncmp(info_in->proj, "ll", 2);
1278#if PROJ_VERSION_MAJOR >= 6
1279 /* PROJ 6+: conversion to radians is not always needed:
1280 * if proj_angular_input(info_trans->pj, dir) == 1
1281 * -> convert from degrees to radians */
1282 if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1283 out_rad2deg = 0;
1284 }
1285#endif
1286 if (info_out->pj) {
1287 METERS_in = info_out->meters;
1288 in_is_ll = !strncmp(info_out->proj, "ll", 2);
1289#if PROJ_VERSION_MAJOR >= 6
1290 /* PROJ 6+: conversion to degrees is not always needed:
1291 * if proj_angular_output(info_trans->pj, dir) == 1
1292 * -> convert from degrees to radians */
1293 if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1294 in_deg2rad = 0;
1295 }
1296#endif
1297 }
1298 else {
1299 METERS_in = 1.0;
1300 in_is_ll = 1;
1301 }
1302 }
1303
1304 if (z == NULL) {
1305 z = G_malloc(sizeof(double) * n);
1306 /* they say memset is only guaranteed for chars ;-( */
1307 for (i = 0; i < n; i++)
1308 z[i] = 0.0;
1309 has_z = 0;
1310 }
1311 ok = 0;
1312 if (in_is_ll) {
1313 c.lpzt.t = 0;
1314 if (out_is_ll) {
1315 /* what is more costly ?
1316 * calling proj_trans for each point
1317 * or having three loops over all points ?
1318 * proj_trans_array() itself calls proj_trans() in a loop
1319 * -> one loop over all points is better than
1320 * three loops over all points
1321 */
1322 for (i = 0; i < n; i++) {
1323 if (in_deg2rad) {
1324 /* convert degrees to radians */
1325 c.lpzt.lam = (*x) / RAD_TO_DEG;
1326 c.lpzt.phi = (*y) / RAD_TO_DEG;
1327 }
1328 else {
1329 c.lpzt.lam = (*x);
1330 c.lpzt.phi = (*y);
1331 }
1332 c.lpzt.z = z[i];
1333 c = proj_trans(info_trans->pj, dir, c);
1334 if ((ok = proj_errno(info_trans->pj)) < 0)
1335 break;
1336 if (out_rad2deg) {
1337 /* convert radians to degrees */
1338 x[i] = c.lp.lam * RAD_TO_DEG;
1339 y[i] = c.lp.phi * RAD_TO_DEG;
1340 }
1341 else {
1342 x[i] = c.lp.lam;
1343 y[i] = c.lp.phi;
1344 }
1345 }
1346 }
1347 else {
1348 for (i = 0; i < n; i++) {
1349 if (in_deg2rad) {
1350 /* convert degrees to radians */
1351 c.lpzt.lam = (*x) / RAD_TO_DEG;
1352 c.lpzt.phi = (*y) / RAD_TO_DEG;
1353 }
1354 else {
1355 c.lpzt.lam = (*x);
1356 c.lpzt.phi = (*y);
1357 }
1358 c.lpzt.z = z[i];
1359 c = proj_trans(info_trans->pj, dir, c);
1360 if ((ok = proj_errno(info_trans->pj)) < 0)
1361 break;
1362
1363 /* convert to map units */
1364 x[i] = c.xy.x / METERS_out;
1365 y[i] = c.xy.y / METERS_out;
1366 }
1367 }
1368 }
1369 else {
1370 c.xyzt.t = 0;
1371 if (out_is_ll) {
1372 for (i = 0; i < n; i++) {
1373 /* convert to meters */
1374 c.xyzt.x = x[i] * METERS_in;
1375 c.xyzt.y = y[i] * METERS_in;
1376 c.xyzt.z = z[i];
1377 c = proj_trans(info_trans->pj, dir, c);
1378 if ((ok = proj_errno(info_trans->pj)) < 0)
1379 break;
1380 if (out_rad2deg) {
1381 /* convert radians to degrees */
1382 x[i] = c.lp.lam * RAD_TO_DEG;
1383 y[i] = c.lp.phi * RAD_TO_DEG;
1384 }
1385 else {
1386 x[i] = c.lp.lam;
1387 y[i] = c.lp.phi;
1388 }
1389 }
1390 }
1391 else {
1392 for (i = 0; i < n; i++) {
1393 /* convert to meters */
1394 c.xyzt.x = x[i] * METERS_in;
1395 c.xyzt.y = y[i] * METERS_in;
1396 c.xyzt.z = z[i];
1397 c = proj_trans(info_trans->pj, dir, c);
1398 if ((ok = proj_errno(info_trans->pj)) < 0)
1399 break;
1400 /* convert to map units */
1401 x[i] = c.xy.x / METERS_out;
1402 y[i] = c.xy.y / METERS_out;
1403 }
1404 }
1405 }
1406 if (!has_z)
1407 G_free(z);
1408
1409 if (ok < 0) {
1410 G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
1411 }
1412#else
1413 /* PROJ 4 variant */
1414 const struct pj_info *p_in, *p_out;
1415
1416 if (dir == PJ_FWD) {
1417 p_in = info_in;
1418 p_out = info_out;
1419 }
1420 else {
1421 p_in = info_out;
1422 p_out = info_in;
1423 }
1424
1425 METERS_in = p_in->meters;
1426 METERS_out = p_out->meters;
1427
1428 if (z == NULL) {
1429 z = G_malloc(sizeof(double) * n);
1430 /* they say memset is only guaranteed for chars ;-( */
1431 for (i = 0; i < n; ++i)
1432 z[i] = 0.0;
1433 has_z = 0;
1434 }
1435 if (strncmp(p_in->proj, "ll", 2) == 0) {
1436 if (strncmp(p_out->proj, "ll", 2) == 0) {
1437 DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
1438 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1439 MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
1440 }
1441 else {
1442 DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
1443 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1444 DIVIDE_LOOP(x, y, n, METERS_out);
1445 }
1446 }
1447 else {
1448 if (strncmp(p_out->proj, "ll", 2) == 0) {
1449 MULTIPLY_LOOP(x, y, n, METERS_in);
1450 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1451 MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
1452 }
1453 else {
1454 MULTIPLY_LOOP(x, y, n, METERS_in);
1455 ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1456 DIVIDE_LOOP(x, y, n, METERS_out);
1457 }
1458 }
1459 if (!has_z)
1460 G_free(z);
1461
1462 if (ok < 0)
1463 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1464#endif
1465
1466 return ok;
1467}
1468
1469/*
1470 * old API, to be deleted
1471 */
1472
1473/**
1474 * \brief Re-project a point between two co-ordinate systems
1475 *
1476 * This function takes pointers to two pj_info structures as arguments,
1477 * and projects a point between the co-ordinate systems represented by them.
1478 * The easting and northing of the point are contained in two pointers passed
1479 * to the function; these will be overwritten by the co-ordinates of the
1480 * re-projected point.
1481 *
1482 * \param x Pointer to a double containing easting or longitude
1483 * \param y Pointer to a double containing northing or latitude
1484 * \param info_in pointer to pj_info struct for input co-ordinate system
1485 * \param info_out pointer to pj_info struct for output co-ordinate system
1486 *
1487 * \return Return value from PROJ proj_trans() function
1488 **/
1489
1490int pj_do_proj(double *x, double *y, const struct pj_info *info_in,
1491 const struct pj_info *info_out)
1492{
1493 int ok;
1494
1495#ifdef HAVE_PROJ_H
1496 struct pj_info info_trans;
1497 PJ_COORD c;
1498
1499 if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
1500 return -1;
1501 }
1502
1503 METERS_in = info_in->meters;
1504 METERS_out = info_out->meters;
1505
1506 if (strncmp(info_in->proj, "ll", 2) == 0) {
1507 /* convert to radians */
1508 c.lpzt.lam = (*x) / RAD_TO_DEG;
1509 c.lpzt.phi = (*y) / RAD_TO_DEG;
1510 c.lpzt.z = 0;
1511 c.lpzt.t = 0;
1512 c = proj_trans(info_trans.pj, PJ_FWD, c);
1513 ok = proj_errno(info_trans.pj);
1514
1515 if (strncmp(info_out->proj, "ll", 2) == 0) {
1516 /* convert to degrees */
1517 *x = c.lp.lam * RAD_TO_DEG;
1518 *y = c.lp.phi * RAD_TO_DEG;
1519 }
1520 else {
1521 /* convert to map units */
1522 *x = c.xy.x / METERS_out;
1523 *y = c.xy.y / METERS_out;
1524 }
1525 }
1526 else {
1527 /* convert to meters */
1528 c.xyzt.x = *x * METERS_in;
1529 c.xyzt.y = *y * METERS_in;
1530 c.xyzt.z = 0;
1531 c.xyzt.t = 0;
1532 c = proj_trans(info_trans.pj, PJ_FWD, c);
1533 ok = proj_errno(info_trans.pj);
1534
1535 if (strncmp(info_out->proj, "ll", 2) == 0) {
1536 /* convert to degrees */
1537 *x = c.lp.lam * RAD_TO_DEG;
1538 *y = c.lp.phi * RAD_TO_DEG;
1539 }
1540 else {
1541 /* convert to map units */
1542 *x = c.xy.x / METERS_out;
1543 *y = c.xy.y / METERS_out;
1544 }
1545 }
1546 proj_destroy(info_trans.pj);
1547
1548 if (ok < 0) {
1549 G_warning(_("proj_trans() failed: %d"), ok);
1550 }
1551#else
1552 double u, v;
1553 double h = 0.0;
1554
1555 METERS_in = info_in->meters;
1556 METERS_out = info_out->meters;
1557
1558 if (strncmp(info_in->proj, "ll", 2) == 0) {
1559 if (strncmp(info_out->proj, "ll", 2) == 0) {
1560 u = (*x) / RAD_TO_DEG;
1561 v = (*y) / RAD_TO_DEG;
1562 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1563 *x = u * RAD_TO_DEG;
1564 *y = v * RAD_TO_DEG;
1565 }
1566 else {
1567 u = (*x) / RAD_TO_DEG;
1568 v = (*y) / RAD_TO_DEG;
1569 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1570 *x = u / METERS_out;
1571 *y = v / METERS_out;
1572 }
1573 }
1574 else {
1575 if (strncmp(info_out->proj, "ll", 2) == 0) {
1576 u = *x * METERS_in;
1577 v = *y * METERS_in;
1578 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1579 *x = u * RAD_TO_DEG;
1580 *y = v * RAD_TO_DEG;
1581 }
1582 else {
1583 u = *x * METERS_in;
1584 v = *y * METERS_in;
1585 ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1586 *x = u / METERS_out;
1587 *y = v / METERS_out;
1588 }
1589 }
1590 if (ok < 0) {
1591 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1592 }
1593#endif
1594 return ok;
1595}
1596
1597/**
1598 * \brief Re-project an array of points between two co-ordinate systems with
1599 * optional ellipsoidal height conversion
1600 *
1601 * This function takes pointers to two pj_info structures as arguments,
1602 * and projects an array of points between the co-ordinate systems
1603 * represented by them. Pointers to the three arrays of easting, northing,
1604 * and ellipsoidal height of the point (this one may be NULL) are passed
1605 * to the function; these will be overwritten by the co-ordinates of the
1606 * re-projected points.
1607 *
1608 * \param count Number of points in the arrays to be transformed
1609 * \param x Pointer to an array of type double containing easting or longitude
1610 * \param y Pointer to an array of type double containing northing or latitude
1611 * \param h Pointer to an array of type double containing ellipsoidal height.
1612 * May be null in which case a two-dimensional re-projection will be
1613 * done
1614 * \param info_in pointer to pj_info struct for input co-ordinate system
1615 * \param info_out pointer to pj_info struct for output co-ordinate system
1616 *
1617 * \return Return value from PROJ proj_trans() function
1618 **/
1619
1620int pj_do_transform(int count, double *x, double *y, double *h,
1621 const struct pj_info *info_in,
1622 const struct pj_info *info_out)
1623{
1624 int ok;
1625 int i;
1626 int has_h = 1;
1627
1628#ifdef HAVE_PROJ_H
1629 struct pj_info info_trans;
1630 PJ_COORD c;
1631
1632 if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
1633 return -1;
1634 }
1635
1636 METERS_in = info_in->meters;
1637 METERS_out = info_out->meters;
1638
1639 if (h == NULL) {
1640 h = G_malloc(sizeof *h * count);
1641 /* they say memset is only guaranteed for chars ;-( */
1642 for (i = 0; i < count; ++i)
1643 h[i] = 0.0;
1644 has_h = 0;
1645 }
1646 ok = 0;
1647 if (strncmp(info_in->proj, "ll", 2) == 0) {
1648 c.lpzt.t = 0;
1649 if (strncmp(info_out->proj, "ll", 2) == 0) {
1650 for (i = 0; i < count; i++) {
1651 /* convert to radians */
1652 c.lpzt.lam = x[i] / RAD_TO_DEG;
1653 c.lpzt.phi = y[i] / RAD_TO_DEG;
1654 c.lpzt.z = h[i];
1655 c = proj_trans(info_trans.pj, PJ_FWD, c);
1656 if ((ok = proj_errno(info_trans.pj)) < 0)
1657 break;
1658 /* convert to degrees */
1659 x[i] = c.lp.lam * RAD_TO_DEG;
1660 y[i] = c.lp.phi * RAD_TO_DEG;
1661 }
1662 }
1663 else {
1664 for (i = 0; i < count; i++) {
1665 /* convert to radians */
1666 c.lpzt.lam = x[i] / RAD_TO_DEG;
1667 c.lpzt.phi = y[i] / RAD_TO_DEG;
1668 c.lpzt.z = h[i];
1669 c = proj_trans(info_trans.pj, PJ_FWD, c);
1670 if ((ok = proj_errno(info_trans.pj)) < 0)
1671 break;
1672 /* convert to map units */
1673 x[i] = c.xy.x / METERS_out;
1674 y[i] = c.xy.y / METERS_out;
1675 }
1676 }
1677 }
1678 else {
1679 c.xyzt.t = 0;
1680 if (strncmp(info_out->proj, "ll", 2) == 0) {
1681 for (i = 0; i < count; i++) {
1682 /* convert to meters */
1683 c.xyzt.x = x[i] * METERS_in;
1684 c.xyzt.y = y[i] * METERS_in;
1685 c.xyzt.z = h[i];
1686 c = proj_trans(info_trans.pj, PJ_FWD, c);
1687 if ((ok = proj_errno(info_trans.pj)) < 0)
1688 break;
1689 /* convert to degrees */
1690 x[i] = c.lp.lam * RAD_TO_DEG;
1691 y[i] = c.lp.phi * RAD_TO_DEG;
1692 }
1693 }
1694 else {
1695 for (i = 0; i < count; i++) {
1696 /* convert to meters */
1697 c.xyzt.x = x[i] * METERS_in;
1698 c.xyzt.y = y[i] * METERS_in;
1699 c.xyzt.z = h[i];
1700 c = proj_trans(info_trans.pj, PJ_FWD, c);
1701 if ((ok = proj_errno(info_trans.pj)) < 0)
1702 break;
1703 /* convert to map units */
1704 x[i] = c.xy.x / METERS_out;
1705 y[i] = c.xy.y / METERS_out;
1706 }
1707 }
1708 }
1709 if (!has_h)
1710 G_free(h);
1711 proj_destroy(info_trans.pj);
1712
1713 if (ok < 0) {
1714 G_warning(_("proj_trans() failed: %d"), ok);
1715 }
1716#else
1717 METERS_in = info_in->meters;
1718 METERS_out = info_out->meters;
1719
1720 if (h == NULL) {
1721 h = G_malloc(sizeof *h * count);
1722 /* they say memset is only guaranteed for chars ;-( */
1723 for (i = 0; i < count; ++i)
1724 h[i] = 0.0;
1725 has_h = 0;
1726 }
1727 if (strncmp(info_in->proj, "ll", 2) == 0) {
1728 if (strncmp(info_out->proj, "ll", 2) == 0) {
1729 DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
1730 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1731 MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
1732 }
1733 else {
1734 DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
1735 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1736 DIVIDE_LOOP(x, y, count, METERS_out);
1737 }
1738 }
1739 else {
1740 if (strncmp(info_out->proj, "ll", 2) == 0) {
1741 MULTIPLY_LOOP(x, y, count, METERS_in);
1742 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1743 MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
1744 }
1745 else {
1746 MULTIPLY_LOOP(x, y, count, METERS_in);
1747 ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1748 DIVIDE_LOOP(x, y, count, METERS_out);
1749 }
1750 }
1751 if (!has_h)
1752 G_free(h);
1753
1754 if (ok < 0) {
1755 G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1756 }
1757#endif
1758 return ok;
1759}
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:1008
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:1490
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:1620
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:1227
#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