GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-df949ca623
get_proj.c
Go to the documentation of this file.
1 /**
2  \file get_proj.c
3 
4  \brief GProj library - Functions for re-projecting point data
5 
6  \author Original Author unknown, probably Soil Conservation Service,
7  Eric Miller, Paul Kelly, Markus Metz
8 
9  (C) 2003-2008, 2018 by the GRASS Development Team
10 
11  This program is free software under the GNU General Public
12  License (>=v2). Read the file COPYING that comes with GRASS
13  for details.
14 **/
15 
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <ctype.h>
19 #include <math.h>
20 #include <string.h>
21 #include <grass/gis.h>
22 #include <grass/gprojects.h>
23 #include <grass/glocale.h>
24 
25 /* Finder function for datum transformation grids */
26 #define FINDERFUNC set_proj_share
27 #define PERMANENT "PERMANENT"
28 #define MAX_PARGS 100
29 
30 static void alloc_options(char *);
31 
32 static char *opt_in[MAX_PARGS];
33 static int nopt;
34 
35 /* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
36 
37 /**
38  * \brief Create a pj_info struct Co-ordinate System definition from a set of
39  * PROJ_INFO / PROJ_UNITS-style key-value pairs
40  *
41  * This function takes a GRASS-style co-ordinate system definition as stored
42  * in the PROJ_INFO and PROJ_UNITS files and processes it to create a pj_info
43  * representation for use in re-projecting with pj_do_proj(). In addition to
44  * the parameters passed to it it may also make reference to the system
45  * ellipse.table and datum.table files if necessary.
46  *
47  * \param info Pointer to a pj_info struct (which must already exist) into
48  * which the co-ordinate system definition will be placed
49  * \param in_proj_keys PROJ_INFO-style key-value pairs
50  * \param in_units_keys PROJ_UNITS-style key-value pairs
51  *
52  * \return -1 on error (unable to initialise PROJ.4)
53  * 2 if "default" 3-parameter datum shift values from datum.table
54  * were used
55  * 3 if an unrecognised datum name was passed on to PROJ.4 (and
56  * initialization was successful)
57  * 1 otherwise
58  **/
59 
60 int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
61  const struct Key_Value *in_units_keys)
62 {
63  const char *str;
64  int i;
65  double a, es, rf;
66  int returnval = 1;
67  char buffa[300], factbuff[50];
68  int deflen;
69  char proj_in[250], *datum, *params;
70 
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  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  else if (strcmp(in_proj_keys->key[i], "proj") == 0) {
145  if (strcmp(in_proj_keys->value[i], "ll") == 0)
146  sprintf(buffa, "proj=longlat");
147  else
148  sprintf(buffa, "proj=%s", in_proj_keys->value[i]);
149 
150  /* 'One-sided' PROJ.4 flags will have the value in
151  * the key-value pair set to 'defined' and only the
152  * key needs to be passed on. */
153  }
154  else if (strcmp(in_proj_keys->value[i], "defined") == 0)
155  sprintf(buffa, "%s", in_proj_keys->key[i]);
156 
157  else
158  sprintf(buffa, "%s=%s", in_proj_keys->key[i],
159  in_proj_keys->value[i]);
160 
161  alloc_options(buffa);
162  }
163 
164  str = G_find_key_value("zone", in_proj_keys);
165  if (str != NULL) {
166  if (sscanf(str, "%d", &(info->zone)) != 1) {
167  G_fatal_error(_("Invalid zone %s specified"), str);
168  }
169  if (info->zone < 0) {
170 
171  /* if zone is negative, write abs(zone) and define south */
172  info->zone = -info->zone;
173 
174  if (G_find_key_value("south", in_proj_keys) == NULL) {
175  sprintf(buffa, "south");
176  alloc_options(buffa);
177  }
178  }
179  sprintf(buffa, "zone=%d", info->zone);
180  alloc_options(buffa);
181  }
182 
183  if ((GPJ__get_ellipsoid_params(in_proj_keys, &a, &es, &rf) == 0) &&
184  (str = G_find_key_value("ellps", in_proj_keys)) != NULL) {
185  /* Default values were returned but an ellipsoid name not recognised
186  * by GRASS is present---perhaps it will be recognised by
187  * PROJ.4 even though it wasn't by GRASS */
188  sprintf(buffa, "ellps=%s", str);
189  alloc_options(buffa);
190  }
191  else {
192  sprintf(buffa, "a=%.16g", a);
193  alloc_options(buffa);
194  /* Cannot use es directly because the OSRImportFromProj4()
195  * function in OGR only accepts b or rf as the 2nd parameter */
196  if (es == 0)
197  sprintf(buffa, "b=%.16g", a);
198  else
199  sprintf(buffa, "rf=%.16g", rf);
200  alloc_options(buffa);
201  }
202  /* Workaround to stop PROJ reading values from defaults file when
203  * rf (and sometimes ellps) is not specified */
204  if (G_find_key_value("no_defs", in_proj_keys) == NULL) {
205  sprintf(buffa, "no_defs");
206  alloc_options(buffa);
207  }
208 
209  /* If datum parameters are present in the PROJ_INFO keys, pass them on */
210  if (GPJ__get_datum_params(in_proj_keys, &datum, &params) == 2) {
211  sprintf(buffa, "%s", params);
212  alloc_options(buffa);
213  G_free(params);
214 
215  /* else if a datum name is present take it and look up the parameters
216  * from the datum.table file */
217  }
218  else if (datum != NULL) {
219 
220  if (GPJ_get_default_datum_params_by_name(datum, &params) > 0) {
221  sprintf(buffa, "%s", params);
222  alloc_options(buffa);
223  returnval = 2;
224  G_free(params);
225 
226  /* else just pass the datum name on and hope it is recognised by
227  * PROJ.4 even though it isn't recognised by GRASS */
228  }
229  else {
230  sprintf(buffa, "datum=%s", datum);
231  alloc_options(buffa);
232  returnval = 3;
233  }
234  /* else there'll be no datum transformation taking place here... */
235  }
236  else {
237  returnval = 4;
238  }
239  G_free(datum);
240 
241 #ifdef HAVE_PROJ_H
242 #if PROJ_VERSION_MAJOR >= 6
243  /* without type=crs, PROJ6 does not recognize what this is,
244  * a crs or some kind of coordinate operation, falling through to
245  * PJ_TYPE_OTHER_COORDINATE_OPERATION */
246  alloc_options("type=crs");
247 #endif
248  pjc = proj_context_create();
249  if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
250 #else
251  /* Set finder function for locating datum conversion tables PK */
252  pj_set_finder(FINDERFUNC);
253 
254  if (!(pj = pj_init(nopt, opt_in))) {
255 #endif
256  strcpy(
257  buffa,
258  _("Unable to initialise PROJ with the following parameter list:"));
259  for (i = 0; i < nopt; i++) {
260  char err[50];
261 
262  sprintf(err, " +%s", opt_in[i]);
263  strcat(buffa, err);
264  }
265  G_warning("%s", buffa);
266 #ifndef HAVE_PROJ_H
267  G_warning(_("The PROJ error message: %s"), pj_strerrno(pj_errno));
268 #endif
269  return -1;
270  }
271 
272 #ifdef HAVE_PROJ_H
273  int perr = proj_errno(pj);
274 
275  if (perr)
276  G_fatal_error("PROJ 5 error %d", perr);
277 
278 #if PROJ_VERSION_MAJOR >= 6
279  if (proj_get_type(pj) == PJ_TYPE_BOUND_CRS) {
280  PJ *source_crs = proj_get_source_crs(pjc, pj);
281  if (source_crs) {
282  proj_destroy(pj);
283  pj = source_crs;
284  }
285  }
286 #endif
287 #endif
288 
289  info->pj = pj;
290 
291  deflen = 0;
292  for (i = 0; i < nopt; i++)
293  deflen += strlen(opt_in[i]) + 2;
294 
295  info->def = G_malloc(deflen + 1);
296 
297  sprintf(buffa, "+%s ", opt_in[0]);
298  strcpy(info->def, buffa);
299  G_free(opt_in[0]);
300 
301  for (i = 1; i < nopt; i++) {
302  sprintf(buffa, "+%s ", opt_in[i]);
303  strcat(info->def, buffa);
304  G_free(opt_in[i]);
305  }
306 
307  return returnval;
308 }
309 
310 static void alloc_options(char *buffa)
311 {
312  int nsize;
313 
314  nsize = strlen(buffa);
315  opt_in[nopt++] = (char *)G_malloc(nsize + 1);
316  sprintf(opt_in[nopt - 1], "%s", buffa);
317  return;
318 }
319 
320 /**
321  * \brief Create a pj_info struct Co-ordinate System definition from a
322  * string with a sequence of key=value pairs
323  *
324  * This function takes a GRASS- or PROJ style co-ordinate system definition
325  * and processes it to create a pj_info representation for use in
326  * re-projecting with pj_do_proj(). In addition to the parameters passed
327  * to it it may also make reference to the system ellipse.table and
328  * datum.table files if necessary.
329  *
330  * \param info Pointer to a pj_info struct (which must already exist) into
331  * which the co-ordinate system definition will be placed
332  * \param str input string with projection definition
333  * \param in_units_keys PROJ_UNITS-style key-value pairs
334  *
335  * \return -1 on error (unable to initialise PROJ.4)
336  * 1 on success
337  **/
338 
339 int pj_get_string(struct pj_info *info, char *str)
340 {
341  char *s;
342  int i, nsize;
343  char zonebuff[50], buffa[300];
344  int deflen;
345 
346 #ifdef HAVE_PROJ_H
347  PJ *pj;
348  PJ_CONTEXT *pjc;
349 #else
350  projPJ *pj;
351 #endif
352 
353  info->zone = 0;
354  info->proj[0] = '\0';
355  info->meters = 1.0;
356  info->def = NULL;
357  info->srid = NULL;
358  info->pj = NULL;
359 
360  nopt = 0;
361 
362  if ((str == NULL) || (str[0] == '\0')) {
363  /* Null Pointer or empty string is supplied for parameters,
364  * implying latlong projection; just need to set proj
365  * parameter and call pj_init PK */
366  sprintf(info->proj, "ll");
367  sprintf(buffa, "proj=latlong ellps=WGS84");
368  alloc_options(buffa);
369  }
370  else {
371  /* Parameters have been provided; parse through them but don't
372  * bother with most of the checks in pj_get_kv; assume the
373  * programmer knows what he / she is doing when using this
374  * function rather than reading a PROJ_INFO file PK */
375  s = str;
376  while (s = strtok(s, " \t\n"), s) {
377  if (strncmp(s, "+unfact=", 8) == 0) {
378  s = s + 8;
379  info->meters = atof(s);
380  }
381  else {
382  if (strncmp(s, "+", 1) == 0)
383  ++s;
384  if (nsize = strlen(s), nsize) {
385  if (nopt >= MAX_PARGS) {
386  fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
388  _("Option input overflowed option table"));
389  }
390 
391  if (strncmp("zone=", s, 5) == 0) {
392  sprintf(zonebuff, "%s", s + 5);
393  sscanf(zonebuff, "%d", &(info->zone));
394  }
395 
396  if (strncmp(s, "init=", 5) == 0) {
397  info->srid = G_store(s + 6);
398  }
399 
400  if (strncmp("proj=", s, 5) == 0) {
401  sprintf(info->proj, "%s", s + 5);
402  if (strcmp(info->proj, "ll") == 0)
403  sprintf(buffa, "proj=latlong");
404  else
405  sprintf(buffa, "%s", s);
406  }
407  else {
408  sprintf(buffa, "%s", s);
409  }
410  alloc_options(buffa);
411  }
412  }
413  s = 0;
414  }
415  }
416 
417 #ifdef HAVE_PROJ_H
418 #if PROJ_VERSION_MAJOR >= 6
419  /* without type=crs, PROJ6 does not recognize what this is,
420  * a crs or some kind of coordinate operation, falling through to
421  * PJ_TYPE_OTHER_COORDINATE_OPERATION */
422  alloc_options("type=crs");
423 #endif
424  pjc = proj_context_create();
425  if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
426  G_warning(_("Unable to initialize pj cause: %s"),
427  proj_errno_string(proj_context_errno(pjc)));
428  return -1;
429  }
430 
431 #if PROJ_VERSION_MAJOR >= 6
432  if (proj_get_type(pj) == PJ_TYPE_BOUND_CRS) {
433  PJ *source_crs = proj_get_source_crs(pjc, pj);
434  if (source_crs) {
435  proj_destroy(pj);
436  pj = source_crs;
437  }
438  }
439 #endif
440 #else
441  /* Set finder function for locating datum conversion tables PK */
442  pj_set_finder(FINDERFUNC);
443 
444  if (!(pj = pj_init(nopt, opt_in))) {
445  G_warning(_("Unable to initialize pj cause: %s"),
446  pj_strerrno(pj_errno));
447  return -1;
448  }
449 #endif
450  info->pj = pj;
451 
452  deflen = 0;
453  for (i = 0; i < nopt; i++)
454  deflen += strlen(opt_in[i]) + 2;
455 
456  info->def = G_malloc(deflen + 1);
457 
458  sprintf(buffa, "+%s ", opt_in[0]);
459  strcpy(info->def, buffa);
460  G_free(opt_in[0]);
461 
462  for (i = 1; i < nopt; i++) {
463  sprintf(buffa, "+%s ", opt_in[i]);
464  strcat(info->def, buffa);
465  G_free(opt_in[i]);
466  }
467 
468  return 1;
469 }
470 
471 #ifndef HAVE_PROJ_H
472 /* GPJ_get_equivalent_latlong(): only available with PROJ 4 API
473  * with the new PROJ 5+ API, use pjold directly with PJ_FWD/PJ_INV
474  * transformation
475  */
476 
477 /**
478  * \brief Define a latitude / longitude co-ordinate system with the same
479  * ellipsoid and datum parameters as an existing projected system
480  *
481  * This function is useful when projected co-ordinates need to be simply
482  * converted to and from latitude / longitude.
483  *
484  * \param pjnew Pointer to pj_info struct for geographic co-ordinate system
485  * that will be created
486  * \param pjold Pointer to pj_info struct for existing projected co-ordinate
487  * system
488  *
489  * \return 1 on success; -1 if there was an error (i.e. if the PROJ.4
490  * pj_latlong_from_proj() function returned NULL)
491  **/
492 
493 int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
494 {
495  char *deftmp;
496 
497  pjnew->meters = 1.;
498  pjnew->zone = 0;
499  pjnew->def = NULL;
500  sprintf(pjnew->proj, "ll");
501  if ((pjnew->pj = pj_latlong_from_proj(pjold->pj)) == NULL)
502  return -1;
503 
504  deftmp = pj_get_def(pjnew->pj, 1);
505  pjnew->def = G_store(deftmp);
506  pj_dalloc(deftmp);
507 
508  return 1;
509 }
510 #endif
511 
512 /* set_proj_share()
513  * 'finder function' for use with PROJ.4 pj_set_finder() function
514  * this is used to find grids, usually in /usr/share/proj
515  * GRASS no longer provides copies of proj grids in GRIDDIR
516  * -> do not use gisbase/GRIDDIR */
517 
518 const char *set_proj_share(const char *name)
519 {
520  static char *buf = NULL;
521  const char *projshare;
522  static size_t buf_len = 0;
523  size_t len;
524 
525  projshare = getenv("GRASS_PROJSHARE");
526  if (!projshare)
527  return NULL;
528 
529  len = strlen(projshare) + strlen(name) + 2;
530 
531  if (buf_len < len) {
532  if (buf != NULL)
533  G_free(buf);
534  buf_len = len + 20;
535  buf = G_malloc(buf_len);
536  }
537 
538  sprintf(buf, "%s/%s", projshare, name);
539 
540  return buf;
541 }
542 
543 /**
544  * \brief Print projection parameters as used by PROJ.4 for input and
545  * output co-ordinate systems
546  *
547  * \param iproj 'Input' co-ordinate system
548  * \param oproj 'Output' co-ordinate system
549  *
550  * \return 1 on success, -1 on error (i.e. if the PROJ-style definition
551  * is NULL for either co-ordinate system)
552  **/
553 
554 int pj_print_proj_params(const struct pj_info *iproj,
555  const struct pj_info *oproj)
556 {
557  char *str;
558 
559  if (iproj) {
560  str = iproj->def;
561  if (str != NULL) {
562  fprintf(stderr, "%s: %s\n", _("Input Projection Parameters"), str);
563  fprintf(stderr, "%s: %.16g\n", _("Input Unit Factor"),
564  iproj->meters);
565  }
566  else
567  return -1;
568  }
569 
570  if (oproj) {
571  str = oproj->def;
572  if (str != NULL) {
573  fprintf(stderr, "%s: %s\n", _("Output Projection Parameters"), str);
574  fprintf(stderr, "%s: %.16g\n", _("Output Unit Factor"),
575  oproj->meters);
576  }
577  else
578  return -1;
579  }
580 
581  return 1;
582 }
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
const char * G_find_key_value(const char *, const struct Key_Value *)
Find given key (case sensitive)
Definition: key_value1.c:85
#define G_malloc(n)
Definition: defs/gis.h:94
char * G_store(const char *)
Copy string to allocated memory.
Definition: strings.c:87
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:86
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
int GPJ__get_ellipsoid_params(const struct Key_Value *, double *, double *, double *)
Get the ellipsoid parameters from proj keys structure.
Definition: ellipse.c:74
int GPJ_get_equivalent_latlong(struct pj_info *, struct pj_info *)
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:60
#define MAX_PARGS
Definition: get_proj.c:28
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:554
const char * set_proj_share(const char *name)
Definition: get_proj.c:518
#define FINDERFUNC
Definition: get_proj.c:26
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:339
#define _(str)
Definition: glocale.h:10
const char * name
Definition: named_colr.c:6
#define strcpy
Definition: parson.c:62
Definition: gis.h:527
char ** value
Definition: gis.h:531
int nitems
Definition: gis.h:528
char ** key
Definition: gis.h:530
char proj[100]
Definition: gprojects.h:79
int zone
Definition: gprojects.h:78
char * def
Definition: gprojects.h:80
PJ * pj
Definition: gprojects.h:73
char * srid
Definition: gprojects.h:81
char * wkt
Definition: gprojects.h:82
double meters
Definition: gprojects.h:77
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:216