GRASS GIS 7 Programmer's Manual  7.7.svn(2018)-r73570
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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, 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 <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 nopt;
35 
36 /* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
37 
38 /**
39  * \brief Create a pj_info struct Co-ordinate System definition from a set of
40  * PROJ_INFO / PROJ_UNITS-style key-value pairs
41  *
42  * This function takes a GRASS-style co-ordinate system definition as stored
43  * in the PROJ_INFO and PROJ_UNITS files and processes it to create a pj_info
44  * representation for use in re-projecting with pj_do_proj(). In addition to
45  * the parameters passed to it it may also make reference to the system
46  * ellipse.table and datum.table files if necessary.
47  *
48  * \param info Pointer to a pj_info struct (which must already exist) into
49  * which the co-ordinate system definition will be placed
50  * \param in_proj_keys PROJ_INFO-style key-value pairs
51  * \param in_units_keys PROJ_UNITS-style key-value pairs
52  *
53  * \return -1 on error (unable to initialise PROJ.4)
54  * 2 if "default" 3-parameter datum shift values from datum.table
55  * were used
56  * 3 if an unrecognised datum name was passed on to PROJ.4 (and
57  * initialization was successful)
58  * 1 otherwise
59  **/
60 
61 int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
62  const struct Key_Value *in_units_keys)
63 {
64  const char *str;
65  int i;
66  double a, es, rf;
67  int returnval = 1;
68  char buffa[300], factbuff[50];
69  int deflen;
70  char proj_in[250], *datum, *params;
71 #ifdef HAVE_PROJ_H
72  PJ *pj;
73  PJ_CONTEXT *pjc;
74 #else
75  projPJ *pj;
76 #endif
77 
78  proj_in[0] = '\0';
79  info->zone = 0;
80  info->meters = 1.0;
81  info->proj[0] = '\0';
82  info->def = NULL;
83  info->pj = NULL;
84 
85  str = G_find_key_value("meters", in_units_keys);
86  if (str != NULL) {
87  strcpy(factbuff, str);
88  if (strlen(factbuff) > 0)
89  sscanf(factbuff, "%lf", &(info->meters));
90  }
91  str = G_find_key_value("name", in_proj_keys);
92  if (str != NULL) {
93  sprintf(proj_in, "%s", str);
94  }
95  str = G_find_key_value("proj", in_proj_keys);
96  if (str != NULL) {
97  sprintf(info->proj, "%s", str);
98  }
99  if (strlen(info->proj) <= 0)
100  sprintf(info->proj, "ll");
101 
102  nopt = 0;
103  for (i = 0; i < in_proj_keys->nitems; i++) {
104  /* the name parameter is just for grasses use */
105  if (strcmp(in_proj_keys->key[i], "name") == 0) {
106  continue;
107 
108  /* zone handled separately at end of loop */
109  }
110  else if (strcmp(in_proj_keys->key[i], "zone") == 0) {
111  continue;
112 
113  /* Datum and ellipsoid-related parameters will be handled
114  * separately after end of this loop PK */
115 
116  }
117  else if (strcmp(in_proj_keys->key[i], "datum") == 0
118  || strcmp(in_proj_keys->key[i], "dx") == 0
119  || strcmp(in_proj_keys->key[i], "dy") == 0
120  || strcmp(in_proj_keys->key[i], "dz") == 0
121  || strcmp(in_proj_keys->key[i], "datumparams") == 0
122  || strcmp(in_proj_keys->key[i], "nadgrids") == 0
123  || strcmp(in_proj_keys->key[i], "towgs84") == 0
124  || strcmp(in_proj_keys->key[i], "ellps") == 0
125  || strcmp(in_proj_keys->key[i], "a") == 0
126  || strcmp(in_proj_keys->key[i], "b") == 0
127  || strcmp(in_proj_keys->key[i], "es") == 0
128  || strcmp(in_proj_keys->key[i], "f") == 0
129  || strcmp(in_proj_keys->key[i], "rf") == 0) {
130  continue;
131 
132  /* PROJ.4 uses longlat instead of ll as 'projection name' */
133 
134  }
135  else if (strcmp(in_proj_keys->key[i], "proj") == 0) {
136  if (strcmp(in_proj_keys->value[i], "ll") == 0)
137  sprintf(buffa, "proj=longlat");
138  else
139  sprintf(buffa, "proj=%s", in_proj_keys->value[i]);
140 
141  /* 'One-sided' PROJ.4 flags will have the value in
142  * the key-value pair set to 'defined' and only the
143  * key needs to be passed on. */
144  }
145  else if (strcmp(in_proj_keys->value[i], "defined") == 0)
146  sprintf(buffa, "%s", in_proj_keys->key[i]);
147 
148  else
149  sprintf(buffa, "%s=%s",
150  in_proj_keys->key[i], in_proj_keys->value[i]);
151 
152  alloc_options(buffa);
153  }
154 
155  str = G_find_key_value("zone", in_proj_keys);
156  if (str != NULL) {
157  if (sscanf(str, "%d", &(info->zone)) != 1) {
158  G_fatal_error(_("Invalid zone %s specified"), str);
159  }
160  if (info->zone < 0) {
161 
162  /* if zone is negative, write abs(zone) and define south */
163  info->zone = -info->zone;
164 
165  if (G_find_key_value("south", in_proj_keys) == NULL) {
166  sprintf(buffa, "south");
167  alloc_options(buffa);
168  }
169  }
170  sprintf(buffa, "zone=%d", info->zone);
171  alloc_options(buffa);
172  }
173 
174  if ((GPJ__get_ellipsoid_params(in_proj_keys, &a, &es, &rf) == 0)
175  && (str = G_find_key_value("ellps", in_proj_keys)) != NULL) {
176  /* Default values were returned but an ellipsoid name not recognised
177  * by GRASS is present---perhaps it will be recognised by
178  * PROJ.4 even though it wasn't by GRASS */
179  sprintf(buffa, "ellps=%s", str);
180  alloc_options(buffa);
181  }
182  else {
183  sprintf(buffa, "a=%.16g", a);
184  alloc_options(buffa);
185  /* Cannot use es directly because the OSRImportFromProj4()
186  * function in OGR only accepts b or rf as the 2nd parameter */
187  if (es == 0)
188  sprintf(buffa, "b=%.16g", a);
189  else
190  sprintf(buffa, "rf=%.16g", rf);
191  alloc_options(buffa);
192 
193  }
194  /* Workaround to stop PROJ reading values from defaults file when
195  * rf (and sometimes ellps) is not specified */
196  if (G_find_key_value("no_defs", in_proj_keys) == NULL) {
197  sprintf(buffa, "no_defs");
198  alloc_options(buffa);
199  }
200 
201  /* If datum parameters are present in the PROJ_INFO keys, pass them on */
202  if (GPJ__get_datum_params(in_proj_keys, &datum, &params) == 2) {
203  sprintf(buffa, "%s", params);
204  alloc_options(buffa);
205  G_free(params);
206 
207  /* else if a datum name is present take it and look up the parameters
208  * from the datum.table file */
209  }
210  else if (datum != NULL) {
211 
212  if (GPJ_get_default_datum_params_by_name(datum, &params) > 0) {
213  sprintf(buffa, "%s", params);
214  alloc_options(buffa);
215  returnval = 2;
216  G_free(params);
217 
218  /* else just pass the datum name on and hope it is recognised by
219  * PROJ.4 even though it isn't recognised by GRASS */
220  }
221  else {
222  sprintf(buffa, "datum=%s", datum);
223  alloc_options(buffa);
224  returnval = 3;
225  }
226  /* else there'll be no datum transformation taking place here... */
227  }
228  else {
229  returnval = 4;
230  }
231  G_free(datum);
232 
233 #ifdef HAVE_PROJ_H
234  pjc = proj_context_create();
235  if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
236 #else
237  /* Set finder function for locating datum conversion tables PK */
238  pj_set_finder(FINDERFUNC);
239 
240  if (!(pj = pj_init(nopt, opt_in))) {
241 #endif
242  strcpy(buffa,
243  _("Unable to initialise PROJ with the following parameter list:"));
244  for (i = 0; i < nopt; i++) {
245  char err[50];
246 
247  sprintf(err, " +%s", opt_in[i]);
248  strcat(buffa, err);
249  }
250  G_warning("%s", buffa);
251 #ifndef HAVE_PROJ_H
252  G_warning(_("The PROJ error message: %s"), pj_strerrno(pj_errno));
253 #endif
254  return -1;
255  }
256 
257 #ifdef HAVE_PROJ_H
258  int perr = proj_errno(pj);
259 
260  if (perr)
261  G_fatal_error("PROJ 5 error %d", perr);
262 #endif
263 
264  info->pj = pj;
265 
266  deflen = 0;
267  for (i = 0; i < nopt; i++)
268  deflen += strlen(opt_in[i]) + 2;
269 
270  info->def = G_malloc(deflen + 1);
271 
272  sprintf(buffa, "+%s ", opt_in[0]);
273  strcpy(info->def, buffa);
274  G_free(opt_in[0]);
275 
276  for (i = 1; i < nopt; i++) {
277  sprintf(buffa, "+%s ", opt_in[i]);
278  strcat(info->def, buffa);
279  G_free(opt_in[i]);
280  }
281 
282  return returnval;
283 }
284 
285 static void alloc_options(char *buffa)
286 {
287  int nsize;
288 
289  nsize = strlen(buffa);
290  opt_in[nopt++] = (char *)G_malloc(nsize + 1);
291  sprintf(opt_in[nopt - 1], "%s", buffa);
292  return;
293 }
294 
295 /**
296  * \brief Create a pj_info struct Co-ordinate System definition from a
297  * string with a sequence of key=value pairs
298  *
299  * This function takes a GRASS- or PROJ style co-ordinate system definition
300  * and processes it to create a pj_info representation for use in
301  * re-projecting with pj_do_proj(). In addition to the parameters passed
302  * to it it may also make reference to the system ellipse.table and
303  * datum.table files if necessary.
304  *
305  * \param info Pointer to a pj_info struct (which must already exist) into
306  * which the co-ordinate system definition will be placed
307  * \param str input string with projection definition
308  * \param in_units_keys PROJ_UNITS-style key-value pairs
309  *
310  * \return -1 on error (unable to initialise PROJ.4)
311  * 1 on success
312  **/
313 
314 int pj_get_string(struct pj_info *info, char *str)
315 {
316  char *s;
317  int i, nsize;
318  char zonebuff[50], buffa[300];
319  int deflen;
320 #ifdef HAVE_PROJ_H
321  PJ *pj;
322  PJ_CONTEXT *pjc;
323 #else
324  projPJ *pj;
325 #endif
326 
327  info->zone = 0;
328  info->proj[0] = '\0';
329  info->meters = 1.0;
330  info->def = NULL;
331  info->pj = NULL;
332 
333  nopt = 0;
334 
335  if ((str == NULL) || (str[0] == '\0')) {
336  /* Null Pointer or empty string is supplied for parameters,
337  * implying latlong projection; just need to set proj
338  * parameter and call pj_init PK */
339  sprintf(info->proj, "ll");
340  sprintf(buffa, "proj=latlong ellps=WGS84");
341  alloc_options(buffa);
342  }
343  else {
344  /* Parameters have been provided; parse through them but don't
345  * bother with most of the checks in pj_get_kv; assume the
346  * programmer knows what he / she is doing when using this
347  * function rather than reading a PROJ_INFO file PK */
348  s = str;
349  while (s = strtok(s, " \t\n"), s) {
350  if (strncmp(s, "+unfact=", 8) == 0) {
351  s = s + 8;
352  info->meters = atof(s);
353  }
354  else {
355  if (strncmp(s, "+", 1) == 0)
356  ++s;
357  if (nsize = strlen(s), nsize) {
358  if (nopt >= MAX_PARGS) {
359  fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
360  G_fatal_error(_("Option input overflowed option table"));
361  }
362 
363  if (strncmp("zone=", s, 5) == 0) {
364  sprintf(zonebuff, "%s", s + 5);
365  sscanf(zonebuff, "%d", &(info->zone));
366  }
367 
368  if (strncmp("proj=", s, 5) == 0) {
369  sprintf(info->proj, "%s", s + 5);
370  if (strcmp(info->proj, "ll") == 0)
371  sprintf(buffa, "proj=latlong");
372  else
373  sprintf(buffa, "%s", s);
374  }
375  else {
376  sprintf(buffa, "%s", s);
377  }
378  alloc_options(buffa);
379  }
380  }
381  s = 0;
382  }
383  }
384 
385 #ifdef HAVE_PROJ_H
386  pjc = proj_context_create();
387  if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
388  G_warning(_("Unable to initialize pj"));
389  return -1;
390  }
391 #else
392  /* Set finder function for locating datum conversion tables PK */
393  pj_set_finder(FINDERFUNC);
394 
395  if (!(pj = pj_init(nopt, opt_in))) {
396  G_warning(_("Unable to initialize pj cause: %s"),
397  pj_strerrno(pj_errno));
398  return -1;
399  }
400 #endif
401  info->pj = pj;
402 
403  deflen = 0;
404  for (i = 0; i < nopt; i++)
405  deflen += strlen(opt_in[i]) + 2;
406 
407  info->def = G_malloc(deflen + 1);
408 
409  sprintf(buffa, "+%s ", opt_in[0]);
410  strcpy(info->def, buffa);
411  G_free(opt_in[0]);
412 
413  for (i = 1; i < nopt; i++) {
414  sprintf(buffa, "+%s ", opt_in[i]);
415  strcat(info->def, buffa);
416  G_free(opt_in[i]);
417  }
418 
419  return 1;
420 }
421 
422 #ifndef HAVE_PROJ_H
423 /* GPJ_get_equivalent_latlong(): only available with PROJ 4 API
424  * with the new PROJ 5+ API, use pjold directly with PJ_FWD/PJ_INV transformation
425 */
426 /**
427  * \brief Define a latitude / longitude co-ordinate system with the same
428  * ellipsoid and datum parameters as an existing projected system
429  *
430  * This function is useful when projected co-ordinates need to be simply
431  * converted to and from latitude / longitude.
432  *
433  * \param pjnew Pointer to pj_info struct for geographic co-ordinate system
434  * that will be created
435  * \param pjold Pointer to pj_info struct for existing projected co-ordinate
436  * system
437  *
438  * \return 1 on success; -1 if there was an error (i.e. if the PROJ.4
439  * pj_latlong_from_proj() function returned NULL)
440  **/
441 
442 int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
443 {
444  char *deftmp;
445 
446  pjnew->meters = 1.;
447  pjnew->zone = 0;
448  pjnew->def = NULL;
449  sprintf(pjnew->proj, "ll");
450  if ((pjnew->pj = pj_latlong_from_proj(pjold->pj)) == NULL)
451  return -1;
452 
453  deftmp = pj_get_def(pjnew->pj, 1);
454  pjnew->def = G_store(deftmp);
455  pj_dalloc(deftmp);
456 
457  return 1;
458 }
459 #endif
460 
461 /* set_proj_lib()
462  * 'finder function' for use with PROJ.4 pj_set_finder() function */
463 
464 const char *set_proj_lib(const char *name)
465 {
466  const char *gisbase = G_gisbase();
467  static char *buf = NULL;
468  static size_t buf_len;
469  size_t len = strlen(gisbase) + sizeof(GRIDDIR) + strlen(name) + 1;
470 
471  if (buf_len < len) {
472  if (buf != NULL)
473  G_free(buf);
474  buf_len = len + 20;
475  buf = G_malloc(buf_len);
476  }
477 
478  sprintf(buf, "%s%s/%s", gisbase, GRIDDIR, name);
479 
480  return buf;
481 }
482 
483 /**
484  * \brief Print projection parameters as used by PROJ.4 for input and
485  * output co-ordinate systems
486  *
487  * \param iproj 'Input' co-ordinate system
488  * \param oproj 'Output' co-ordinate system
489  *
490  * \return 1 on success, -1 on error (i.e. if the PROJ-style definition
491  * is NULL for either co-ordinate system)
492  **/
493 
494 int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *oproj)
495 {
496  char *str;
497 
498  if (iproj) {
499  str = iproj->def;
500  if (str != NULL) {
501  fprintf(stderr, "%s: %s\n", _("Input Projection Parameters"),
502  str);
503  fprintf(stderr, "%s: %.16g\n", _("Input Unit Factor"),
504  iproj->meters);
505  }
506  else
507  return -1;
508  }
509 
510  if (oproj) {
511  str = oproj->def;
512  if (str != NULL) {
513  fprintf(stderr, "%s: %s\n", _("Output Projection Parameters"),
514  str);
515  fprintf(stderr, "%s: %.16g\n", _("Output Unit Factor"),
516  oproj->meters);
517  }
518  else
519  return -1;
520  }
521 
522  return 1;
523 }
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
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:442
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:494
char proj[100]
Definition: gprojects.h:52
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
int zone
Definition: gprojects.h:51
#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:160
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:464
projPJ pj
Definition: gprojects.h:48
char * def
Definition: gprojects.h:53
char ** value
Definition: gis.h:484
char ** key
Definition: gis.h:483
Definition: gis.h:479
double meters
Definition: gprojects.h:50
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:61
#define MAX_PARGS
Definition: get_proj.c:29
int GPJ_get_default_datum_params_by_name(const char *name, char **params)
&quot;Last resort&quot; function to retrieve a &quot;default&quot; set of datum parameters for a datum (N...
Definition: proj/datum.c:85
#define GRIDDIR
Definition: gprojects.h:40
#define _(str)
Definition: glocale.h:13
#define FINDERFUNC
Definition: get_proj.c:27
int nitems
Definition: gis.h:481
const char * name
Definition: named_colr.c:7
int pj_get_string(struct pj_info *info, char *str)
Create a pj_info struct Co-ordinate System definition from a string with a sequence of key=value pair...
Definition: get_proj.c:314
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:204