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