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