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