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