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