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