GRASS GIS 7 Programmer's Manual  7.8.0(2019)-exported
do_proj.c
Go to the documentation of this file.
1 
2 /**
3  \file do_proj.c
4 
5  \brief GProj library - Functions for re-projecting point data
6 
7  \author Original Author unknown, probably Soil Conservation Service
8  Eric Miller, Paul Kelly, Markus Metz
9 
10  (C) 2003-2008,2018 by the GRASS Development Team
11 
12  This program is free software under the GNU General Public
13  License (>=v2). Read the file COPYING that comes with GRASS
14  for details.
15 **/
16 
17 #include <stdio.h>
18 #include <string.h>
19 #include <ctype.h>
20 
21 #include <grass/gis.h>
22 #include <grass/gprojects.h>
23 #include <grass/glocale.h>
24 
25 /* 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 /**
45  * \brief Create a PROJ transformation object to transform coordinates
46  * from an input SRS to an output SRS
47  *
48  * After the transformation has been initialized with this function,
49  * coordinates can be transformed from input SRS to output SRS with
50  * GPJ_transform() and direction = PJ_FWD, and back from output SRS to
51  * input SRS with direction = OJ_INV.
52  * If coordinates should be transformed between the input SRS and its
53  * latlong equivalent, an uninitialized info_out with
54  * info_out->pj = NULL can be passed to the function. In this case,
55  * coordinates will be transformed between the input SRS and its
56  * latlong equivalent, and for PROJ 5+, the transformation object is
57  * created accordingly, while for PROJ 4, the output SRS is created as
58  * latlong equivalent of the input SRS
59  *
60 * PROJ 5+:
61  * info_in->pj must not be null
62  * if info_out->pj is null, assume info_out to be the ll equivalent
63  * of info_in
64  * create info_trans as conversion from info_in to its ll equivalent
65  * NOTE: this is the inverse of the logic of PROJ 5 which by default
66  * converts from ll to a given SRS, not from a given SRS to ll
67  * thus PROJ 5+ itself uses an inverse transformation in the
68  * first step of the pipeline for proj_create_crs_to_crs()
69  * if info_trans->def is not NULL, this pipeline definition will be
70  * used to create a transformation object
71  * PROJ 4:
72  * info_in->pj must not be null
73  * if info_out->pj is null, create info_out as ll equivalent
74  * else do nothing, info_trans is not used
75  *
76  * \param info_in pointer to pj_info struct for input co-ordinate system
77  * \param info_out pointer to pj_info struct for output co-ordinate system
78  * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
79  *
80  * \return 1 on success, -1 on failure
81  **/
82 int GPJ_init_transform(const struct pj_info *info_in,
83  const struct pj_info *info_out,
84  struct pj_info *info_trans)
85 {
86  if (info_in->pj == NULL)
87  G_fatal_error(_("Input coordinate system is NULL"));
88 
89 #ifdef HAVE_PROJ_H
90  info_trans->pj = NULL;
91  if (!info_trans->def) {
92  if (info_in->srid && info_out->pj && info_out->srid) {
93  char *insrid, *outsrid;
94 
95 #if PROJ_VERSION_MAJOR >= 6
96  /* PROJ6+: EPSG must uppercase EPSG */
97  if (strncmp(info_in->srid, "epsg", 4) == 0)
98  insrid = G_store_upper(info_in->srid);
99  else
100  insrid = G_store(info_in->srid);
101 
102  if (strncmp(info_out->srid, "epsg", 4) == 0)
103  outsrid = G_store_upper(info_out->srid);
104  else
105  outsrid = G_store(info_out->srid);
106 
107  /* PROJ6+: enforce axis order easting, northing
108  * +axis=enu (works with proj-4.8+) */
109 
110 #else
111  /* PROJ5: EPSG must lowercase epsg */
112  if (strncmp(info_in->srid, "EPSG", 4) == 0)
113  insrid = G_store_lower(info_in->srid);
114  else
115  insrid = G_store(info_in->srid);
116 
117  if (strncmp(info_out->srid, "EPSG", 4) == 0)
118  outsrid = G_store_lower(info_out->srid);
119  else
120  outsrid = G_store(info_out->srid);
121 
122 #endif
123  /* ask PROJ for the best pipeline */
124  info_trans->pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
125  insrid,
126  outsrid,
127  NULL);
128 
129  if (info_trans->pj == NULL) {
130  G_warning(_("proj_create_crs_to_crs() failed for '%s' and '%s'"),
131  insrid, outsrid);
132  }
133 #if PROJ_VERSION_MAJOR >= 6
134  else {
135  const char *str = proj_as_proj_string(NULL, info_trans->pj,
136  PJ_PROJ_5, NULL);
137 
138  if (str)
139  info_trans->def = G_store(str);
140  }
141 #endif
142  G_free(insrid);
143  G_free(outsrid);
144  }
145  if (info_trans->pj == NULL) {
146  /* PROJ6+: enforce axis order easting, northing
147  * +axis=enu (works with proj-4.8+) */
148  /* PROJ6+: what should we do with +towgs ?
149  * +towgs works only if WGS84 is used as pivot datum on both sides */
150  if (info_out->pj != NULL && info_out->def != NULL)
151  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s +step %s",
152  info_in->def, info_out->def);
153  else
154  /* assume info_out to be ll equivalent of info_in */
155  G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",
156  info_in->def);
157  }
158  }
159 
160  if (info_trans->pj == NULL)
161  info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
162  if (info_trans->pj == NULL) {
163  G_warning(_("proj_create() failed for '%s'"), info_trans->def);
164  return -1;
165  }
166 
167  info_trans->meters = 1.;
168  info_trans->zone = 0;
169  sprintf(info_trans->proj, "pipeline");
170 #else
171  if (info_out->pj == NULL) {
172  if (GPJ_get_equivalent_latlong(info_out, info_in) < 0) {
173  G_warning(_("Unable to create latlong equivalent for '%s'"),
174  info_in->def);
175  return -1;
176  }
177  }
178 #endif
179 
180  return 1;
181 }
182 
183 /* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
184 
185 /**
186  * \brief Re-project a point between two co-ordinate systems using a
187  * transformation object prepared with GPJ_prepare_pj()
188  *
189  * This function takes pointers to three pj_info structures as arguments,
190  * and projects a point between the input and output co-ordinate system.
191  * The pj_info structure info_trans must have been initialized with
192  * GPJ_init_transform().
193  * The direction determines if a point is projected from input CRS to
194  * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
195  * The easting, northing, and height of the point are contained in the
196  * pointers passed to the function; these will be overwritten by the
197  * coordinates of the transformed point.
198  *
199  * \param info_in pointer to pj_info struct for input co-ordinate system
200  * \param info_out pointer to pj_info struct for output co-ordinate system
201  * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
202  * \param dir direction of the transformation (PJ_FWD or PJ_INV)
203  * \param x Pointer to a double containing easting or longitude
204  * \param y Pointer to a double containing northing or latitude
205  * \param z Pointer to a double containing height, or NULL
206  *
207  * \return Return value from PROJ proj_trans() function
208  **/
209 
210 int GPJ_transform(const struct pj_info *info_in,
211  const struct pj_info *info_out,
212  const struct pj_info *info_trans, int dir,
213  double *x, double *y, double *z)
214 {
215  int ok = 0;
216 
217 #ifdef HAVE_PROJ_H
218  /* PROJ 5+ variant */
219  int in_is_ll, out_is_ll;
220  PJ_COORD c;
221 
222  if (info_in->pj == NULL)
223  G_fatal_error(_("No input projection"));
224 
225  if (info_trans->pj == NULL)
226  G_fatal_error(_("No transformation object"));
227 
228  if (dir == PJ_FWD) {
229  /* info_in -> info_out */
230  METERS_in = info_in->meters;
231  in_is_ll = !strncmp(info_in->proj, "ll", 2);
232  if (info_out->pj) {
233  METERS_out = info_out->meters;
234  out_is_ll = !strncmp(info_out->proj, "ll", 2);
235  }
236  else {
237  METERS_out = 1.0;
238  out_is_ll = 1;
239  }
240  }
241  else {
242  /* info_out -> info_in */
243  METERS_out = info_in->meters;
244  out_is_ll = !strncmp(info_in->proj, "ll", 2);
245  if (info_out->pj) {
246  METERS_in = info_out->meters;
247  in_is_ll = !strncmp(info_out->proj, "ll", 2);
248  }
249  else {
250  METERS_in = 1.0;
251  in_is_ll = 1;
252  }
253  }
254 
255  /* prepare */
256  if (in_is_ll) {
257  /* convert to radians */
258  /* PROJ 6: conversion to radians is not always needed:
259  * if proj_angular_input(info_trans->pj, dir) == 1
260  * -> convert from degrees to radians */
261  c.lpzt.lam = (*x) / RAD_TO_DEG;
262  c.lpzt.phi = (*y) / RAD_TO_DEG;
263  c.lpzt.z = 0;
264  if (z)
265  c.lpzt.z = *z;
266  c.lpzt.t = 0;
267  }
268  else {
269  /* convert to meters */
270  c.xyzt.x = *x * METERS_in;
271  c.xyzt.y = *y * METERS_in;
272  c.xyzt.z = 0;
273  if (z)
274  c.xyzt.z = *z;
275  c.xyzt.t = 0;
276  }
277 
278  /* transform */
279  c = proj_trans(info_trans->pj, dir, c);
280  ok = proj_errno(info_trans->pj);
281 
282  if (ok < 0) {
283  G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
284  return ok;
285  }
286 
287  /* output */
288  if (out_is_ll) {
289  /* convert to degrees */
290  /* PROJ 6: conversion to radians is not always needed:
291  * if proj_angular_output(info_trans->pj, dir) == 1
292  * -> convert from radians to degrees */
293  *x = c.lpzt.lam * RAD_TO_DEG;
294  *y = c.lpzt.phi * RAD_TO_DEG;
295  if (z)
296  *z = c.lpzt.z;
297  }
298  else {
299  /* convert to map units */
300  *x = c.xyzt.x / METERS_out;
301  *y = c.xyzt.y / METERS_out;
302  if (z)
303  *z = c.xyzt.z;
304  }
305 #else
306  /* PROJ 4 variant */
307  double u, v;
308  double h = 0.0;
309  const struct pj_info *p_in, *p_out;
310 
311  if (info_out == NULL)
312  G_fatal_error(_("No output projection"));
313 
314  if (dir == PJ_FWD) {
315  p_in = info_in;
316  p_out = info_out;
317  }
318  else {
319  p_in = info_out;
320  p_out = info_in;
321  }
322 
323  METERS_in = p_in->meters;
324  METERS_out = p_out->meters;
325 
326  if (z)
327  h = *z;
328 
329  if (strncmp(p_in->proj, "ll", 2) == 0) {
330  u = (*x) / RAD_TO_DEG;
331  v = (*y) / RAD_TO_DEG;
332  }
333  else {
334  u = *x * METERS_in;
335  v = *y * METERS_in;
336  }
337 
338  ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
339 
340  if (ok < 0) {
341  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
342  return ok;
343  }
344 
345  if (strncmp(p_out->proj, "ll", 2) == 0) {
346  *x = u * RAD_TO_DEG;
347  *y = v * RAD_TO_DEG;
348  }
349  else {
350  *x = u / METERS_out;
351  *y = v / METERS_out;
352  }
353  if (z)
354  *z = h;
355 #endif
356 
357  return ok;
358 }
359 
360 /**
361  * \brief Re-project an array of points between two co-ordinate systems
362  * using a transformation object prepared with GPJ_prepare_pj()
363  *
364  * This function takes pointers to three pj_info structures as arguments,
365  * and projects an array of pointd between the input and output
366  * co-ordinate system. The pj_info structure info_trans must have been
367  * initialized with GPJ_init_transform().
368  * The direction determines if a point is projected from input CRS to
369  * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
370  * The easting, northing, and height of the point are contained in the
371  * pointers passed to the function; these will be overwritten by the
372  * coordinates of the transformed point.
373  *
374  * \param info_in pointer to pj_info struct for input co-ordinate system
375  * \param info_out pointer to pj_info struct for output co-ordinate system
376  * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
377  * \param dir direction of the transformation (PJ_FWD or PJ_INV)
378  * \param x pointer to an array of type double containing easting or longitude
379  * \param y pointer to an array of type double containing northing or latitude
380  * \param z pointer to an array of type double containing height, or NULL
381  * \param n number of points in the arrays to be transformed
382  *
383  * \return Return value from PROJ proj_trans() function
384  **/
385 
386 int GPJ_transform_array(const struct pj_info *info_in,
387  const struct pj_info *info_out,
388  const struct pj_info *info_trans, int dir,
389  double *x, double *y, double *z, int n)
390 {
391  int ok;
392  int i;
393  int has_z = 1;
394 
395 #ifdef HAVE_PROJ_H
396  /* PROJ 5+ variant */
397  int in_is_ll, out_is_ll;
398  PJ_COORD c;
399 
400  if (info_trans->pj == NULL)
401  G_fatal_error(_("No transformation object"));
402 
403  if (dir == PJ_FWD) {
404  /* info_in -> info_out */
405  METERS_in = info_in->meters;
406  in_is_ll = !strncmp(info_in->proj, "ll", 2);
407  if (info_out->pj) {
408  METERS_out = info_out->meters;
409  out_is_ll = !strncmp(info_out->proj, "ll", 2);
410  }
411  else {
412  METERS_out = 1.0;
413  out_is_ll = 1;
414  }
415  }
416  else {
417  /* info_out -> info_in */
418  METERS_out = info_in->meters;
419  out_is_ll = !strncmp(info_in->proj, "ll", 2);
420  if (info_out->pj) {
421  METERS_in = info_out->meters;
422  in_is_ll = !strncmp(info_out->proj, "ll", 2);
423  }
424  else {
425  METERS_in = 1.0;
426  in_is_ll = 1;
427  }
428  }
429 
430  if (z == NULL) {
431  z = G_malloc(sizeof(double) * n);
432  /* they say memset is only guaranteed for chars ;-( */
433  for (i = 0; i < n; i++)
434  z[i] = 0.0;
435  has_z = 0;
436  }
437  ok = 0;
438  if (in_is_ll) {
439  c.lpzt.t = 0;
440  if (out_is_ll) {
441  /* what is more costly ?
442  * calling proj_trans for each point
443  * or having three loops over all points ?
444  * proj_trans_array() itself calls proj_trans() in a loop
445  * -> one loop over all points is better than
446  * three loops over all points
447  */
448  for (i = 0; i < n; i++) {
449  /* convert to radians */
450  c.lpzt.lam = x[i] / RAD_TO_DEG;
451  c.lpzt.phi = y[i] / RAD_TO_DEG;
452  c.lpzt.z = z[i];
453  c = proj_trans(info_trans->pj, dir, c);
454  if ((ok = proj_errno(info_trans->pj)) < 0)
455  break;
456  /* convert to degrees */
457  x[i] = c.lp.lam * RAD_TO_DEG;
458  y[i] = c.lp.phi * RAD_TO_DEG;
459  }
460  }
461  else {
462  for (i = 0; i < n; i++) {
463  /* convert to radians */
464  c.lpzt.lam = x[i] / RAD_TO_DEG;
465  c.lpzt.phi = y[i] / RAD_TO_DEG;
466  c.lpzt.z = z[i];
467  c = proj_trans(info_trans->pj, dir, c);
468  if ((ok = proj_errno(info_trans->pj)) < 0)
469  break;
470  /* convert to map units */
471  x[i] = c.xy.x / METERS_out;
472  y[i] = c.xy.y / METERS_out;
473  }
474  }
475  }
476  else {
477  c.xyzt.t = 0;
478  if (out_is_ll) {
479  for (i = 0; i < n; i++) {
480  /* convert to meters */
481  c.xyzt.x = x[i] * METERS_in;
482  c.xyzt.y = y[i] * METERS_in;
483  c.xyzt.z = z[i];
484  c = proj_trans(info_trans->pj, dir, c);
485  if ((ok = proj_errno(info_trans->pj)) < 0)
486  break;
487  /* convert to degrees */
488  x[i] = c.lp.lam * RAD_TO_DEG;
489  y[i] = c.lp.phi * RAD_TO_DEG;
490  }
491  }
492  else {
493  for (i = 0; i < n; i++) {
494  /* convert to meters */
495  c.xyzt.x = x[i] * METERS_in;
496  c.xyzt.y = y[i] * METERS_in;
497  c.xyzt.z = z[i];
498  c = proj_trans(info_trans->pj, dir, c);
499  if ((ok = proj_errno(info_trans->pj)) < 0)
500  break;
501  /* convert to map units */
502  x[i] = c.xy.x / METERS_out;
503  y[i] = c.xy.y / METERS_out;
504  }
505  }
506  }
507  if (!has_z)
508  G_free(z);
509 
510  if (ok < 0) {
511  G_warning(_("proj_trans() failed: %s"), proj_errno_string(ok));
512  }
513 #else
514  /* PROJ 4 variant */
515  const struct pj_info *p_in, *p_out;
516 
517  if (dir == PJ_FWD) {
518  p_in = info_in;
519  p_out = info_out;
520  }
521  else {
522  p_in = info_out;
523  p_out = info_in;
524  }
525 
526  METERS_in = p_in->meters;
527  METERS_out = p_out->meters;
528 
529  if (z == NULL) {
530  z = G_malloc(sizeof(double) * n);
531  /* they say memset is only guaranteed for chars ;-( */
532  for (i = 0; i < n; ++i)
533  z[i] = 0.0;
534  has_z = 0;
535  }
536  if (strncmp(p_in->proj, "ll", 2) == 0) {
537  if (strncmp(p_out->proj, "ll", 2) == 0) {
538  DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
539  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
540  MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
541  }
542  else {
543  DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
544  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
545  DIVIDE_LOOP(x, y, n, METERS_out);
546  }
547  }
548  else {
549  if (strncmp(p_out->proj, "ll", 2) == 0) {
550  MULTIPLY_LOOP(x, y, n, METERS_in);
551  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
552  MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
553  }
554  else {
555  MULTIPLY_LOOP(x, y, n, METERS_in);
556  ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
557  DIVIDE_LOOP(x, y, n, METERS_out);
558  }
559  }
560  if (!has_z)
561  G_free(z);
562 
563  if (ok < 0)
564  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
565 #endif
566 
567  return ok;
568 }
569 
570 /*
571  * old API, to be deleted
572  */
573 
574 /**
575  * \brief Re-project a point between two co-ordinate systems
576  *
577  * This function takes pointers to two pj_info structures as arguments,
578  * and projects a point between the co-ordinate systems represented by them.
579  * The easting and northing of the point are contained in two pointers passed
580  * to the function; these will be overwritten by the co-ordinates of the
581  * re-projected point.
582  *
583  * \param x Pointer to a double containing easting or longitude
584  * \param y Pointer to a double containing northing or latitude
585  * \param info_in pointer to pj_info struct for input co-ordinate system
586  * \param info_out pointer to pj_info struct for output co-ordinate system
587  *
588  * \return Return value from PROJ proj_trans() function
589  **/
590 
591 int pj_do_proj(double *x, double *y,
592  const struct pj_info *info_in, const struct pj_info *info_out)
593 {
594  int ok;
595 #ifdef HAVE_PROJ_H
596  struct pj_info info_trans;
597  PJ_COORD c;
598 
599  if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
600  return -1;
601  }
602 
603  METERS_in = info_in->meters;
604  METERS_out = info_out->meters;
605 
606  if (strncmp(info_in->proj, "ll", 2) == 0) {
607  /* convert to radians */
608  c.lpzt.lam = (*x) / RAD_TO_DEG;
609  c.lpzt.phi = (*y) / RAD_TO_DEG;
610  c.lpzt.z = 0;
611  c.lpzt.t = 0;
612  c = proj_trans(info_trans.pj, PJ_FWD, c);
613  ok = proj_errno(info_trans.pj);
614 
615  if (strncmp(info_out->proj, "ll", 2) == 0) {
616  /* convert to degrees */
617  *x = c.lp.lam * RAD_TO_DEG;
618  *y = c.lp.phi * RAD_TO_DEG;
619  }
620  else {
621  /* convert to map units */
622  *x = c.xy.x / METERS_out;
623  *y = c.xy.y / METERS_out;
624  }
625  }
626  else {
627  /* convert to meters */
628  c.xyzt.x = *x * METERS_in;
629  c.xyzt.y = *y * METERS_in;
630  c.xyzt.z = 0;
631  c.xyzt.t = 0;
632  c = proj_trans(info_trans.pj, PJ_FWD, c);
633  ok = proj_errno(info_trans.pj);
634 
635  if (strncmp(info_out->proj, "ll", 2) == 0) {
636  /* convert to degrees */
637  *x = c.lp.lam * RAD_TO_DEG;
638  *y = c.lp.phi * RAD_TO_DEG;
639  }
640  else {
641  /* convert to map units */
642  *x = c.xy.x / METERS_out;
643  *y = c.xy.y / METERS_out;
644  }
645  }
646  proj_destroy(info_trans.pj);
647 
648  if (ok < 0) {
649  G_warning(_("proj_trans() failed: %d"), ok);
650  }
651 #else
652  double u, v;
653  double h = 0.0;
654 
655  METERS_in = info_in->meters;
656  METERS_out = info_out->meters;
657 
658  if (strncmp(info_in->proj, "ll", 2) == 0) {
659  if (strncmp(info_out->proj, "ll", 2) == 0) {
660  u = (*x) / RAD_TO_DEG;
661  v = (*y) / RAD_TO_DEG;
662  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
663  *x = u * RAD_TO_DEG;
664  *y = v * RAD_TO_DEG;
665  }
666  else {
667  u = (*x) / RAD_TO_DEG;
668  v = (*y) / RAD_TO_DEG;
669  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
670  *x = u / METERS_out;
671  *y = v / METERS_out;
672  }
673  }
674  else {
675  if (strncmp(info_out->proj, "ll", 2) == 0) {
676  u = *x * METERS_in;
677  v = *y * METERS_in;
678  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
679  *x = u * RAD_TO_DEG;
680  *y = v * RAD_TO_DEG;
681  }
682  else {
683  u = *x * METERS_in;
684  v = *y * METERS_in;
685  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
686  *x = u / METERS_out;
687  *y = v / METERS_out;
688  }
689  }
690  if (ok < 0) {
691  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
692  }
693 #endif
694  return ok;
695 }
696 
697 /**
698  * \brief Re-project an array of points between two co-ordinate systems with
699  * optional ellipsoidal height conversion
700  *
701  * This function takes pointers to two pj_info structures as arguments,
702  * and projects an array of points between the co-ordinate systems
703  * represented by them. Pointers to the three arrays of easting, northing,
704  * and ellipsoidal height of the point (this one may be NULL) are passed
705  * to the function; these will be overwritten by the co-ordinates of the
706  * re-projected points.
707  *
708  * \param count Number of points in the arrays to be transformed
709  * \param x Pointer to an array of type double containing easting or longitude
710  * \param y Pointer to an array of type double containing northing or latitude
711  * \param h Pointer to an array of type double containing ellipsoidal height.
712  * May be null in which case a two-dimensional re-projection will be
713  * done
714  * \param info_in pointer to pj_info struct for input co-ordinate system
715  * \param info_out pointer to pj_info struct for output co-ordinate system
716  *
717  * \return Return value from PROJ proj_trans() function
718  **/
719 
720 int pj_do_transform(int count, double *x, double *y, double *h,
721  const struct pj_info *info_in, const struct pj_info *info_out)
722 {
723  int ok;
724  int i;
725  int has_h = 1;
726 #ifdef HAVE_PROJ_H
727  struct pj_info info_trans;
728  PJ_COORD c;
729 
730  if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
731  return -1;
732  }
733 
734  METERS_in = info_in->meters;
735  METERS_out = info_out->meters;
736 
737  if (h == NULL) {
738  h = G_malloc(sizeof *h * count);
739  /* they say memset is only guaranteed for chars ;-( */
740  for (i = 0; i < count; ++i)
741  h[i] = 0.0;
742  has_h = 0;
743  }
744  ok = 0;
745  if (strncmp(info_in->proj, "ll", 2) == 0) {
746  c.lpzt.t = 0;
747  if (strncmp(info_out->proj, "ll", 2) == 0) {
748  for (i = 0; i < count; i++) {
749  /* convert to radians */
750  c.lpzt.lam = x[i] / RAD_TO_DEG;
751  c.lpzt.phi = y[i] / RAD_TO_DEG;
752  c.lpzt.z = h[i];
753  c = proj_trans(info_trans.pj, PJ_FWD, c);
754  if ((ok = proj_errno(info_trans.pj)) < 0)
755  break;
756  /* convert to degrees */
757  x[i] = c.lp.lam * RAD_TO_DEG;
758  y[i] = c.lp.phi * RAD_TO_DEG;
759  }
760  }
761  else {
762  for (i = 0; i < count; i++) {
763  /* convert to radians */
764  c.lpzt.lam = x[i] / RAD_TO_DEG;
765  c.lpzt.phi = y[i] / RAD_TO_DEG;
766  c.lpzt.z = h[i];
767  c = proj_trans(info_trans.pj, PJ_FWD, c);
768  if ((ok = proj_errno(info_trans.pj)) < 0)
769  break;
770  /* convert to map units */
771  x[i] = c.xy.x / METERS_out;
772  y[i] = c.xy.y / METERS_out;
773  }
774  }
775  }
776  else {
777  c.xyzt.t = 0;
778  if (strncmp(info_out->proj, "ll", 2) == 0) {
779  for (i = 0; i < count; i++) {
780  /* convert to meters */
781  c.xyzt.x = x[i] * METERS_in;
782  c.xyzt.y = y[i] * METERS_in;
783  c.xyzt.z = h[i];
784  c = proj_trans(info_trans.pj, PJ_FWD, c);
785  if ((ok = proj_errno(info_trans.pj)) < 0)
786  break;
787  /* convert to degrees */
788  x[i] = c.lp.lam * RAD_TO_DEG;
789  y[i] = c.lp.phi * RAD_TO_DEG;
790  }
791  }
792  else {
793  for (i = 0; i < count; i++) {
794  /* convert to meters */
795  c.xyzt.x = x[i] * METERS_in;
796  c.xyzt.y = y[i] * METERS_in;
797  c.xyzt.z = h[i];
798  c = proj_trans(info_trans.pj, PJ_FWD, c);
799  if ((ok = proj_errno(info_trans.pj)) < 0)
800  break;
801  /* convert to map units */
802  x[i] = c.xy.x / METERS_out;
803  y[i] = c.xy.y / METERS_out;
804  }
805  }
806  }
807  if (!has_h)
808  G_free(h);
809  proj_destroy(info_trans.pj);
810 
811  if (ok < 0) {
812  G_warning(_("proj_trans() failed: %d"), ok);
813  }
814 #else
815  METERS_in = info_in->meters;
816  METERS_out = info_out->meters;
817 
818  if (h == NULL) {
819  h = G_malloc(sizeof *h * count);
820  /* they say memset is only guaranteed for chars ;-( */
821  for (i = 0; i < count; ++i)
822  h[i] = 0.0;
823  has_h = 0;
824  }
825  if (strncmp(info_in->proj, "ll", 2) == 0) {
826  if (strncmp(info_out->proj, "ll", 2) == 0) {
827  DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
828  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
829  MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
830  }
831  else {
832  DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
833  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
834  DIVIDE_LOOP(x, y, count, METERS_out);
835  }
836  }
837  else {
838  if (strncmp(info_out->proj, "ll", 2) == 0) {
839  MULTIPLY_LOOP(x, y, count, METERS_in);
840  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
841  MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
842  }
843  else {
844  MULTIPLY_LOOP(x, y, count, METERS_in);
845  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
846  DIVIDE_LOOP(x, y, count, METERS_out);
847  }
848  }
849  if (!has_h)
850  G_free(h);
851 
852  if (ok < 0) {
853  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
854  }
855 #endif
856  return ok;
857 }
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:720
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:591
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:458
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:82
int count
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
#define DIVIDE_LOOP(x, y, c, m)
Definition: do_proj.c:34
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_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:386
char * G_store_upper(const char *s)
Copy string to allocated memory and convert copied string to upper case.
Definition: strings.c:116
#define MULTIPLY_LOOP(x, y, c, m)
Definition: do_proj.c:26
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:210
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204