GRASS GIS 8 Programmer's Manual  8.5.0dev(2025)-eb16a2cfc9
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, struct Cell_head *wind,
54  const struct Key_Value *proj_info,
55  const struct Key_Value *proj_units)
56 {
57  char path[GPATH_MAX];
58 
59  /* check if location name is legal */
60  if (G_legal_filename(location_name) != 1)
61  return -3;
62 
63  /* Try to create the location directory, under the gisdbase. */
64  snprintf(path, sizeof(path), "%s/%s", G_gisdbase(), location_name);
65  if (G_mkdir(path) != 0)
66  return -1;
67 
68  /* Make the PERMANENT mapset. */
69  snprintf(path, sizeof(path), "%s/%s/%s", G_gisdbase(), location_name,
70  "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, struct Cell_head *wind,
128  const struct Key_Value *proj_info,
129  const struct Key_Value *proj_units,
130  const struct Key_Value *proj_epsg)
131 {
132  int ret;
133 
134  ret = G_make_location(location_name, wind, proj_info, proj_units);
135 
136  if (ret != 0)
137  return ret;
138 
139  /* Write out the PROJ_EPSG if available. */
140  if (proj_epsg != NULL) {
141  char path[GPATH_MAX];
142 
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, struct Cell_head *wind,
187  const struct Key_Value *proj_info,
188  const struct Key_Value *proj_units,
189  const char *proj_srid, const char *proj_wkt)
190 {
191  int ret;
192 
193  ret = G_make_location(location_name, wind, proj_info, proj_units);
194 
195  if (ret != 0)
196  return ret;
197 
198  /* Write out PROJ_SRID if srid is available. */
199  if (proj_srid != NULL) {
200  G_write_projsrid(location_name, proj_srid);
201  }
202 
203  /* Write out PROJ_WKT if WKT is available. */
204  if (proj_wkt != NULL) {
205  G_write_projwkt(location_name, proj_wkt);
206  }
207 
208  return 0;
209 }
210 
211 /*!
212  * \brief Compare projections including units
213  *
214  * \param proj_info1 projection info to compare
215  * \param proj_units1 projection units to compare
216  * \param proj_info2 projection info to compare
217  * \param proj_units2 projection units to compare
218 
219  * \return -1 if not the same projection
220  * \return -2 if linear unit translation to meters fails
221  * \return -3 if not the same datum,
222  * \return -4 if not the same ellipsoid,
223  * \return -5 if UTM zone differs
224  * \return -6 if UTM hemisphere differs,
225  * \return -7 if false easting differs
226  * \return -8 if false northing differs,
227  * \return -9 if center longitude differs,
228  * \return -10 if center latitude differs,
229  * \return -11 if standard parallels differ,
230  * \return 1 if projections match.
231  */
232 int G_compare_projections(const struct Key_Value *proj_info1,
233  const struct Key_Value *proj_units1,
234  const struct Key_Value *proj_info2,
235  const struct Key_Value *proj_units2)
236 {
237  const char *proj1, *proj2;
238 
239  if (proj_info1 == NULL && proj_info2 == NULL)
240  return TRUE;
241 
242  /* -------------------------------------------------------------------- */
243  /* Are they both in the same projection? */
244  /* -------------------------------------------------------------------- */
245  /* prevent seg fault in G_find_key_value */
246  if (proj_info1 == NULL || proj_info2 == NULL)
247  return -1;
248 
249  proj1 = G_find_key_value("proj", proj_info1);
250  proj2 = G_find_key_value("proj", proj_info2);
251 
252  if (proj1 == NULL || proj2 == NULL || strcmp(proj1, proj2))
253  return -1;
254 
255  /* -------------------------------------------------------------------- */
256  /* Verify that the linear unit translation to meters is OK. */
257  /* -------------------------------------------------------------------- */
258  /* prevent seg fault in G_find_key_value */
259  if (proj_units1 == NULL && proj_units2 == NULL)
260  return 1;
261 
262  if (proj_units1 == NULL || proj_units2 == NULL)
263  return -2;
264 
265  {
266  double a1 = 0, a2 = 0;
267 
268  if (G_find_key_value("meters", proj_units1) != NULL)
269  a1 = atof(G_find_key_value("meters", proj_units1));
270  if (G_find_key_value("meters", proj_units2) != NULL)
271  a2 = atof(G_find_key_value("meters", proj_units2));
272 
273  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
274  return -2;
275  }
276  /* compare unit name only if there is no to meter conversion factor */
277  if (G_find_key_value("meters", proj_units1) == NULL ||
278  G_find_key_value("meters", proj_units2) == NULL) {
279  const char *u_1 = NULL, *u_2 = NULL;
280 
281  u_1 = G_find_key_value("unit", proj_units1);
282  u_2 = G_find_key_value("unit", proj_units2);
283 
284  if ((u_1 && !u_2) || (!u_1 && u_2))
285  return -2;
286 
287  /* the unit name can be arbitrary: the following can be the same
288  * us-ft (proj.4 keyword)
289  * U.S. Surveyor's Foot (proj.4 name)
290  * US survey foot (WKT)
291  * Foot_US (WKT)
292  */
293  if (u_1 && u_2 && G_strcasecmp(u_1, u_2))
294  return -2;
295  }
296 
297  /* -------------------------------------------------------------------- */
298  /* Do they both have the same datum? */
299  /* -------------------------------------------------------------------- */
300  {
301  const char *d_1 = NULL, *d_2 = NULL;
302 
303  d_1 = G_find_key_value("datum", proj_info1);
304  d_2 = G_find_key_value("datum", proj_info2);
305 
306  if ((d_1 && !d_2) || (!d_1 && d_2))
307  return -3;
308 
309  if (d_1 && d_2 && strcmp(d_1, d_2)) {
310  /* different datum short names can mean the same datum,
311  * see lib/gis/datum.table */
312  G_debug(1, "Different datum names");
313  }
314  }
315 
316  /* -------------------------------------------------------------------- */
317  /* Do they both have the same ellipsoid? */
318  /* -------------------------------------------------------------------- */
319  {
320  const char *e_1 = NULL, *e_2 = NULL;
321 
322  e_1 = G_find_key_value("ellps", proj_info1);
323  e_2 = G_find_key_value("ellps", proj_info2);
324 
325  if (e_1 && e_2 && strcmp(e_1, e_2))
326  return -4;
327 
328  if (e_1 == NULL || e_2 == NULL) {
329  double a1 = 0, a2 = 0;
330  double es1 = 0, es2 = 0;
331 
332  /* it may happen that one proj_info has ellps,
333  * while the other has a, es: translate ellps to a, es */
334  if (e_1)
335  G_get_ellipsoid_by_name(e_1, &a1, &es1);
336  else {
337  if (G_find_key_value("a", proj_info1) != NULL)
338  a1 = atof(G_find_key_value("a", proj_info1));
339  if (G_find_key_value("es", proj_info1) != NULL)
340  es1 = atof(G_find_key_value("es", proj_info1));
341  }
342 
343  if (e_2)
344  G_get_ellipsoid_by_name(e_2, &a2, &es2);
345  else {
346  if (G_find_key_value("a", proj_info2) != NULL)
347  a2 = atof(G_find_key_value("a", proj_info2));
348  if (G_find_key_value("es", proj_info2) != NULL)
349  es2 = atof(G_find_key_value("es", proj_info2));
350  }
351 
352  /* it should be an error if a = 0 */
353  if ((a1 == 0 && a2 != 0) || (a1 != 0 && a2 == 0))
354  return -4;
355 
356  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
357  return -4;
358 
359  if ((es1 == 0 && es2 != 0) || (es1 != 0 && es2 == 0))
360  return -4;
361 
362  if (es1 && es2 && (fabs(es2 - es1) > 0.000001))
363  return -4;
364  }
365  }
366 
367  /* -------------------------------------------------------------------- */
368  /* Zone check specially for UTM */
369  /* -------------------------------------------------------------------- */
370  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm") &&
371  atof(G_find_key_value("zone", proj_info1)) !=
372  atof(G_find_key_value("zone", proj_info2)))
373  return -5;
374 
375  /* -------------------------------------------------------------------- */
376  /* Hemisphere check specially for UTM */
377  /* -------------------------------------------------------------------- */
378  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm") &&
379  !!G_find_key_value("south", proj_info1) !=
380  !!G_find_key_value("south", proj_info2))
381  return -6;
382 
383  /* -------------------------------------------------------------------- */
384  /* Do they both have the same false easting? */
385  /* -------------------------------------------------------------------- */
386 
387  {
388  const char *x_0_1 = NULL, *x_0_2 = NULL;
389 
390  x_0_1 = G_find_key_value("x_0", proj_info1);
391  x_0_2 = G_find_key_value("x_0", proj_info2);
392 
393  if ((x_0_1 && !x_0_2) || (!x_0_1 && x_0_2))
394  return -7;
395 
396  if (x_0_1 && x_0_2 && (fabs(atof(x_0_1) - atof(x_0_2)) > 0.000001))
397  return -7;
398  }
399 
400  /* -------------------------------------------------------------------- */
401  /* Do they both have the same false northing? */
402  /* -------------------------------------------------------------------- */
403 
404  {
405  const char *y_0_1 = NULL, *y_0_2 = NULL;
406 
407  y_0_1 = G_find_key_value("y_0", proj_info1);
408  y_0_2 = G_find_key_value("y_0", proj_info2);
409 
410  if ((y_0_1 && !y_0_2) || (!y_0_1 && y_0_2))
411  return -8;
412 
413  if (y_0_1 && y_0_2 && (fabs(atof(y_0_1) - atof(y_0_2)) > 0.000001))
414  return -8;
415  }
416 
417  /* -------------------------------------------------------------------- */
418  /* Do they have the same center longitude? */
419  /* -------------------------------------------------------------------- */
420 
421  {
422  const char *l_1 = NULL, *l_2 = NULL;
423 
424  l_1 = G_find_key_value("lon_0", proj_info1);
425  l_2 = G_find_key_value("lon_0", proj_info2);
426 
427  if ((l_1 && !l_2) || (!l_1 && l_2))
428  return -9;
429 
430  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
431  return -9;
432 
433  /* --------------------------------------------------------------------
434  */
435  /* Do they have the same center latitude? */
436  /* --------------------------------------------------------------------
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  */
450  /* Do they have the same standard parallels? */
451  /* --------------------------------------------------------------------
452  */
453 
454  l_1 = G_find_key_value("lat_1", proj_info1);
455  l_2 = G_find_key_value("lat_1", proj_info2);
456 
457  if ((l_1 && !l_2) || (!l_1 && l_2))
458  return -11;
459 
460  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
461  /* lat_1 differ */
462  /* check for swapped lat_1, lat_2 */
463  l_2 = G_find_key_value("lat_2", proj_info2);
464 
465  if (!l_2)
466  return -11;
467  if (fabs(atof(l_1) - atof(l_2)) > 0.000001) {
468  return -11;
469  }
470  }
471 
472  l_1 = G_find_key_value("lat_2", proj_info1);
473  l_2 = G_find_key_value("lat_2", proj_info2);
474 
475  if ((l_1 && !l_2) || (!l_1 && l_2))
476  return -11;
477 
478  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
479  /* lat_2 differ */
480  /* check for swapped lat_1, lat_2 */
481  l_2 = G_find_key_value("lat_1", proj_info2);
482 
483  if (!l_2)
484  return -11;
485  if (fabs(atof(l_1) - atof(l_2)) > 0.000001) {
486  return -11;
487  }
488  }
489  }
490 
491  /* towgs84 ? */
492 
493  /* -------------------------------------------------------------------- */
494  /* Add more details in later. */
495  /* -------------------------------------------------------------------- */
496 
497  return 1;
498 }
499 
500 /*!
501  \brief Write WKT definition to file
502 
503  Any WKT string and version recognized by PROJ is supported.
504 
505  \param location_name name of the location to write the WKT definition
506  \param wktstring pointer to WKT string
507 
508  \return 0 success
509  \return -1 error writing
510  */
511 
512 int G_write_projwkt(const char *location_name, const char *wktstring)
513 {
514  FILE *fp;
515  char path[GPATH_MAX];
516  int err, n;
517 
518  if (!wktstring)
519  return 0;
520 
521  if (location_name && *location_name)
522  snprintf(path, sizeof(path), "%s/%s/%s/%s", G_gisdbase(), location_name,
523  "PERMANENT", WKT_FILE);
524  else
525  G_file_name(path, "", WKT_FILE, "PERMANENT");
526 
527  fp = fopen(path, "w");
528 
529  if (!fp)
530  G_fatal_error(_("Unable to open output file <%s>: %s"), path,
531  strerror(errno));
532 
533  err = 0;
534  n = strlen(wktstring);
535  if (wktstring[n - 1] != '\n') {
536  if (n != fprintf(fp, "%s\n", wktstring))
537  err = -1;
538  }
539  else {
540  if (n != fprintf(fp, "%s", wktstring))
541  err = -1;
542  }
543 
544  if (fclose(fp) != 0)
545  G_fatal_error(_("Error closing output file <%s>: %s"), path,
546  strerror(errno));
547 
548  return err;
549 }
550 
551 /*!
552  \brief Write srid (spatial reference id) to file
553 
554  A srid consists of an authority name and code and must be known to
555  PROJ.
556 
557  \param location_name name of the location to write the srid
558  \param sridstring pointer to srid string
559 
560  \return 0 success
561  \return -1 error writing
562  */
563 
564 int G_write_projsrid(const char *location_name, const char *sridstring)
565 {
566  FILE *fp;
567  char path[GPATH_MAX];
568  int err, n;
569 
570  if (!sridstring)
571  return 0;
572 
573  if (location_name && *location_name)
574  snprintf(path, sizeof(path), "%s/%s/%s/%s", G_gisdbase(), location_name,
575  "PERMANENT", SRID_FILE);
576  else
577  G_file_name(path, "", SRID_FILE, "PERMANENT");
578 
579  fp = fopen(path, "w");
580 
581  if (!fp)
582  G_fatal_error(_("Unable to open output file <%s>: %s"), path,
583  strerror(errno));
584 
585  err = 0;
586  n = strlen(sridstring);
587  if (sridstring[n - 1] != '\n') {
588  if (n != fprintf(fp, "%s\n", sridstring))
589  err = -1;
590  }
591  else {
592  if (n != fprintf(fp, "%s", sridstring))
593  err = -1;
594  }
595 
596  if (fclose(fp) != 0)
597  G_fatal_error(_("Error closing output file <%s>: %s"), path,
598  strerror(errno));
599 
600  return err;
601 }
#define NULL
Definition: ccmath.h:32
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
char * G_file_name(char *, const char *, const char *, const char *)
Builds full path names to GIS data files.
Definition: file_name.c:61
const char * G_gisdbase(void)
Get name of top level database directory.
Definition: gisdbase.c:26
int G_legal_filename(const char *)
Check for legal database file name.
Definition: legal_name.c:34
const char * G_find_key_value(const char *, const struct Key_Value *)
Find given key (case sensitive)
Definition: key_value1.c:85
int G_get_ellipsoid_by_name(const char *, double *, double *)
Get ellipsoid parameters by name.
Definition: get_ellipse.c:104
void G_write_key_value_file(const char *, const struct Key_Value *)
Write key/value pairs to file.
Definition: key_value3.c:28
int int G_strcasecmp(const char *, const char *)
String compare ignoring case (upper or lower)
Definition: strings.c:47
int G_put_element_window(const struct Cell_head *, const char *, const char *)
Write the region.
Definition: put_window.c:74
int G_mkdir(const char *)
Creates a new directory.
Definition: paths.c:27
int G_debug(int, const char *,...) __attribute__((format(printf
void G_setenv_nogisrc(const char *, const char *)
Set environment name to value (doesn't update .gisrc)
Definition: env.c:472
#define WKT_FILE
Definition: gis.h:136
#define GPATH_MAX
Definition: gis.h:199
#define SRID_FILE
Definition: gis.h:137
#define TRUE
Definition: gis.h:78
#define _(str)
Definition: glocale.h:10
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_write_projsrid(const char *location_name, const char *sridstring)
Write srid (spatial reference id) to file.
Definition: make_loc.c:564
int G_write_projwkt(const char *location_name, const char *wktstring)
Write WKT definition to file.
Definition: make_loc.c:512
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:232
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_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
2D/3D raster map header (used also for region)
Definition: gis.h:446
Definition: gis.h:534
Definition: path.h:15
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:216