GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72090
 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
9 
10  (C) 2003-2008 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 Re-project a point between two co-ordinate systems
48  *
49  * This function takes pointers to two pj_info structures as arguments,
50  * and projects a point between the co-ordinate systems represented by them.
51  * The easting and northing of the point are contained in two pointers passed
52  * to the function; these will be overwritten by the co-ordinates of the
53  * re-projected point.
54  *
55  * \param x Pointer to a double containing easting or longitude
56  * \param y Pointer to a double containing northing or latitude
57  * \param info_in pointer to pj_info struct for input co-ordinate system
58  * \param info_out pointer to pj_info struct for output co-ordinate system
59  *
60  * \return Return value from PROJ pj_transform() function
61  **/
62 
63 int pj_do_proj(double *x, double *y,
64  const struct pj_info *info_in, const struct pj_info *info_out)
65 {
66  int ok;
67  double u, v;
68  double h = 0.0;
69 
70  METERS_in = info_in->meters;
71  METERS_out = info_out->meters;
72 
73  if (strncmp(info_in->proj, "ll", 2) == 0) {
74  if (strncmp(info_out->proj, "ll", 2) == 0) {
75  u = (*x) / RAD_TO_DEG;
76  v = (*y) / RAD_TO_DEG;
77  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
78  *x = u * RAD_TO_DEG;
79  *y = v * RAD_TO_DEG;
80  }
81  else {
82  u = (*x) / RAD_TO_DEG;
83  v = (*y) / RAD_TO_DEG;
84  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
85  *x = u / METERS_out;
86  *y = v / METERS_out;
87  }
88  }
89  else {
90  if (strncmp(info_out->proj, "ll", 2) == 0) {
91  u = *x * METERS_in;
92  v = *y * METERS_in;
93  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
94  *x = u * RAD_TO_DEG;
95  *y = v * RAD_TO_DEG;
96  }
97  else {
98  u = *x * METERS_in;
99  v = *y * METERS_in;
100  ok = pj_transform(info_in->pj, info_out->pj, 1, 0, &u, &v, &h);
101  *x = u / METERS_out;
102  *y = v / METERS_out;
103  }
104  }
105  if (ok < 0) {
106  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
107  }
108  return ok;
109 }
110 
111 /**
112  * \brief Re-project an array of points between two co-ordinate systems with
113  * optional ellipsoidal height conversion
114  *
115  * This function takes pointers to two pj_info structures as arguments,
116  * and projects an array of points between the co-ordinate systems
117  * represented by them. Pointers to the three arrays of easting, northing,
118  * and ellipsoidal height of the point (this one may be NULL) are passed
119  * to the function; these will be overwritten by the co-ordinates of the
120  * re-projected points.
121  *
122  * \param count Number of points in the arrays to be transformed
123  * \param x Pointer to an array of type double containing easting or longitude
124  * \param y Pointer to an array of type double containing northing or latitude
125  * \param h Pointer to an array of type double containing ellipsoidal height.
126  * May be null in which case a two-dimensional re-projection will be
127  * done
128  * \param info_in pointer to pj_info struct for input co-ordinate system
129  * \param info_out pointer to pj_info struct for output co-ordinate system
130  *
131  * \return Return value from PROJ pj_transform() function
132  **/
133 
134 int pj_do_transform(int count, double *x, double *y, double *h,
135  const struct pj_info *info_in, const struct pj_info *info_out)
136 {
137  int ok;
138  int has_h = 1;
139 
140  METERS_in = info_in->meters;
141  METERS_out = info_out->meters;
142 
143  if (h == NULL) {
144  int i;
145 
146  h = G_malloc(sizeof *h * count);
147  /* they say memset is only guaranteed for chars ;-( */
148  for (i = 0; i < count; ++i)
149  h[i] = 0.0;
150  has_h = 0;
151  }
152  if (strncmp(info_in->proj, "ll", 2) == 0) {
153  if (strncmp(info_out->proj, "ll", 2) == 0) {
154  DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
155  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
156  MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
157  }
158  else {
159  DIVIDE_LOOP(x, y, count, RAD_TO_DEG);
160  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
161  DIVIDE_LOOP(x, y, count, METERS_out);
162  }
163  }
164  else {
165  if (strncmp(info_out->proj, "ll", 2) == 0) {
166  MULTIPLY_LOOP(x, y, count, METERS_in);
167  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
168  MULTIPLY_LOOP(x, y, count, RAD_TO_DEG);
169  }
170  else {
171  MULTIPLY_LOOP(x, y, count, METERS_in);
172  ok = pj_transform(info_in->pj, info_out->pj, count, 1, x, y, h);
173  DIVIDE_LOOP(x, y, count, METERS_out);
174  }
175  }
176  if (!has_h)
177  G_free(h);
178 
179  if (ok < 0) {
180  G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
181  }
182  return ok;
183 }
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:134
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:63
char proj[100]
Definition: gprojects.h:38
int count
#define DIVIDE_LOOP(x, y, c, m)
Definition: do_proj.c:35
#define NULL
Definition: ccmath.h:32
#define x
projPJ pj
Definition: gprojects.h:35
double meters
Definition: gprojects.h:36
#define MULTIPLY_LOOP(x, y, c, m)
Definition: do_proj.c:26
#define _(str)
Definition: glocale.h:13
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203