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