GRASS GIS 7 Programmer's Manual  7.2.0(2016)-exported
get_proj.c
Go to the documentation of this file.
1 
2 /**
3  \file get_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 <stdlib.h>
19 #include <ctype.h>
20 #include <math.h>
21 #include <string.h>
22 #include <grass/gis.h>
23 #include <grass/gprojects.h>
24 #include <grass/glocale.h>
25 
26 /* Finder function for datum conversion lookup tables */
27 #define FINDERFUNC set_proj_lib
28 #define PERMANENT "PERMANENT"
29 #define MAX_PARGS 100
30 
31 static void alloc_options(char *);
32 
33 static char *opt_in[MAX_PARGS];
34 static int nopt1;
35 
36 /**
37  * \brief Create a pj_info struct Co-ordinate System definition from a set of
38  * PROJ_INFO / PROJ_UNITS-style key-value pairs
39  *
40  * This function takes a GRASS-style co-ordinate system definition as stored
41  * in the PROJ_INFO and PROJ_UNITS files and processes it to create a pj_info
42  * representation for use in re-projecting with pj_do_proj(). In addition to
43  * the parameters passed to it it may also make reference to the system
44  * ellipse.table and datum.table files if necessary.
45  *
46  * \param info Pointer to a pj_info struct (which must already exist) into
47  * which the co-ordinate system definition will be placed
48  * \param in_proj_keys PROJ_INFO-style key-value pairs
49  * \param in_units_keys PROJ_UNITS-style key-value pairs
50  *
51  * \return -1 on error (unable to initialise PROJ.4)
52  * 2 if "default" 3-parameter datum shift values from datum.table
53  * were used
54  * 3 if an unrecognised datum name was passed on to PROJ.4 (and
55  * initialization was successful
56  * 1 otherwise
57  **/
58 
59 int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
60  const struct Key_Value *in_units_keys)
61 {
62  const char *str;
63  int i;
64  double a, es, rf;
65  int returnval = 1;
66  char buffa[300], factbuff[50];
67  char proj_in[50], *datum, *params;
68  projPJ *pj;
69 
70  proj_in[0] = '\0';
71  info->zone = 0;
72  info->meters = 1.0;
73  info->proj[0] = '\0';
74 
75  str = G_find_key_value("meters", in_units_keys);
76  if (str != NULL) {
77  strcpy(factbuff, str);
78  if (strlen(factbuff) > 0)
79  sscanf(factbuff, "%lf", &(info->meters));
80  }
81  str = G_find_key_value("name", in_proj_keys);
82  if (str != NULL) {
83  sprintf(proj_in, "%s", str);
84  }
85  str = G_find_key_value("proj", in_proj_keys);
86  if (str != NULL) {
87  sprintf(info->proj, "%s", str);
88  }
89  if (strlen(info->proj) <= 0)
90  sprintf(info->proj, "ll");
91 
92  nopt1 = 0;
93  for (i = 0; i < in_proj_keys->nitems; i++) {
94  /* the name parameter is just for grasses use */
95  if (strcmp(in_proj_keys->key[i], "name") == 0) {
96  continue;
97 
98  /* zone handled separately at end of loop */
99  }
100  else if (strcmp(in_proj_keys->key[i], "zone") == 0) {
101  continue;
102 
103  /* Datum and ellipsoid-related parameters will be handled
104  * separately after end of this loop PK */
105 
106  }
107  else if (strcmp(in_proj_keys->key[i], "datum") == 0
108  || strcmp(in_proj_keys->key[i], "dx") == 0
109  || strcmp(in_proj_keys->key[i], "dy") == 0
110  || strcmp(in_proj_keys->key[i], "dz") == 0
111  || strcmp(in_proj_keys->key[i], "datumparams") == 0
112  || strcmp(in_proj_keys->key[i], "nadgrids") == 0
113  || strcmp(in_proj_keys->key[i], "towgs84") == 0
114  || strcmp(in_proj_keys->key[i], "ellps") == 0
115  || strcmp(in_proj_keys->key[i], "a") == 0
116  || strcmp(in_proj_keys->key[i], "b") == 0
117  || strcmp(in_proj_keys->key[i], "es") == 0
118  || strcmp(in_proj_keys->key[i], "f") == 0
119  || strcmp(in_proj_keys->key[i], "rf") == 0) {
120  continue;
121 
122  /* PROJ.4 uses longlat instead of ll as 'projection name' */
123 
124  }
125  else if (strcmp(in_proj_keys->key[i], "proj") == 0) {
126  if (strcmp(in_proj_keys->value[i], "ll") == 0)
127  sprintf(buffa, "proj=longlat");
128  else
129  sprintf(buffa, "proj=%s", in_proj_keys->value[i]);
130 
131  /* 'One-sided' PROJ.4 flags will have the value in
132  * the key-value pair set to 'defined' and only the
133  * key needs to be passed on. */
134  }
135  else if (strcmp(in_proj_keys->value[i], "defined") == 0)
136  sprintf(buffa, "%s", in_proj_keys->key[i]);
137 
138  else
139  sprintf(buffa, "%s=%s",
140  in_proj_keys->key[i], in_proj_keys->value[i]);
141 
142  alloc_options(buffa);
143  }
144 
145  str = G_find_key_value("zone", in_proj_keys);
146  if (str != NULL) {
147  if (sscanf(str, "%d", &(info->zone)) != 1) {
148  G_fatal_error(_("Invalid zone %s specified"), str);
149  }
150  if (info->zone < 0) {
151 
152  /* if zone is negative, write abs(zone) and define south */
153  info->zone = -info->zone;
154 
155  if (G_find_key_value("south", in_proj_keys) == NULL) {
156  sprintf(buffa, "south");
157  alloc_options(buffa);
158  }
159  }
160  sprintf(buffa, "zone=%d", info->zone);
161  alloc_options(buffa);
162  }
163 
164  if ((GPJ__get_ellipsoid_params(in_proj_keys, &a, &es, &rf) == 0)
165  && (str = G_find_key_value("ellps", in_proj_keys)) != NULL) {
166  /* Default values were returned but an ellipsoid name not recognised
167  * by GRASS is present---perhaps it will be recognised by
168  * PROJ.4 even though it wasn't by GRASS */
169  sprintf(buffa, "ellps=%s", str);
170  alloc_options(buffa);
171  }
172  else {
173  sprintf(buffa, "a=%.16g", a);
174  alloc_options(buffa);
175  /* Cannot use es directly because the OSRImportFromProj4()
176  * function in OGR only accepts b or rf as the 2nd parameter */
177  if (es == 0)
178  sprintf(buffa, "b=%.16g", a);
179  else
180  sprintf(buffa, "rf=%.16g", rf);
181  alloc_options(buffa);
182 
183  }
184  /* Workaround to stop PROJ reading values from defaults file when
185  * rf (and sometimes ellps) is not specified */
186  if (G_find_key_value("no_defs", in_proj_keys) == NULL) {
187  sprintf(buffa, "no_defs");
188  alloc_options(buffa);
189  }
190 
191  /* If datum parameters are present in the PROJ_INFO keys, pass them on */
192  if (GPJ__get_datum_params(in_proj_keys, &datum, &params) == 2) {
193  sprintf(buffa, "%s", params);
194  alloc_options(buffa);
195  G_free(params);
196 
197  /* else if a datum name is present take it and look up the parameters
198  * from the datum.table file */
199  }
200  else if (datum != NULL) {
201 
202  if (GPJ_get_default_datum_params_by_name(datum, &params) > 0) {
203  sprintf(buffa, "%s", params);
204  alloc_options(buffa);
205  returnval = 2;
206  G_free(params);
207 
208  /* else just pass the datum name on and hope it is recognised by
209  * PROJ.4 even though it isn't recognised by GRASS */
210  }
211  else {
212  sprintf(buffa, "datum=%s", datum);
213  alloc_options(buffa);
214  returnval = 3;
215  }
216  /* else there'll be no datum transformation taking place here... */
217  }
218  else {
219  returnval = 4;
220  }
221  G_free(datum);
222 
223  /* Set finder function for locating datum conversion tables PK */
224  pj_set_finder(FINDERFUNC);
225 
226  if (!(pj = pj_init(nopt1, opt_in))) {
227  strcpy(buffa,
228  _("Unable to initialise PROJ.4 with the following parameter list:"));
229  for (i = 0; i < nopt1; i++) {
230  char err[50];
231 
232  sprintf(err, " +%s", opt_in[i]);
233  strcat(buffa, err);
234  }
235  G_warning("%s", buffa);
236  G_warning(_("The error message: %s"), pj_strerrno(pj_errno));
237  return -1;
238  }
239  info->pj = pj;
240 
241  return returnval;
242 }
243 
244 static void alloc_options(char *buffa)
245 {
246  int nsize;
247 
248  nsize = strlen(buffa);
249  opt_in[nopt1++] = (char *)G_malloc(nsize + 1);
250  sprintf(opt_in[nopt1 - 1], "%s", buffa);
251  return;
252 }
253 
254 int pj_get_string(struct pj_info *info, char *str)
255 {
256  char *opt_in[MAX_PARGS];
257  char *s;
258  int nopt = 0;
259  int nsize;
260  char zonebuff[50], buffa[300];
261  projPJ *pj;
262 
263  info->zone = 0;
264  info->proj[0] = '\0';
265  info->meters = 1.0;
266 
267  if ((str == NULL) || (str[0] == '\0')) {
268  /* Null Pointer or empty string is supplied for parameters,
269  * implying latlong projection; just need to set proj
270  * parameter and call pj_init PK */
271  sprintf(info->proj, "ll");
272  sprintf(buffa, "proj=latlong ellps=WGS84");
273  nsize = strlen(buffa);
274  opt_in[nopt] = (char *)G_malloc(nsize + 1);
275  sprintf(opt_in[nopt++], "%s", buffa);
276  }
277  else {
278  /* Parameters have been provided; parse through them but don't
279  * bother with most of the checks in pj_get_kv; assume the
280  * programmer knows what he / she is doing when using this
281  * function rather than reading a PROJ_INFO file PK */
282  s = str;
283  while (s = strtok(s, " \t\n"), s) {
284  if (strncmp(s, "+unfact=", 8) == 0) {
285  s = s + 8;
286  info->meters = atof(s);
287  }
288  else {
289  if (strncmp(s, "+", 1) == 0)
290  ++s;
291  if (nsize = strlen(s), nsize) {
292  if (nopt >= MAX_PARGS) {
293  fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
294  G_fatal_error(_("Option input overflowed option table"));
295  }
296 
297  if (strncmp("zone=", s, 5) == 0) {
298  sprintf(zonebuff, "%s", s + 5);
299  sscanf(zonebuff, "%d", &(info->zone));
300  }
301 
302  if (strncmp("proj=", s, 5) == 0) {
303  sprintf(info->proj, "%s", s + 5);
304  if (strcmp(info->proj, "ll") == 0)
305  sprintf(buffa, "proj=latlong");
306  else
307  sprintf(buffa, "%s", s);
308  }
309  else {
310  sprintf(buffa, "%s", s);
311  }
312  nsize = strlen(buffa);
313  opt_in[nopt] = (char *)G_malloc(nsize + 1);
314  sprintf(opt_in[nopt++], "%s", buffa);
315  }
316  }
317  s = 0;
318  }
319  }
320 
321  /* Set finder function for locating datum conversion tables PK */
322  pj_set_finder(FINDERFUNC);
323 
324  if (!(pj = pj_init(nopt, opt_in))) {
325  G_warning(_("Unable to initialize pj cause: %s"),
326  pj_strerrno(pj_errno));
327  return -1;
328  }
329  info->pj = pj;
330 
331  return 1;
332 }
333 
334 /**
335  * \brief Define a latitude / longitude co-ordinate system with the same
336  * ellipsoid and datum parameters as an existing projected system
337  *
338  * This function is useful when projected co-ordinates need to be simply
339  * converted to and from latitude / longitude.
340  *
341  * \param pjnew Pointer to pj_info struct for geographic co-ordinate system
342  * that will be created
343  * \param pjold Pointer to pj_info struct for existing projected co-ordinate
344  * system
345  *
346  * \return 1 on success; -1 if there was an error (i.e. if the PROJ.4
347  * pj_latlong_from_proj() function returned NULL)
348  **/
349 
350 int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
351 {
352  pjnew->meters = 1.;
353  pjnew->zone = 0;
354  sprintf(pjnew->proj, "ll");
355  if ((pjnew->pj = pj_latlong_from_proj(pjold->pj)) == NULL)
356  return -1;
357  else
358  return 1;
359 }
360 
361 /* set_proj_lib()
362  * 'finder function' for use with PROJ.4 pj_set_finder() function */
363 
364 const char *set_proj_lib(const char *name)
365 {
366  const char *gisbase = G_gisbase();
367  static char *buf = NULL;
368  static size_t buf_len;
369  size_t len = strlen(gisbase) + sizeof(GRIDDIR) + strlen(name) + 1;
370 
371  if (buf_len < len) {
372  if (buf != NULL)
373  G_free(buf);
374  buf_len = len + 20;
375  buf = G_malloc(buf_len);
376  }
377 
378  sprintf(buf, "%s%s/%s", gisbase, GRIDDIR, name);
379 
380  return buf;
381 }
382 
383 /**
384  * \brief Print projection parameters as used by PROJ.4 for input and
385  * output co-ordinate systems
386  *
387  * \param iproj 'Input' co-ordinate system
388  * \param oproj 'Output' co-ordinate system
389  *
390  * \return 1 on success, -1 on error (i.e. if PROJ.4 pj_get_def() function
391  * returned NULL for either co-ordinate system)
392  **/
393 
394 int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *oproj)
395 {
396  char *str;
397 
398  if (iproj) {
399  str = pj_get_def(iproj->pj, 1);
400  if (str != NULL) {
401  fprintf(stderr, "%s: %s\n", _("Input Projection Parameters"),
402  str);
403  pj_dalloc(str);
404  fprintf(stderr, "%s: %.16g\n", _("Input Unit Factor"),
405  iproj->meters);
406  }
407  else
408  return -1;
409  }
410 
411  if (oproj) {
412  str = pj_get_def(oproj->pj, 1);
413  if (str != NULL) {
414  fprintf(stderr, "%s: %s\n", _("Output Projection Parameters"),
415  str);
416  pj_dalloc(str);
417  fprintf(stderr, "%s: %.16g\n", _("Output Unit Factor"),
418  oproj->meters);
419  }
420  else
421  return -1;
422  }
423 
424  return 1;
425 }
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
int GPJ__get_datum_params(const struct Key_Value *projinfo, char **datumname, char **params)
Extract the datum transformation-related parameters from a set of general PROJ_INFO parameters...
Definition: proj/datum.c:173
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:350
int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *oproj)
Print projection parameters as used by PROJ.4 for input and output co-ordinate systems.
Definition: get_proj.c:394
#define NULL
Definition: ccmath.h:32
int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a, double *e2, double *rf)
Get the ellipsoid parameters from proj keys structure.
Definition: ellipse.c:73
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
const char * set_proj_lib(const char *name)
Definition: get_proj.c:364
int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys, const struct Key_Value *in_units_keys)
Create a pj_info struct Co-ordinate System definition from a set of PROJ_INFO / PROJ_UNITS-style key-...
Definition: get_proj.c:59
#define MAX_PARGS
Definition: get_proj.c:29
int GPJ_get_default_datum_params_by_name(const char *name, char **params)
"Last resort" function to retrieve a "default" set of datum parameters for a datum (N...
Definition: proj/datum.c:85
#define FINDERFUNC
Definition: get_proj.c:27
const char * name
Definition: named_colr.c:7
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
int pj_get_string(struct pj_info *info, char *str)
Definition: get_proj.c:254
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203