GRASS GIS 7 Programmer's Manual  7.8.1(2019)-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  PJ *source_crs;
382 
383  /* Even Rouault:
384  * if info_in->def contains a +towgs84/+nadgrids clause,
385  * this pipeline would apply it, whereas you probably only want
386  * the reverse projection, and no datum shift.
387  * The easiest would probably to mess up with the PROJ string.
388  * Otherwise with the PROJ API, you could
389  * instanciate a PJ object from the string,
390  * check if it is a BoundCRS with proj_get_source_crs(),
391  * and in that case, take the source CRS with proj_get_source_crs(),
392  * and do the inverse transform on it */
393 
394  source_crs = proj_get_source_crs(NULL, info_in->pj);
395  if (source_crs) {
396  projstr = proj_as_proj_string(NULL, source_crs, PJ_PROJ_5, NULL);
397  if (projstr)
398  indef = G_store(projstr);
399  proj_destroy(source_crs);
400  }
401  if (indef == NULL)
402  indef = G_store(info_in->def);
403 
404  /* what about axis order?
405  * is it always enu?
406  * probably yes, as long as there is no +proj=axisswap step */
407  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",
408  indef);
409  info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
410  if (info_trans->pj == NULL) {
411  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
412  G_free(indef);
413 
414  return -1;
415  }
416  projstr = proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);
417  if (projstr == NULL) {
418  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
419  G_free(indef);
420 
421  return -1;
422  }
423  G_free(indef);
424  }
425  /* input and output CRS are available */
426  else if (info_in->def && info_out->pj && info_out->def) {
427  char *indef = NULL, *outdef = NULL;
428  char *indefcrs = NULL, *outdefcrs = NULL;
429  char *insrid = NULL, *outsrid = NULL;
430  int use_insrid = 0, use_outsrid = 0;
431  PJ *source_crs, *target_crs;
432  PJ_AREA *pj_area = NULL;
433  double xmin, xmax, ymin, ymax;
434  int op_count = 0;
435 
436  /* remove any +towgs84/+nadgrids clause, see above
437  * does not always remove +towgs84=0.000,0.000,0.000 ??? */
438  G_debug(1, "source proj string: %s", info_in->def);
439  G_debug(1, "source type: %s", get_pj_type_string(info_in->pj));
440  indefcrs = info_in->def;
441  source_crs = proj_get_source_crs(NULL, info_in->pj);
442  if (source_crs) {
443  const char *projstr;
444 
445  projstr = proj_as_proj_string(NULL, source_crs, PJ_PROJ_5, NULL);
446  if (projstr) {
447  indefcrs = G_store(projstr);
448  G_debug(1, "Input CRS definition converted from '%s' to '%s'",
449  info_in->def, indefcrs);
450  }
451  proj_destroy(source_crs);
452  source_crs = NULL;
453  }
454 
455  G_debug(1, "target proj string: %s", info_out->def);
456  G_debug(1, "target type: %s", get_pj_type_string(info_out->pj));
457  outdefcrs = info_out->def;
458  target_crs = proj_get_source_crs(NULL, info_out->pj);
459  if (target_crs) {
460  const char *projstr;
461 
462  projstr = proj_as_proj_string(NULL, target_crs, PJ_PROJ_5, NULL);
463  if (projstr) {
464  outdefcrs = G_store(projstr);
465  G_debug(1, "Output CRS definition converted from '%s' to '%s'",
466  info_out->def, outdefcrs);
467  }
468  proj_destroy(target_crs);
469  target_crs = NULL;
470  }
471 
472  /* PROJ6+: EPSG must be uppercase EPSG */
473  if (info_in->srid) {
474  if (strncmp(info_in->srid, "epsg", 4) == 0)
475  insrid = G_store_upper(info_in->srid);
476  else
477  insrid = G_store(info_in->srid);
478  }
479 
480  if (info_out->srid) {
481  if (strncmp(info_out->srid, "epsg", 4) == 0)
482  outsrid = G_store_upper(info_out->srid);
483  else
484  outsrid = G_store(info_out->srid);
485  }
486 
487  if (insrid) {
488  G_asprintf(&indef, "%s", insrid);
489  use_insrid = 1;
490  }
491  else {
492  G_asprintf(&indef, "%s", indefcrs);
493  }
494  G_debug(1, "Input CRS definition: %s", indef);
495 
496  if (outsrid) {
497  G_asprintf(&outdef, "%s", outsrid);
498  use_outsrid = 1;
499  }
500  else {
501  G_asprintf(&outdef, "%s", outdefcrs);
502  }
503  G_debug(1, "Output CRS definition: %s", outdef);
504 
505  /* check number of operations */
506  source_crs = proj_create(NULL, indef);
507  target_crs = proj_create(NULL, outdef);
508 
509  /* get pj_area */
510  if (get_pj_area(&xmin, &xmax, &ymin, &ymax)) {
511  pj_area = proj_area_create();
512  proj_area_set_bbox(pj_area, xmin, ymin, xmax, ymax);
513  }
514 
515  if (source_crs && target_crs) {
516  PJ_OPERATION_FACTORY_CONTEXT *operation_ctx;
517  PJ_OBJ_LIST *op_list;
518 
519  operation_ctx = proj_create_operation_factory_context(NULL, NULL);
520  /* proj_create_operations() works only if both source_crs
521  * and target_crs are found in the proj db
522  * if any is not found, proj can not get a list of operations
523  * and we have to take care of datumshift manually */
524  /* constrain by area ? */
525  op_list = proj_create_operations(NULL,
526  source_crs,
527  target_crs,
528  operation_ctx);
529 
530  op_count = 0;
531  if (op_list)
532  op_count = proj_list_get_count(op_list);
533  if (op_count > 1) {
534  int i;
535 
536  G_warning(_("Found %d possible transformations"), op_count);
537  for (i = 0; i < op_count; i++) {
538  const char *str;
539  const char *area_of_use, *projstr;
540  double e, w, s, n;
541  PJ_PROJ_INFO pj_info;
542  PJ *op, *op_norm;
543 
544  op = proj_list_get(NULL, op_list, i);
545  op_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, op);
546 
547  if (!op_norm) {
548  G_warning(_("proj_normalize_for_visualization() failed for operation %d"),
549  i + 1);
550  }
551  else {
552  proj_destroy(op);
553  op = op_norm;
554  }
555 
556  projstr = proj_as_proj_string(NULL, op,
557  PJ_PROJ_5, NULL);
558  pj_info = proj_pj_info(op);
559  proj_get_area_of_use(NULL, op, &w, &s, &e, &n, &area_of_use);
560  if (projstr) {
561  G_important_message("************************");
562  G_important_message(_("Operation %d:"), i + 1);
563  G_important_message(_("Description: %s"),
564  pj_info.description);
565  G_important_message(" ");
566  G_important_message(_("Area of use: %s"),
567  area_of_use);
568  if (pj_info.accuracy > 0) {
569  G_important_message(" ");
570  G_important_message(_("Accuracy within area of use: %g m"),
571  pj_info.accuracy);
572  }
573 #if PROJ_VERSION_NUM >= 6020000
574  str = proj_get_remarks(op);
575  if (str && *str) {
576  G_important_message(" ");
577  G_important_message(_("Remarks: %s"), str);
578  }
579  str = proj_get_scope(op);
580  if (str && *str) {
581  G_important_message(" ");
582  G_important_message(_("Scope: %s"), str);
583  }
584 #endif
585  G_important_message(" ");
586  G_important_message(_("PROJ string:"));
587  G_important_message("%s", projstr);
588  }
589  proj_destroy(op);
590  }
591  G_important_message("************************");
592 
593  G_important_message(_("See also output of:"));
594  G_important_message("projinfo -o PROJ -s \"%s\" -t \"%s\"",
595  indef, outdef);
596  G_important_message(_("Please provide the appropriate PROJ string with the %s option"),
597  "pipeline");
598  G_important_message("************************");
599  }
600 
601  if (op_list)
602  proj_list_destroy(op_list);
603  proj_operation_factory_context_destroy(operation_ctx);
604  }
605  if (source_crs)
606  proj_destroy(source_crs);
607  if (target_crs)
608  proj_destroy(target_crs);
609 
610  /* try proj_create_crs_to_crs() */
611  G_debug(1, "trying %s to %s", indef, outdef);
612 
613  info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
614  indef,
615  outdef,
616  pj_area);
617 
618  if (info_trans->pj) {
619  const char *projstr = NULL;
620 
621  projstr = proj_as_proj_string(NULL, info_trans->pj,
622  PJ_PROJ_5, NULL);
623  if (projstr == NULL) {
624  G_debug(1, "proj_create_crs_to_crs() failed with PROJ%d for input \"%s\", output \"%s\"",
625  PROJ_VERSION_MAJOR, indef, outdef);
626 
627  G_asprintf(&indef, "%s", indefcrs);
628  G_asprintf(&outdef, "%s", outdefcrs);
629 
630  G_debug(1, "trying %s to %s", indef, outdef);
631 
632  /* try proj_create_crs_to_crs() */
633  proj_destroy(info_trans->pj);
634  info_trans->pj = NULL;
635  info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
636  indef,
637  outdef,
638  NULL);
639  }
640  else {
641  /* PROJ will do the unit conversion if set up from srid
642  * -> disable unit conversion for GPJ_transform */
643  /* ugly hack */
644  if (use_insrid) {
645  ((struct pj_info *)info_in)->meters = 1;
646  }
647  if (use_outsrid) {
648  ((struct pj_info *)info_out)->meters = 1;
649  }
650  }
651  }
652 
653  if (info_trans->pj) {
654  const char *projstr;
655  PJ *pj_norm = NULL;
656 
657  G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ%d",
658  PROJ_VERSION_MAJOR);
659 
660  projstr = proj_as_proj_string(NULL, info_trans->pj,
661  PJ_PROJ_5, NULL);
662 
663  info_trans->def = G_store(projstr);
664 
665  if (projstr) {
666  /* make sure axis order is easting, northing
667  * proj_normalize_for_visualization() requires
668  * source and target CRS
669  * -> does not work with ll equivalent of input:
670  * no target CRS in +proj=pipeline +step +inv %s */
671  pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, info_trans->pj);
672 
673  if (!pj_norm) {
674  G_warning(_("proj_normalize_for_visualization() failed for '%s'"),
675  info_trans->def);
676  }
677  else {
678  proj_destroy(info_trans->pj);
679  info_trans->pj = pj_norm;
680  projstr = proj_as_proj_string(NULL, info_trans->pj,
681  PJ_PROJ_5, NULL);
682  info_trans->def = G_store(projstr);
683  }
684  if (op_count > 1) {
685  G_important_message(_("Selected pipeline:"));
686  G_important_message(_("%s"), info_trans->def);
687  G_important_message("************************");
688  }
689 
690  G_debug(1, "proj_create_crs_to_crs() pipeline: %s", info_trans->def);
691  }
692  else {
693  proj_destroy(info_trans->pj);
694  info_trans->pj = NULL;
695  }
696  }
697  /* last try with proj_create() */
698  if (info_trans->pj == NULL) {
699  G_debug(1, "proj_create_crs_to_crs() failed with PROJ%d for input \"%s\", output \"%s\"",
700  PROJ_VERSION_MAJOR, indef, outdef);
701 
702  G_warning("GPJ_init_transform(): falling back to proj_create()");
703 
704  if (insrid) {
705  G_free(indef);
706  }
707  G_asprintf(&indef, "%s", info_in->def);
708  if (outsrid) {
709  G_free(outdef);
710  }
711  G_asprintf(&outdef, "%s", info_out->def);
712  /* try proj_create() with +proj=pipeline +step +inv %s +step %s" */
713  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s +step %s",
714  indef, outdef);
715 
716  info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
717  }
718 
719  if (pj_area)
720  proj_area_destroy(pj_area);
721 
722  if (insrid)
723  G_free(insrid);
724  if (outsrid)
725  G_free(outsrid);
726  G_free(indef);
727  G_free(outdef);
728  }
729  if (info_trans->pj == NULL) {
730  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
731 
732  return -1;
733  }
734 #else /* PROJ 5 */
735  info_trans->pj = NULL;
736  info_trans->meters = 1.;
737  info_trans->zone = 0;
738  sprintf(info_trans->proj, "pipeline");
739 
740  /* user-provided pipeline */
741  if (info_trans->def) {
742  /* create a pj from user-defined transformation pipeline */
743  info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
744  if (info_trans->pj == NULL) {
745  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
746 
747  return -1;
748  }
749  }
750  /* if no output CRS is defined,
751  * assume info_out to be ll equivalent of info_in */
752  else if (info_out->pj == NULL) {
753  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",
754  info_in->def);
755  info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
756  if (info_trans->pj == NULL) {
757  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
758 
759  return -1;
760  }
761  }
762  else if (info_in->def && info_out->pj && info_out->def) {
763  char *indef = NULL, *outdef = NULL;
764  char *insrid = NULL, *outsrid = NULL;
765 
766  /* PROJ5: EPSG must be lowercase epsg */
767  if (info_in->srid) {
768  if (strncmp(info_in->srid, "EPSG", 4) == 0)
769  insrid = G_store_lower(info_in->srid);
770  else
771  insrid = G_store(info_in->srid);
772  }
773 
774  if (info_out->pj && info_out->srid) {
775  if (strncmp(info_out->srid, "EPSG", 4) == 0)
776  outsrid = G_store_lower(info_out->srid);
777  else
778  outsrid = G_store(info_out->srid);
779  }
780 
781  info_trans->pj = NULL;
782 
783  if (insrid && outsrid) {
784  G_asprintf(&indef, "%s", insrid);
785  G_asprintf(&outdef, "%s", outsrid);
786 
787  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv +init=%s +step +init=%s",
788  indef, outdef);
789 
790  /* try proj_create_crs_to_crs() */
791  info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
792  indef,
793  outdef,
794  NULL);
795  }
796 
797  if (info_trans->pj) {
798  G_debug(1, "proj_create_crs_to_crs() succeeded with PROJ5");
799  }
800  else {
801  if (indef) {
802  G_free(indef);
803  indef = NULL;
804  }
805  if (insrid) {
806  G_asprintf(&indef, "+init=%s", insrid);
807  }
808  else {
809  G_asprintf(&indef, "%s", info_in->def);
810  }
811 
812  if (outdef) {
813  G_free(outdef);
814  outdef = NULL;
815  }
816  if (outsrid) {
817  G_asprintf(&outdef, "+init=%s", outsrid);
818  }
819  else {
820  G_asprintf(&outdef, "%s", info_out->def);
821  }
822 
823  /* try proj_create() with +proj=pipeline +step +inv %s +step %s" */
824  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s +step %s",
825  indef, outdef);
826 
827  info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
828  }
829  if (insrid)
830  G_free(insrid);
831  if (outsrid)
832  G_free(outsrid);
833  G_free(indef);
834  G_free(outdef);
835  }
836  if (info_trans->pj == NULL) {
837  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
838 
839  return -1;
840  }
841 
842 #endif
843 #else /* PROJ 4 */
844  if (info_out->pj == NULL) {
845  if (GPJ_get_equivalent_latlong(info_out, info_in) < 0) {
846  G_warning(_("Unable to create latlong equivalent for '%s'"),
847  info_in->def);
848 
849  return -1;
850  }
851  }
852 #endif
853 
854  return 1;
855 }
856 
857 /* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
858 
859 /**
860  * \brief Re-project a point between two co-ordinate systems using a
861  * transformation object prepared with GPJ_prepare_pj()
862  *
863  * This function takes pointers to three pj_info structures as arguments,
864  * and projects a point between the input and output co-ordinate system.
865  * The pj_info structure info_trans must have been initialized with
866  * GPJ_init_transform().
867  * The direction determines if a point is projected from input CRS to
868  * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
869  * The easting, northing, and height of the point are contained in the
870  * pointers passed to the function; these will be overwritten by the
871  * coordinates of the transformed point.
872  *
873  * \param info_in pointer to pj_info struct for input co-ordinate system
874  * \param info_out pointer to pj_info struct for output co-ordinate system
875  * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
876  * \param dir direction of the transformation (PJ_FWD or PJ_INV)
877  * \param x Pointer to a double containing easting or longitude
878  * \param y Pointer to a double containing northing or latitude
879  * \param z Pointer to a double containing height, or NULL
880  *
881  * \return Return value from PROJ proj_trans() function
882  **/
883 
884 int GPJ_transform(const struct pj_info *info_in,
885  const struct pj_info *info_out,
886  const struct pj_info *info_trans, int dir,
887  double *x, double *y, double *z)
888 {
889  int ok = 0;
890 
891 #ifdef HAVE_PROJ_H
892  /* PROJ 5+ variant */
893  int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
894  PJ_COORD c;
895 
896  if (info_in->pj == NULL)
897  G_fatal_error(_("No input projection"));
898 
899  if (info_trans->pj == NULL)
900  G_fatal_error(_("No transformation object"));
901 
902  in_deg2rad = out_rad2deg = 1;
903  if (dir == PJ_FWD) {
904  /* info_in -> info_out */
905  METERS_in = info_in->meters;
906  in_is_ll = !strncmp(info_in->proj, "ll", 2);
907 #if PROJ_VERSION_MAJOR >= 6
908  /* PROJ 6+: conversion to radians is not always needed:
909  * if proj_angular_input(info_trans->pj, dir) == 1
910  * -> convert from degrees to radians */
911  if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
912  in_deg2rad = 0;
913  }
914 #endif
915  if (info_out->pj) {
916  METERS_out = info_out->meters;
917  out_is_ll = !strncmp(info_out->proj, "ll", 2);
918 #if PROJ_VERSION_MAJOR >= 6
919  /* PROJ 6+: conversion to radians is not always needed:
920  * if proj_angular_input(info_trans->pj, dir) == 1
921  * -> convert from degrees to radians */
922  if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
923  out_rad2deg = 0;
924  }
925 #endif
926  }
927  else {
928  METERS_out = 1.0;
929  out_is_ll = 1;
930  }
931  }
932  else {
933  /* info_out -> info_in */
934  METERS_out = info_in->meters;
935  out_is_ll = !strncmp(info_in->proj, "ll", 2);
936 #if PROJ_VERSION_MAJOR >= 6
937  /* PROJ 6+: conversion to radians is not always needed:
938  * if proj_angular_input(info_trans->pj, dir) == 1
939  * -> convert from degrees to radians */
940  if (out_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
941  out_rad2deg = 0;
942  }
943 #endif
944  if (info_out->pj) {
945  METERS_in = info_out->meters;
946  in_is_ll = !strncmp(info_out->proj, "ll", 2);
947 #if PROJ_VERSION_MAJOR >= 6
948  /* PROJ 6+: conversion to radians is not always needed:
949  * if proj_angular_input(info_trans->pj, dir) == 1
950  * -> convert from degrees to radians */
951  if (in_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
952  in_deg2rad = 0;
953  }
954 #endif
955  }
956  else {
957  METERS_in = 1.0;
958  in_is_ll = 1;
959  }
960  }
961 
962  /* prepare */
963  if (in_is_ll) {
964  if (in_deg2rad) {
965  /* convert degrees to radians */
966  c.lpzt.lam = (*x) / RAD_TO_DEG;
967  c.lpzt.phi = (*y) / RAD_TO_DEG;
968  }
969  else {
970  c.lpzt.lam = (*x);
971  c.lpzt.phi = (*y);
972  }
973  c.lpzt.z = 0;
974  if (z)
975  c.lpzt.z = *z;
976  c.lpzt.t = 0;
977  }
978  else {
979  /* convert to meters */
980  c.xyzt.x = *x * METERS_in;
981  c.xyzt.y = *y * METERS_in;
982  c.xyzt.z = 0;
983  if (z)
984  c.xyzt.z = *z;
985  c.xyzt.t = 0;
986  }
987 
988  /* transform */
989  c = proj_trans(info_trans->pj, dir, c);
990  ok = proj_errno(info_trans->pj);
991 
992  if (ok < 0) {
993  G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
994  return ok;
995  }
996 
997  /* output */
998  if (out_is_ll) {
999  /* convert to degrees */
1000  if (out_rad2deg) {
1001  /* convert radians to degrees */
1002  *x = c.lp.lam * RAD_TO_DEG;
1003  *y = c.lp.phi * RAD_TO_DEG;
1004  }
1005  else {
1006  *x = c.lp.lam;
1007  *y = c.lp.phi;
1008  }
1009  if (z)
1010  *z = c.lpzt.z;
1011  }
1012  else {
1013  /* convert to map units */
1014  *x = c.xyzt.x / METERS_out;
1015  *y = c.xyzt.y / METERS_out;
1016  if (z)
1017  *z = c.xyzt.z;
1018  }
1019 #else
1020  /* PROJ 4 variant */
1021  double u, v;
1022  double h = 0.0;
1023  const struct pj_info *p_in, *p_out;
1024 
1025  if (info_out == NULL)
1026  G_fatal_error(_("No output projection"));
1027 
1028  if (dir == PJ_FWD) {
1029  p_in = info_in;
1030  p_out = info_out;
1031  }
1032  else {
1033  p_in = info_out;
1034  p_out = info_in;
1035  }
1036 
1037  METERS_in = p_in->meters;
1038  METERS_out = p_out->meters;
1039 
1040  if (z)
1041  h = *z;
1042 
1043  if (strncmp(p_in->proj, "ll", 2) == 0) {
1044  u = (*x) / RAD_TO_DEG;
1045  v = (*y) / RAD_TO_DEG;
1046  }
1047  else {
1048  u = *x * METERS_in;
1049  v = *y * METERS_in;
1050  }
1051 
1052  ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
1053 
1054  if (ok < 0) {
1055  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1056  return ok;
1057  }
1058 
1059  if (strncmp(p_out->proj, "ll", 2) == 0) {
1060  *x = u * RAD_TO_DEG;
1061  *y = v * RAD_TO_DEG;
1062  }
1063  else {
1064  *x = u / METERS_out;
1065  *y = v / METERS_out;
1066  }
1067  if (z)
1068  *z = h;
1069 #endif
1070 
1071  return ok;
1072 }
1073 
1074 /**
1075  * \brief Re-project an array of points between two co-ordinate systems
1076  * using a transformation object prepared with GPJ_prepare_pj()
1077  *
1078  * This function takes pointers to three pj_info structures as arguments,
1079  * and projects an array of pointd between the input and output
1080  * co-ordinate system. The pj_info structure info_trans must have been
1081  * initialized with GPJ_init_transform().
1082  * The direction determines if a point is projected from input CRS to
1083  * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
1084  * The easting, northing, and height of the point are contained in the
1085  * pointers passed to the function; these will be overwritten by the
1086  * coordinates of the transformed point.
1087  *
1088  * \param info_in pointer to pj_info struct for input co-ordinate system
1089  * \param info_out pointer to pj_info struct for output co-ordinate system
1090  * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
1091  * \param dir direction of the transformation (PJ_FWD or PJ_INV)
1092  * \param x pointer to an array of type double containing easting or longitude
1093  * \param y pointer to an array of type double containing northing or latitude
1094  * \param z pointer to an array of type double containing height, or NULL
1095  * \param n number of points in the arrays to be transformed
1096  *
1097  * \return Return value from PROJ proj_trans() function
1098  **/
1099 
1100 int GPJ_transform_array(const struct pj_info *info_in,
1101  const struct pj_info *info_out,
1102  const struct pj_info *info_trans, int dir,
1103  double *x, double *y, double *z, int n)
1104 {
1105  int ok;
1106  int i;
1107  int has_z = 1;
1108 
1109 #ifdef HAVE_PROJ_H
1110  /* PROJ 5+ variant */
1111  int in_is_ll, out_is_ll, in_deg2rad, out_rad2deg;
1112  PJ_COORD c;
1113 
1114  if (info_trans->pj == NULL)
1115  G_fatal_error(_("No transformation object"));
1116 
1117  in_deg2rad = out_rad2deg = 1;
1118  if (dir == PJ_FWD) {
1119  /* info_in -> info_out */
1120  METERS_in = info_in->meters;
1121  in_is_ll = !strncmp(info_in->proj, "ll", 2);
1122 #if PROJ_VERSION_MAJOR >= 6
1123  /* PROJ 6+: conversion to radians is not always needed:
1124  * if proj_angular_input(info_trans->pj, dir) == 1
1125  * -> convert from degrees to radians */
1126  if (in_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1127  in_deg2rad = 0;
1128  }
1129 #endif
1130  if (info_out->pj) {
1131  METERS_out = info_out->meters;
1132  out_is_ll = !strncmp(info_out->proj, "ll", 2);
1133 #if PROJ_VERSION_MAJOR >= 6
1134  /* PROJ 6+: conversion to radians is not always needed:
1135  * if proj_angular_input(info_trans->pj, dir) == 1
1136  * -> convert from degrees to radians */
1137  if (out_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1138  out_rad2deg = 0;
1139  }
1140 #endif
1141  }
1142  else {
1143  METERS_out = 1.0;
1144  out_is_ll = 1;
1145  }
1146  }
1147  else {
1148  /* info_out -> info_in */
1149  METERS_out = info_in->meters;
1150  out_is_ll = !strncmp(info_in->proj, "ll", 2);
1151 #if PROJ_VERSION_MAJOR >= 6
1152  /* PROJ 6+: conversion to radians is not always needed:
1153  * if proj_angular_input(info_trans->pj, dir) == 1
1154  * -> convert from degrees to radians */
1155  if (out_is_ll && proj_angular_input(info_trans->pj, dir) == 0) {
1156  out_rad2deg = 0;
1157  }
1158 #endif
1159  if (info_out->pj) {
1160  METERS_in = info_out->meters;
1161  in_is_ll = !strncmp(info_out->proj, "ll", 2);
1162 #if PROJ_VERSION_MAJOR >= 6
1163  /* PROJ 6+: conversion to degrees is not always needed:
1164  * if proj_angular_output(info_trans->pj, dir) == 1
1165  * -> convert from degrees to radians */
1166  if (in_is_ll && proj_angular_output(info_trans->pj, dir) == 0) {
1167  in_deg2rad = 0;
1168  }
1169 #endif
1170  }
1171  else {
1172  METERS_in = 1.0;
1173  in_is_ll = 1;
1174  }
1175  }
1176 
1177  if (z == NULL) {
1178  z = G_malloc(sizeof(double) * n);
1179  /* they say memset is only guaranteed for chars ;-( */
1180  for (i = 0; i < n; i++)
1181  z[i] = 0.0;
1182  has_z = 0;
1183  }
1184  ok = 0;
1185  if (in_is_ll) {
1186  c.lpzt.t = 0;
1187  if (out_is_ll) {
1188  /* what is more costly ?
1189  * calling proj_trans for each point
1190  * or having three loops over all points ?
1191  * proj_trans_array() itself calls proj_trans() in a loop
1192  * -> one loop over all points is better than
1193  * three loops over all points
1194  */
1195  for (i = 0; i < n; i++) {
1196  if (in_deg2rad) {
1197  /* convert degrees to radians */
1198  c.lpzt.lam = (*x) / RAD_TO_DEG;
1199  c.lpzt.phi = (*y) / RAD_TO_DEG;
1200  }
1201  else {
1202  c.lpzt.lam = (*x);
1203  c.lpzt.phi = (*y);
1204  }
1205  c.lpzt.z = z[i];
1206  c = proj_trans(info_trans->pj, dir, c);
1207  if ((ok = proj_errno(info_trans->pj)) < 0)
1208  break;
1209  if (out_rad2deg) {
1210  /* convert radians to degrees */
1211  x[i] = c.lp.lam * RAD_TO_DEG;
1212  y[i] = c.lp.phi * RAD_TO_DEG;
1213  }
1214  else {
1215  x[i] = c.lp.lam;
1216  y[i] = c.lp.phi;
1217  }
1218  }
1219  }
1220  else {
1221  for (i = 0; i < n; i++) {
1222  if (in_deg2rad) {
1223  /* convert degrees to radians */
1224  c.lpzt.lam = (*x) / RAD_TO_DEG;
1225  c.lpzt.phi = (*y) / RAD_TO_DEG;
1226  }
1227  else {
1228  c.lpzt.lam = (*x);
1229  c.lpzt.phi = (*y);
1230  }
1231  c.lpzt.z = z[i];
1232  c = proj_trans(info_trans->pj, dir, c);
1233  if ((ok = proj_errno(info_trans->pj)) < 0)
1234  break;
1235 
1236  /* convert to map units */
1237  x[i] = c.xy.x / METERS_out;
1238  y[i] = c.xy.y / METERS_out;
1239  }
1240  }
1241  }
1242  else {
1243  c.xyzt.t = 0;
1244  if (out_is_ll) {
1245  for (i = 0; i < n; i++) {
1246  /* convert to meters */
1247  c.xyzt.x = x[i] * METERS_in;
1248  c.xyzt.y = y[i] * METERS_in;
1249  c.xyzt.z = z[i];
1250  c = proj_trans(info_trans->pj, dir, c);
1251  if ((ok = proj_errno(info_trans->pj)) < 0)
1252  break;
1253  if (out_rad2deg) {
1254  /* convert radians to degrees */
1255  x[i] = c.lp.lam * RAD_TO_DEG;
1256  y[i] = c.lp.phi * RAD_TO_DEG;
1257  }
1258  else {
1259  x[i] = c.lp.lam;
1260  y[i] = c.lp.phi;
1261  }
1262  }
1263  }
1264  else {
1265  for (i = 0; i < n; i++) {
1266  /* convert to meters */
1267  c.xyzt.x = x[i] * METERS_in;
1268  c.xyzt.y = y[i] * METERS_in;
1269  c.xyzt.z = z[i];
1270  c = proj_trans(info_trans->pj, dir, c);
1271  if ((ok = proj_errno(info_trans->pj)) < 0)
1272  break;
1273  /* convert to map units */
1274  x[i] = c.xy.x / METERS_out;
1275  y[i] = c.xy.y / METERS_out;
1276  }
1277  }
1278  }
1279  if (!has_z)
1280  G_free(z);
1281 
1282  if (ok < 0) {
1283  G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
1284  }
1285 #else
1286  /* PROJ 4 variant */
1287  const struct pj_info *p_in, *p_out;
1288 
1289  if (dir == PJ_FWD) {
1290  p_in = info_in;
1291  p_out = info_out;
1292  }
1293  else {
1294  p_in = info_out;
1295  p_out = info_in;
1296  }
1297 
1298  METERS_in = p_in->meters;
1299  METERS_out = p_out->meters;
1300 
1301  if (z == NULL) {
1302  z = G_malloc(sizeof(double) * n);
1303  /* they say memset is only guaranteed for chars ;-( */
1304  for (i = 0; i < n; ++i)
1305  z[i] = 0.0;
1306  has_z = 0;
1307  }
1308  if (strncmp(p_in->proj, "ll", 2) == 0) {
1309  if (strncmp(p_out->proj, "ll", 2) == 0) {
1310  DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
1311  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1312  MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
1313  }
1314  else {
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  DIVIDE_LOOP(x, y, n, METERS_out);
1318  }
1319  }
1320  else {
1321  if (strncmp(p_out->proj, "ll", 2) == 0) {
1322  MULTIPLY_LOOP(x, y, n, METERS_in);
1323  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1324  MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
1325  }
1326  else {
1327  MULTIPLY_LOOP(x, y, n, METERS_in);
1328  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
1329  DIVIDE_LOOP(x, y, n, METERS_out);
1330  }
1331  }
1332  if (!has_z)
1333  G_free(z);
1334 
1335  if (ok < 0)
1336  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1337 #endif
1338 
1339  return ok;
1340 }
1341 
1342 /*
1343  * old API, to be deleted
1344  */
1345 
1346 /**
1347  * \brief Re-project a point between two co-ordinate systems
1348  *
1349  * This function takes pointers to two pj_info structures as arguments,
1350  * and projects a point between the co-ordinate systems represented by them.
1351  * The easting and northing of the point are contained in two pointers passed
1352  * to the function; these will be overwritten by the co-ordinates of the
1353  * re-projected point.
1354  *
1355  * \param x Pointer to a double containing easting or longitude
1356  * \param y Pointer to a double containing northing or latitude
1357  * \param info_in pointer to pj_info struct for input co-ordinate system
1358  * \param info_out pointer to pj_info struct for output co-ordinate system
1359  *
1360  * \return Return value from PROJ proj_trans() function
1361  **/
1362 
1363 int pj_do_proj(double *x, double *y,
1364  const struct pj_info *info_in, const struct pj_info *info_out)
1365 {
1366  int ok;
1367 #ifdef HAVE_PROJ_H
1368  struct pj_info info_trans;
1369  PJ_COORD c;
1370 
1371  if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
1372  return -1;
1373  }
1374 
1375  METERS_in = info_in->meters;
1376  METERS_out = info_out->meters;
1377 
1378  if (strncmp(info_in->proj, "ll", 2) == 0) {
1379  /* convert to radians */
1380  c.lpzt.lam = (*x) / RAD_TO_DEG;
1381  c.lpzt.phi = (*y) / RAD_TO_DEG;
1382  c.lpzt.z = 0;
1383  c.lpzt.t = 0;
1384  c = proj_trans(info_trans.pj, PJ_FWD, c);
1385  ok = proj_errno(info_trans.pj);
1386 
1387  if (strncmp(info_out->proj, "ll", 2) == 0) {
1388  /* convert to degrees */
1389  *x = c.lp.lam * RAD_TO_DEG;
1390  *y = c.lp.phi * RAD_TO_DEG;
1391  }
1392  else {
1393  /* convert to map units */
1394  *x = c.xy.x / METERS_out;
1395  *y = c.xy.y / METERS_out;
1396  }
1397  }
1398  else {
1399  /* convert to meters */
1400  c.xyzt.x = *x * METERS_in;
1401  c.xyzt.y = *y * METERS_in;
1402  c.xyzt.z = 0;
1403  c.xyzt.t = 0;
1404  c = proj_trans(info_trans.pj, PJ_FWD, c);
1405  ok = proj_errno(info_trans.pj);
1406 
1407  if (strncmp(info_out->proj, "ll", 2) == 0) {
1408  /* convert to degrees */
1409  *x = c.lp.lam * RAD_TO_DEG;
1410  *y = c.lp.phi * RAD_TO_DEG;
1411  }
1412  else {
1413  /* convert to map units */
1414  *x = c.xy.x / METERS_out;
1415  *y = c.xy.y / METERS_out;
1416  }
1417  }
1418  proj_destroy(info_trans.pj);
1419 
1420  if (ok < 0) {
1421  G_warning(_("proj_trans() failed: %d"), ok);
1422  }
1423 #else
1424  double u, v;
1425  double h = 0.0;
1426 
1427  METERS_in = info_in->meters;
1428  METERS_out = info_out->meters;
1429 
1430  if (strncmp(info_in->proj, "ll", 2) == 0) {
1431  if (strncmp(info_out->proj, "ll", 2) == 0) {
1432  u = (*x) / RAD_TO_DEG;
1433  v = (*y) / RAD_TO_DEG;
1434  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1435  *x = u * RAD_TO_DEG;
1436  *y = v * RAD_TO_DEG;
1437  }
1438  else {
1439  u = (*x) / RAD_TO_DEG;
1440  v = (*y) / RAD_TO_DEG;
1441  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1442  *x = u / METERS_out;
1443  *y = v / METERS_out;
1444  }
1445  }
1446  else {
1447  if (strncmp(info_out->proj, "ll", 2) == 0) {
1448  u = *x * METERS_in;
1449  v = *y * METERS_in;
1450  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1451  *x = u * RAD_TO_DEG;
1452  *y = v * RAD_TO_DEG;
1453  }
1454  else {
1455  u = *x * METERS_in;
1456  v = *y * METERS_in;
1457  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
1458  *x = u / METERS_out;
1459  *y = v / METERS_out;
1460  }
1461  }
1462  if (ok < 0) {
1463  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1464  }
1465 #endif
1466  return ok;
1467 }
1468 
1469 /**
1470  * \brief Re-project an array of points between two co-ordinate systems with
1471  * optional ellipsoidal height conversion
1472  *
1473  * This function takes pointers to two pj_info structures as arguments,
1474  * and projects an array of points between the co-ordinate systems
1475  * represented by them. Pointers to the three arrays of easting, northing,
1476  * and ellipsoidal height of the point (this one may be NULL) are passed
1477  * to the function; these will be overwritten by the co-ordinates of the
1478  * re-projected points.
1479  *
1480  * \param count Number of points in the arrays to be transformed
1481  * \param x Pointer to an array of type double containing easting or longitude
1482  * \param y Pointer to an array of type double containing northing or latitude
1483  * \param h Pointer to an array of type double containing ellipsoidal height.
1484  * May be null in which case a two-dimensional re-projection will be
1485  * done
1486  * \param info_in pointer to pj_info struct for input co-ordinate system
1487  * \param info_out pointer to pj_info struct for output co-ordinate system
1488  *
1489  * \return Return value from PROJ proj_trans() function
1490  **/
1491 
1492 int pj_do_transform(int count, double *x, double *y, double *h,
1493  const struct pj_info *info_in, const struct pj_info *info_out)
1494 {
1495  int ok;
1496  int i;
1497  int has_h = 1;
1498 #ifdef HAVE_PROJ_H
1499  struct pj_info info_trans;
1500  PJ_COORD c;
1501 
1502  if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
1503  return -1;
1504  }
1505 
1506  METERS_in = info_in->meters;
1507  METERS_out = info_out->meters;
1508 
1509  if (h == NULL) {
1510  h = G_malloc(sizeof *h * count);
1511  /* they say memset is only guaranteed for chars ;-( */
1512  for (i = 0; i < count; ++i)
1513  h[i] = 0.0;
1514  has_h = 0;
1515  }
1516  ok = 0;
1517  if (strncmp(info_in->proj, "ll", 2) == 0) {
1518  c.lpzt.t = 0;
1519  if (strncmp(info_out->proj, "ll", 2) == 0) {
1520  for (i = 0; i < count; i++) {
1521  /* convert to radians */
1522  c.lpzt.lam = x[i] / RAD_TO_DEG;
1523  c.lpzt.phi = y[i] / RAD_TO_DEG;
1524  c.lpzt.z = h[i];
1525  c = proj_trans(info_trans.pj, PJ_FWD, c);
1526  if ((ok = proj_errno(info_trans.pj)) < 0)
1527  break;
1528  /* convert to degrees */
1529  x[i] = c.lp.lam * RAD_TO_DEG;
1530  y[i] = c.lp.phi * RAD_TO_DEG;
1531  }
1532  }
1533  else {
1534  for (i = 0; i < count; i++) {
1535  /* convert to radians */
1536  c.lpzt.lam = x[i] / RAD_TO_DEG;
1537  c.lpzt.phi = y[i] / RAD_TO_DEG;
1538  c.lpzt.z = h[i];
1539  c = proj_trans(info_trans.pj, PJ_FWD, c);
1540  if ((ok = proj_errno(info_trans.pj)) < 0)
1541  break;
1542  /* convert to map units */
1543  x[i] = c.xy.x / METERS_out;
1544  y[i] = c.xy.y / METERS_out;
1545  }
1546  }
1547  }
1548  else {
1549  c.xyzt.t = 0;
1550  if (strncmp(info_out->proj, "ll", 2) == 0) {
1551  for (i = 0; i < count; i++) {
1552  /* convert to meters */
1553  c.xyzt.x = x[i] * METERS_in;
1554  c.xyzt.y = y[i] * METERS_in;
1555  c.xyzt.z = h[i];
1556  c = proj_trans(info_trans.pj, PJ_FWD, c);
1557  if ((ok = proj_errno(info_trans.pj)) < 0)
1558  break;
1559  /* convert to degrees */
1560  x[i] = c.lp.lam * RAD_TO_DEG;
1561  y[i] = c.lp.phi * RAD_TO_DEG;
1562  }
1563  }
1564  else {
1565  for (i = 0; i < count; i++) {
1566  /* convert to meters */
1567  c.xyzt.x = x[i] * METERS_in;
1568  c.xyzt.y = y[i] * METERS_in;
1569  c.xyzt.z = h[i];
1570  c = proj_trans(info_trans.pj, PJ_FWD, c);
1571  if ((ok = proj_errno(info_trans.pj)) < 0)
1572  break;
1573  /* convert to map units */
1574  x[i] = c.xy.x / METERS_out;
1575  y[i] = c.xy.y / METERS_out;
1576  }
1577  }
1578  }
1579  if (!has_h)
1580  G_free(h);
1581  proj_destroy(info_trans.pj);
1582 
1583  if (ok < 0) {
1584  G_warning(_("proj_trans() failed: %d"), ok);
1585  }
1586 #else
1587  METERS_in = info_in->meters;
1588  METERS_out = info_out->meters;
1589 
1590  if (h == NULL) {
1591  h = G_malloc(sizeof *h * count);
1592  /* they say memset is only guaranteed for chars ;-( */
1593  for (i = 0; i < count; ++i)
1594  h[i] = 0.0;
1595  has_h = 0;
1596  }
1597  if (strncmp(info_in->proj, "ll", 2) == 0) {
1598  if (strncmp(info_out->proj, "ll", 2) == 0) {
1599  DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
1600  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1601  MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
1602  }
1603  else {
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  DIVIDE_LOOP(x, y, count, METERS_out);
1607  }
1608  }
1609  else {
1610  if (strncmp(info_out->proj, "ll", 2) == 0) {
1611  MULTIPLY_LOOP(x, y, count, METERS_in);
1612  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1613  MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
1614  }
1615  else {
1616  MULTIPLY_LOOP(x, y, count, METERS_in);
1617  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
1618  DIVIDE_LOOP(x, y, count, METERS_out);
1619  }
1620  }
1621  if (!has_h)
1622  G_free(h);
1623 
1624  if (ok < 0) {
1625  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
1626  }
1627 #endif
1628  return ok;
1629 }
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:1492
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:1363
void G_get_window(struct Cell_head *window)
Get the current region.
Definition: get_window.c:47
void G_important_message(const char *msg,...)
Print a message to stderr even in brief mode (verbosity=1)
Definition: gis/error.c:131
struct Key_Value * G_get_projunits(void)
Gets units information for location.
Definition: get_projinfo.c:30
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
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
int count
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
#define DIVIDE_LOOP(x, y, c, m)
Definition: do_proj.c:42
int G_asprintf(char **out, const char *fmt,...)
Definition: asprintf.c:70
char * G_store_lower(const char *s)
Copy string to allocated memory and convert copied string to lower case.
Definition: strings.c:140
#define NULL
Definition: ccmath.h:32
#define x
void G_unset_window()
Unset current region.
Definition: get_window.c:132
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition: key_value1.c:103
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
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:1100
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
Definition: get_projinfo.c:59
char * G_store_upper(const char *s)
Copy string to allocated memory and convert copied string to upper case.
Definition: strings.c:116
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
#define MULTIPLY_LOOP(x, y, c, m)
Definition: do_proj.c:34
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:884
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204