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