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