GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
30static void alloc_options(char *);
31
32static char *opt_in[MAX_PARGS];
33static 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
60int 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
310static 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
339int 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
493int 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
518const 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
554int 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}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#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:74
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
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:493
#define FINDERFUNC
Definition get_proj.c:26
const char * set_proj_share(const char *name)
Definition get_proj.c:518
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
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition key_value1.c:85
const char * name
Definition named_colr.c:6
#define strcpy
Definition parson.c:62
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:86
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)