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