GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
make_loc.c
Go to the documentation of this file.
1 /*!
2  * \file lib/gis/make_loc.c
3  *
4  * \brief GIS Library - Functions to create a new location
5  *
6  * Creates a new location automatically given a "Cell_head", PROJ_INFO
7  * and PROJ_UNITS information.
8  *
9  * (C) 2000-2013 by the GRASS Development Team
10  *
11  * This program is free software under the GNU General Public License
12  * (>=v2). Read the file COPYING that comes with GRASS for details.
13  *
14  * \author Frank Warmerdam
15  */
16 
17 #include <grass/gis.h>
18 
19 #include <stdlib.h>
20 #include <string.h>
21 #include <errno.h>
22 #include <unistd.h>
23 #include <sys/stat.h>
24 #include <math.h>
25 #include <grass/glocale.h>
26 
27 /*!
28  * \brief Create a new location
29  *
30  * This function creates a new location in the current database,
31  * initializes the projection, default window and current window.
32  *
33  * \param location_name Name of the new location. Should not include
34  * the full path, the location will be created within
35  * the current database.
36  * \param wind default window setting for the new location.
37  * All fields should be set in this
38  * structure, and care should be taken to ensure that
39  * the proj/zone fields match the definition in the
40  * proj_info parameter(see G_set_cellhd_from_projinfo()).
41  *
42  * \param proj_info projection definition suitable to write to the
43  * PROJ_INFO file, or NULL for PROJECTION_XY.
44  *
45  * \param proj_units projection units suitable to write to the PROJ_UNITS
46  * file, or NULL.
47  *
48  * \return 0 on success
49  * \return -1 to indicate a system error (check errno).
50  * \return -2 failed to create projection file (currently not used)
51  * \return -3 illegal name
52  */
53 int G_make_location(const char *location_name,
54  struct Cell_head *wind,
55  const struct Key_Value *proj_info,
56  const struct Key_Value *proj_units)
57 {
58  char path[GPATH_MAX];
59 
60  /* check if location name is legal */
61  if (G_legal_filename(location_name) != 1)
62  return -3;
63 
64  /* Try to create the location directory, under the gisdbase. */
65  sprintf(path, "%s/%s", G_gisdbase(), location_name);
66  if (G_mkdir(path) != 0)
67  return -1;
68 
69  /* Make the PERMANENT mapset. */
70  sprintf(path, "%s/%s/%s", G_gisdbase(), location_name, "PERMANENT");
71  if (G_mkdir(path) != 0) {
72  return -1;
73  }
74 
75  /* make these the new current location and mapset */
76  G_setenv_nogisrc("LOCATION_NAME", location_name);
77  G_setenv_nogisrc("MAPSET", "PERMANENT");
78 
79  /* Create the default, and current window files */
80  G_put_element_window(wind, "", "DEFAULT_WIND");
81  G_put_element_window(wind, "", "WIND");
82 
83  /* Write out the PROJ_INFO, and PROJ_UNITS if available. */
84  if (proj_info != NULL) {
85  G_file_name(path, "", "PROJ_INFO", "PERMANENT");
86  G_write_key_value_file(path, proj_info);
87  }
88 
89  if (proj_units != NULL) {
90  G_file_name(path, "", "PROJ_UNITS", "PERMANENT");
91  G_write_key_value_file(path, proj_units);
92  }
93 
94  return 0;
95 }
96 
97 /*!
98  * \brief Create a new location
99  *
100  * This function creates a new location in the current database,
101  * initializes the projection, default window and current window,
102  * and sets the EPSG code if present
103  *
104  * \param location_name Name of the new location. Should not include
105  * the full path, the location will be created within
106  * the current database.
107  * \param wind default window setting for the new location.
108  * All fields should be set in this
109  * structure, and care should be taken to ensure that
110  * the proj/zone fields match the definition in the
111  * proj_info parameter(see G_set_cellhd_from_projinfo()).
112  *
113  * \param proj_info projection definition suitable to write to the
114  * PROJ_INFO file, or NULL for PROJECTION_XY.
115  *
116  * \param proj_units projection units suitable to write to the PROJ_UNITS
117  * file, or NULL.
118  *
119  * \param proj_epsg EPSG code suitable to write to the PROJ_EPSG
120  * file, or NULL.
121  *
122  * \return 0 on success
123  * \return -1 to indicate a system error (check errno).
124  * \return -2 failed to create projection file (currently not used)
125  * \return -3 illegal name
126  */
127 int G_make_location_epsg(const char *location_name,
128  struct Cell_head *wind,
129  const struct Key_Value *proj_info,
130  const struct Key_Value *proj_units,
131  const struct Key_Value *proj_epsg)
132 {
133  int ret;
134 
135  ret = G_make_location(location_name, wind, proj_info, proj_units);
136 
137  if (ret != 0)
138  return ret;
139 
140  /* Write out the PROJ_EPSG if available. */
141  if (proj_epsg != NULL) {
142  char path[GPATH_MAX];
143  G_file_name(path, "", "PROJ_EPSG", "PERMANENT");
144  G_write_key_value_file(path, proj_epsg);
145  }
146 
147  return 0;
148 }
149 
150 /*!
151  * \brief Create a new location
152  *
153  * This function creates a new location in the current database,
154  * initializes the projection, default window and current window,
155  * and sets WKT, srid, and EPSG code if present
156  *
157  * \param location_name Name of the new location. Should not include
158  * the full path, the location will be created within
159  * the current database.
160  * \param wind default window setting for the new location.
161  * All fields should be set in this
162  * structure, and care should be taken to ensure that
163  * the proj/zone fields match the definition in the
164  * proj_info parameter(see G_set_cellhd_from_projinfo()).
165  *
166  * \param proj_info projection definition suitable to write to the
167  * PROJ_INFO file, or NULL for PROJECTION_XY.
168  *
169  * \param proj_units projection units suitable to write to the PROJ_UNITS
170  * file, or NULL.
171  *
172  * \param proj_epsg EPSG code suitable to write to the PROJ_EPSG
173  * file, or NULL.
174  *
175  * \param proj_wkt WKT definition suitable to write to the PROJ_WKT
176  * file, or NULL.
177  *
178  * \param proj_srid Spatial reference ID suitable to write to the PROJ_SRID
179  * file, or NULL.
180  *
181  * \return 0 on success
182  * \return -1 to indicate a system error (check errno).
183  * \return -2 failed to create projection file (currently not used)
184  * \return -3 illegal name
185  */
186 int G_make_location_crs(const char *location_name,
187  struct Cell_head *wind,
188  const struct Key_Value *proj_info,
189  const struct Key_Value *proj_units,
190  const char *proj_srid,
191  const char *proj_wkt)
192 {
193  int ret;
194 
195  ret = G_make_location(location_name, wind, proj_info, proj_units);
196 
197  if (ret != 0)
198  return ret;
199 
200  /* Write out PROJ_SRID if srid is available. */
201  if (proj_srid != NULL) {
202  G_write_projsrid(location_name, proj_srid);
203  }
204 
205  /* Write out PROJ_WKT if WKT is available. */
206  if (proj_wkt != NULL) {
207  G_write_projwkt(location_name, proj_wkt);
208  }
209 
210  return 0;
211 }
212 
213 /*!
214  * \brief Compare projections including units
215  *
216  * \param proj_info1 projection info to compare
217  * \param proj_units1 projection units to compare
218  * \param proj_info2 projection info to compare
219  * \param proj_units2 projection units to compare
220 
221  * \return -1 if not the same projection
222  * \return -2 if linear unit translation to meters fails
223  * \return -3 if not the same datum,
224  * \return -4 if not the same ellipsoid,
225  * \return -5 if UTM zone differs
226  * \return -6 if UTM hemisphere differs,
227  * \return -7 if false easting differs
228  * \return -8 if false northing differs,
229  * \return -9 if center longitude differs,
230  * \return -10 if center latitude differs,
231  * \return -11 if standard parallels differ,
232  * \return 1 if projections match.
233  */
234 int G_compare_projections(const struct Key_Value *proj_info1,
235  const struct Key_Value *proj_units1,
236  const struct Key_Value *proj_info2,
237  const struct Key_Value *proj_units2)
238 {
239  const char *proj1, *proj2;
240 
241  if (proj_info1 == NULL && proj_info2 == NULL)
242  return TRUE;
243 
244  /* -------------------------------------------------------------------- */
245  /* Are they both in the same projection? */
246  /* -------------------------------------------------------------------- */
247  /* prevent seg fault in G_find_key_value */
248  if (proj_info1 == NULL || proj_info2 == NULL)
249  return -1;
250 
251  proj1 = G_find_key_value("proj", proj_info1);
252  proj2 = G_find_key_value("proj", proj_info2);
253 
254  if (proj1 == NULL || proj2 == NULL || strcmp(proj1, proj2))
255  return -1;
256 
257  /* -------------------------------------------------------------------- */
258  /* Verify that the linear unit translation to meters is OK. */
259  /* -------------------------------------------------------------------- */
260  /* prevent seg fault in G_find_key_value */
261  if (proj_units1 == NULL && proj_units2 == NULL)
262  return 1;
263 
264  if (proj_units1 == NULL || proj_units2 == NULL)
265  return -2;
266 
267  {
268  double a1 = 0, a2 = 0;
269 
270  if (G_find_key_value("meters", proj_units1) != NULL)
271  a1 = atof(G_find_key_value("meters", proj_units1));
272  if (G_find_key_value("meters", proj_units2) != NULL)
273  a2 = atof(G_find_key_value("meters", proj_units2));
274 
275  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
276  return -2;
277  }
278  /* compare unit name only if there is no to meter conversion factor */
279  if (G_find_key_value("meters", proj_units1) == NULL ||
280  G_find_key_value("meters", proj_units2) == NULL) {
281  const char *u_1 = NULL, *u_2 = NULL;
282 
283  u_1 = G_find_key_value("unit", proj_units1);
284  u_2 = G_find_key_value("unit", proj_units2);
285 
286  if ((u_1 && !u_2) || (!u_1 && u_2))
287  return -2;
288 
289  /* the unit name can be arbitrary: the following can be the same
290  * us-ft (proj.4 keyword)
291  * U.S. Surveyor's Foot (proj.4 name)
292  * US survey foot (WKT)
293  * Foot_US (WKT)
294  */
295  if (u_1 && u_2 && G_strcasecmp(u_1, u_2))
296  return -2;
297  }
298 
299  /* -------------------------------------------------------------------- */
300  /* Do they both have the same datum? */
301  /* -------------------------------------------------------------------- */
302  {
303  const char *d_1 = NULL, *d_2 = NULL;
304 
305  d_1 = G_find_key_value("datum", proj_info1);
306  d_2 = G_find_key_value("datum", proj_info2);
307 
308  if ((d_1 && !d_2) || (!d_1 && d_2))
309  return -3;
310 
311  if (d_1 && d_2 && strcmp(d_1, d_2)) {
312  /* different datum short names can mean the same datum,
313  * see lib/gis/datum.table */
314  G_debug(1, "Different datum names");
315  }
316  }
317 
318  /* -------------------------------------------------------------------- */
319  /* Do they both have the same ellipsoid? */
320  /* -------------------------------------------------------------------- */
321  {
322  const char *e_1 = NULL, *e_2 = NULL;
323 
324  e_1 = G_find_key_value("ellps", proj_info1);
325  e_2 = G_find_key_value("ellps", proj_info2);
326 
327  if (e_1 && e_2 && strcmp(e_1, e_2))
328  return -4;
329 
330  if (e_1 == NULL || e_2 == NULL) {
331  double a1 = 0, a2 = 0;
332  double es1 = 0, es2 = 0;
333 
334  /* it may happen that one proj_info has ellps,
335  * while the other has a, es: translate ellps to a, es */
336  if (e_1)
337  G_get_ellipsoid_by_name(e_1, &a1, &es1);
338  else {
339  if (G_find_key_value("a", proj_info1) != NULL)
340  a1 = atof(G_find_key_value("a", proj_info1));
341  if (G_find_key_value("es", proj_info1) != NULL)
342  es1 = atof(G_find_key_value("es", proj_info1));
343  }
344 
345  if (e_2)
346  G_get_ellipsoid_by_name(e_2, &a2, &es2);
347  else {
348  if (G_find_key_value("a", proj_info2) != NULL)
349  a2 = atof(G_find_key_value("a", proj_info2));
350  if (G_find_key_value("es", proj_info2) != NULL)
351  es2 = atof(G_find_key_value("es", proj_info2));
352  }
353 
354  /* it should be an error if a = 0 */
355  if ((a1 == 0 && a2 != 0) || (a1 != 0 && a2 == 0))
356  return -4;
357 
358  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
359  return -4;
360 
361  if ((es1 == 0 && es2 != 0) || (es1 != 0 && es2 == 0))
362  return -4;
363 
364  if (es1 && es2 && (fabs(es2 - es1) > 0.000001))
365  return -4;
366  }
367  }
368 
369  /* -------------------------------------------------------------------- */
370  /* Zone check specially for UTM */
371  /* -------------------------------------------------------------------- */
372  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
373  && atof(G_find_key_value("zone", proj_info1))
374  != atof(G_find_key_value("zone", proj_info2)))
375  return -5;
376 
377  /* -------------------------------------------------------------------- */
378  /* Hemisphere check specially for UTM */
379  /* -------------------------------------------------------------------- */
380  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
381  && !!G_find_key_value("south", proj_info1)
382  != !!G_find_key_value("south", proj_info2))
383  return -6;
384 
385  /* -------------------------------------------------------------------- */
386  /* Do they both have the same false easting? */
387  /* -------------------------------------------------------------------- */
388 
389  {
390  const char *x_0_1 = NULL, *x_0_2 = NULL;
391 
392  x_0_1 = G_find_key_value("x_0", proj_info1);
393  x_0_2 = G_find_key_value("x_0", proj_info2);
394 
395  if ((x_0_1 && !x_0_2) || (!x_0_1 && x_0_2))
396  return -7;
397 
398  if (x_0_1 && x_0_2 && (fabs(atof(x_0_1) - atof(x_0_2)) > 0.000001))
399  return -7;
400  }
401 
402  /* -------------------------------------------------------------------- */
403  /* Do they both have the same false northing? */
404  /* -------------------------------------------------------------------- */
405 
406  {
407  const char *y_0_1 = NULL, *y_0_2 = NULL;
408 
409  y_0_1 = G_find_key_value("y_0", proj_info1);
410  y_0_2 = G_find_key_value("y_0", proj_info2);
411 
412  if ((y_0_1 && !y_0_2) || (!y_0_1 && y_0_2))
413  return -8;
414 
415  if (y_0_1 && y_0_2 && (fabs(atof(y_0_1) - atof(y_0_2)) > 0.000001))
416  return -8;
417  }
418 
419  /* -------------------------------------------------------------------- */
420  /* Do they have the same center longitude? */
421  /* -------------------------------------------------------------------- */
422 
423  {
424  const char *l_1 = NULL, *l_2 = NULL;
425 
426  l_1 = G_find_key_value("lon_0", proj_info1);
427  l_2 = G_find_key_value("lon_0", proj_info2);
428 
429  if ((l_1 && !l_2) || (!l_1 && l_2))
430  return -9;
431 
432  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
433  return -9;
434 
435  /* -------------------------------------------------------------------- */
436  /* Do they have the same center latitude? */
437  /* -------------------------------------------------------------------- */
438 
439  l_1 = G_find_key_value("lat_0", proj_info1);
440  l_2 = G_find_key_value("lat_0", proj_info2);
441 
442  if ((l_1 && !l_2) || (!l_1 && l_2))
443  return -10;
444 
445  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
446  return -10;
447 
448  /* -------------------------------------------------------------------- */
449  /* Do they have the same standard parallels? */
450  /* -------------------------------------------------------------------- */
451 
452  l_1 = G_find_key_value("lat_1", proj_info1);
453  l_2 = G_find_key_value("lat_1", proj_info2);
454 
455  if ((l_1 && !l_2) || (!l_1 && l_2))
456  return -11;
457 
458  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
459  /* lat_1 differ */
460  /* check for swapped lat_1, lat_2 */
461  l_2 = G_find_key_value("lat_2", proj_info2);
462 
463  if (!l_2)
464  return -11;
465  if (fabs(atof(l_1) - atof(l_2)) > 0.000001) {
466  return -11;
467  }
468  }
469 
470  l_1 = G_find_key_value("lat_2", proj_info1);
471  l_2 = G_find_key_value("lat_2", proj_info2);
472 
473  if ((l_1 && !l_2) || (!l_1 && l_2))
474  return -11;
475 
476  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
477  /* lat_2 differ */
478  /* check for swapped lat_1, lat_2 */
479  l_2 = G_find_key_value("lat_1", proj_info2);
480 
481  if (!l_2)
482  return -11;
483  if (fabs(atof(l_1) - atof(l_2)) > 0.000001) {
484  return -11;
485  }
486  }
487  }
488 
489  /* towgs84 ? */
490 
491  /* -------------------------------------------------------------------- */
492  /* Add more details in later. */
493  /* -------------------------------------------------------------------- */
494 
495  return 1;
496 }
497 
498 /*!
499  \brief Write WKT definition to file
500 
501  Any WKT string and version recognized by PROJ is supported.
502 
503  \param location_name name of the location to write the WKT definition
504  \param wktstring pointer to WKT string
505 
506  \return 0 success
507  \return -1 error writing
508  */
509 
510 int G_write_projwkt(const char *location_name, const char *wktstring)
511 {
512  FILE *fp;
513  char path[GPATH_MAX];
514  int err, n;
515 
516  if (!wktstring)
517  return 0;
518 
519  if (location_name && *location_name)
520  sprintf(path, "%s/%s/%s/%s", G_gisdbase(), location_name, "PERMANENT", WKT_FILE);
521  else
522  G_file_name(path, "", WKT_FILE, "PERMANENT");
523 
524  fp = fopen(path, "w");
525 
526  if (!fp)
527  G_fatal_error(_("Unable to open output file <%s>: %s"), path, strerror(errno));
528 
529  err = 0;
530  n = strlen(wktstring);
531  if (wktstring[n - 1] != '\n') {
532  if (n != fprintf(fp, "%s\n", wktstring))
533  err = -1;
534  }
535  else {
536  if (n != fprintf(fp, "%s", wktstring))
537  err = -1;
538  }
539 
540  if (fclose(fp) != 0)
541  G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
542 
543  return err;
544 }
545 
546 
547 /*!
548  \brief Write srid (spatial reference id) to file
549 
550  A srid consists of an authority name and code and must be known to
551  PROJ.
552 
553  \param location_name name of the location to write the srid
554  \param sridstring pointer to srid string
555 
556  \return 0 success
557  \return -1 error writing
558  */
559 
560 int G_write_projsrid(const char *location_name, const char *sridstring)
561 {
562  FILE *fp;
563  char path[GPATH_MAX];
564  int err, n;
565 
566  if (!sridstring)
567  return 0;
568 
569  if (location_name && *location_name)
570  sprintf(path, "%s/%s/%s/%s", G_gisdbase(), location_name, "PERMANENT", SRID_FILE);
571  else
572  G_file_name(path, "", SRID_FILE, "PERMANENT");
573 
574  fp = fopen(path, "w");
575 
576  if (!fp)
577  G_fatal_error(_("Unable to open output file <%s>: %s"), path, strerror(errno));
578 
579  err = 0;
580  n = strlen(sridstring);
581  if (sridstring[n - 1] != '\n') {
582  if (n != fprintf(fp, "%s\n", sridstring))
583  err = -1;
584  }
585  else {
586  if (n != fprintf(fp, "%s", sridstring))
587  err = -1;
588  }
589 
590  if (fclose(fp) != 0)
591  G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
592 
593  return err;
594 }
#define TRUE
Definition: gis.h:59
char * G_file_name(char *, const char *, const char *, const char *)
Builds full path names to GIS data files.
Definition: file_name.c:61
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_setenv_nogisrc(const char *, const char *)
Set environment name to value (doesn&#39;t update .gisrc)
Definition: env.c:465
2D/3D raster map header (used also for region)
Definition: gis.h:423
const char * G_gisdbase(void)
Get name of top level database directory.
Definition: gisdbase.c:26
#define NULL
Definition: ccmath.h:32
int G_compare_projections(const struct Key_Value *proj_info1, const struct Key_Value *proj_units1, const struct Key_Value *proj_info2, const struct Key_Value *proj_units2)
Compare projections including units.
Definition: make_loc.c:234
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
int int G_strcasecmp(const char *, const char *)
String compare ignoring case (upper or lower)
Definition: strings.c:47
int G_get_ellipsoid_by_name(const char *, double *, double *)
Get ellipsoid parameters by name.
Definition: get_ellipse.c:105
int G_legal_filename(const char *)
Check for legal database file name.
Definition: legal_name.c:34
#define GPATH_MAX
Definition: gis.h:180
int G_make_location_epsg(const char *location_name, struct Cell_head *wind, const struct Key_Value *proj_info, const struct Key_Value *proj_units, const struct Key_Value *proj_epsg)
Create a new location.
Definition: make_loc.c:127
int G_write_projsrid(const char *location_name, const char *sridstring)
Write srid (spatial reference id) to file.
Definition: make_loc.c:560
int G_make_location_crs(const char *location_name, struct Cell_head *wind, const struct Key_Value *proj_info, const struct Key_Value *proj_units, const char *proj_srid, const char *proj_wkt)
Create a new location.
Definition: make_loc.c:186
int G_put_element_window(const struct Cell_head *, const char *, const char *)
Write the region.
Definition: put_window.c:74
void G_write_key_value_file(const char *, const struct Key_Value *)
Write key/value pairs to file.
Definition: key_value3.c:28
Definition: gis.h:512
#define SRID_FILE
Definition: gis.h:124
Definition: path.h:16
int G_mkdir(const char *)
Creates a new directory.
Definition: paths.c:27
#define WKT_FILE
Definition: gis.h:123
#define _(str)
Definition: glocale.h:10
int G_make_location(const char *location_name, struct Cell_head *wind, const struct Key_Value *proj_info, const struct Key_Value *proj_units)
Create a new location.
Definition: make_loc.c:53
int G_debug(int, 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:84
int G_write_projwkt(const char *location_name, const char *wktstring)
Write WKT definition to file.
Definition: make_loc.c:510