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