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