GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-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
31static void alloc_options(char *);
32
33static char *opt_in[MAX_PARGS];
34static 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
61int 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
302static 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
331int 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
471int 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
496const 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
532int 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}
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
#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
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 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
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:471
#define FINDERFUNC
Definition: get_proj.c:27
const char * set_proj_share(const char *name)
Definition: get_proj.c:496
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
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
const char * name
Definition: named_colr.c:7
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
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
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:87
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220