GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-8fe13636b
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  char path[GPATH_MAX];
135 
136  ret = G_make_location(location_name, wind, proj_info, proj_units);
137 
138  if (ret != 0)
139  return ret;
140 
141  /* Write out the PROJ_EPSG if available. */
142  if (proj_epsg != NULL) {
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  char path[GPATH_MAX];
195 
196  ret = G_make_location(location_name, wind, proj_info, proj_units);
197 
198  if (ret != 0)
199  return ret;
200 
201  /* Write out PROJ_SRID if srid is available. */
202  if (proj_srid != NULL) {
203  G_write_projsrid(location_name, proj_srid);
204  }
205 
206  /* Write out PROJ_WKT if WKT is available. */
207  if (proj_wkt != NULL) {
208  G_write_projwkt(location_name, proj_wkt);
209  }
210 
211  return 0;
212 }
213 
214 /*!
215  * \brief Compare projections including units
216  *
217  * \param proj_info1 projection info to compare
218  * \param proj_units1 projection units to compare
219  * \param proj_info2 projection info to compare
220  * \param proj_units2 projection units to compare
221 
222  * \return -1 if not the same projection
223  * \return -2 if linear unit translation to meters fails
224  * \return -3 if not the same datum,
225  * \return -4 if not the same ellipsoid,
226  * \return -5 if UTM zone differs
227  * \return -6 if UTM hemisphere differs,
228  * \return -7 if false easting differs
229  * \return -8 if false northing differs,
230  * \return -9 if center longitude differs,
231  * \return -10 if center latitude differs,
232  * \return -11 if standard parallels differ,
233  * \return 1 if projections match.
234  */
235 int G_compare_projections(const struct Key_Value *proj_info1,
236  const struct Key_Value *proj_units1,
237  const struct Key_Value *proj_info2,
238  const struct Key_Value *proj_units2)
239 {
240  const char *proj1, *proj2;
241 
242  if (proj_info1 == NULL && proj_info2 == NULL)
243  return TRUE;
244 
245  /* -------------------------------------------------------------------- */
246  /* Are they both in the same projection? */
247  /* -------------------------------------------------------------------- */
248  /* prevent seg fault in G_find_key_value */
249  if (proj_info1 == NULL || proj_info2 == NULL)
250  return -1;
251 
252  proj1 = G_find_key_value("proj", proj_info1);
253  proj2 = G_find_key_value("proj", proj_info2);
254 
255  if (proj1 == NULL || proj2 == NULL || strcmp(proj1, proj2))
256  return -1;
257 
258  /* -------------------------------------------------------------------- */
259  /* Verify that the linear unit translation to meters is OK. */
260  /* -------------------------------------------------------------------- */
261  /* prevent seg fault in G_find_key_value */
262  if (proj_units1 == NULL && proj_units2 == NULL)
263  return 1;
264 
265  if (proj_units1 == NULL || proj_units2 == NULL)
266  return -2;
267 
268  {
269  double a1 = 0, a2 = 0;
270 
271  if (G_find_key_value("meters", proj_units1) != NULL)
272  a1 = atof(G_find_key_value("meters", proj_units1));
273  if (G_find_key_value("meters", proj_units2) != NULL)
274  a2 = atof(G_find_key_value("meters", proj_units2));
275 
276  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
277  return -2;
278  }
279  /* compare unit name only if there is no to meter conversion factor */
280  if (G_find_key_value("meters", proj_units1) == NULL ||
281  G_find_key_value("meters", proj_units2) == NULL) {
282  const char *u_1 = NULL, *u_2 = NULL;
283 
284  u_1 = G_find_key_value("unit", proj_units1);
285  u_2 = G_find_key_value("unit", proj_units2);
286 
287  if ((u_1 && !u_2) || (!u_1 && u_2))
288  return -2;
289 
290  /* the unit name can be arbitrary: the following can be the same
291  * us-ft (proj.4 keyword)
292  * U.S. Surveyor's Foot (proj.4 name)
293  * US survey foot (WKT)
294  * Foot_US (WKT)
295  */
296  if (u_1 && u_2 && G_strcasecmp(u_1, u_2))
297  return -2;
298  }
299 
300  /* -------------------------------------------------------------------- */
301  /* Do they both have the same datum? */
302  /* -------------------------------------------------------------------- */
303  {
304  const char *d_1 = NULL, *d_2 = NULL;
305 
306  d_1 = G_find_key_value("datum", proj_info1);
307  d_2 = G_find_key_value("datum", proj_info2);
308 
309  if ((d_1 && !d_2) || (!d_1 && d_2))
310  return -3;
311 
312  if (d_1 && d_2 && strcmp(d_1, d_2)) {
313  /* different datum short names can mean the same datum,
314  * see lib/gis/datum.table */
315  G_debug(1, "Different datum names");
316  }
317  }
318 
319  /* -------------------------------------------------------------------- */
320  /* Do they both have the same ellipsoid? */
321  /* -------------------------------------------------------------------- */
322  {
323  const char *e_1 = NULL, *e_2 = NULL;
324 
325  e_1 = G_find_key_value("ellps", proj_info1);
326  e_2 = G_find_key_value("ellps", proj_info2);
327 
328  if (e_1 && e_2 && strcmp(e_1, e_2))
329  return -4;
330 
331  if (e_1 == NULL || e_2 == NULL) {
332  double a1 = 0, a2 = 0;
333  double es1 = 0, es2 = 0;
334 
335  /* it may happen that one proj_info has ellps,
336  * while the other has a, es: translate ellps to a, es */
337  if (e_1)
338  G_get_ellipsoid_by_name(e_1, &a1, &es1);
339  else {
340  if (G_find_key_value("a", proj_info1) != NULL)
341  a1 = atof(G_find_key_value("a", proj_info1));
342  if (G_find_key_value("es", proj_info1) != NULL)
343  es1 = atof(G_find_key_value("es", proj_info1));
344  }
345 
346  if (e_2)
347  G_get_ellipsoid_by_name(e_2, &a2, &es2);
348  else {
349  if (G_find_key_value("a", proj_info2) != NULL)
350  a2 = atof(G_find_key_value("a", proj_info2));
351  if (G_find_key_value("es", proj_info2) != NULL)
352  es2 = atof(G_find_key_value("es", proj_info2));
353  }
354 
355  /* it should be an error if a = 0 */
356  if ((a1 == 0 && a2 != 0) || (a1 != 0 && a2 == 0))
357  return -4;
358 
359  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
360  return -4;
361 
362  if ((es1 == 0 && es2 != 0) || (es1 != 0 && es2 == 0))
363  return -4;
364 
365  if (es1 && es2 && (fabs(es2 - es1) > 0.000001))
366  return -4;
367  }
368  }
369 
370  /* -------------------------------------------------------------------- */
371  /* Zone check specially for UTM */
372  /* -------------------------------------------------------------------- */
373  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
374  && atof(G_find_key_value("zone", proj_info1))
375  != atof(G_find_key_value("zone", proj_info2)))
376  return -5;
377 
378  /* -------------------------------------------------------------------- */
379  /* Hemisphere check specially for UTM */
380  /* -------------------------------------------------------------------- */
381  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
382  && !!G_find_key_value("south", proj_info1)
383  != !!G_find_key_value("south", proj_info2))
384  return -6;
385 
386  /* -------------------------------------------------------------------- */
387  /* Do they both have the same false easting? */
388  /* -------------------------------------------------------------------- */
389 
390  {
391  const char *x_0_1 = NULL, *x_0_2 = NULL;
392 
393  x_0_1 = G_find_key_value("x_0", proj_info1);
394  x_0_2 = G_find_key_value("x_0", proj_info2);
395 
396  if ((x_0_1 && !x_0_2) || (!x_0_1 && x_0_2))
397  return -7;
398 
399  if (x_0_1 && x_0_2 && (fabs(atof(x_0_1) - atof(x_0_2)) > 0.000001))
400  return -7;
401  }
402 
403  /* -------------------------------------------------------------------- */
404  /* Do they both have the same false northing? */
405  /* -------------------------------------------------------------------- */
406 
407  {
408  const char *y_0_1 = NULL, *y_0_2 = NULL;
409 
410  y_0_1 = G_find_key_value("y_0", proj_info1);
411  y_0_2 = G_find_key_value("y_0", proj_info2);
412 
413  if ((y_0_1 && !y_0_2) || (!y_0_1 && y_0_2))
414  return -8;
415 
416  if (y_0_1 && y_0_2 && (fabs(atof(y_0_1) - atof(y_0_2)) > 0.000001))
417  return -8;
418  }
419 
420  /* -------------------------------------------------------------------- */
421  /* Do they have the same center longitude? */
422  /* -------------------------------------------------------------------- */
423 
424  {
425  const char *l_1 = NULL, *l_2 = NULL;
426 
427  l_1 = G_find_key_value("lon_0", proj_info1);
428  l_2 = G_find_key_value("lon_0", proj_info2);
429 
430  if ((l_1 && !l_2) || (!l_1 && l_2))
431  return -9;
432 
433  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
434  return -9;
435 
436  /* -------------------------------------------------------------------- */
437  /* Do they have the same center latitude? */
438  /* -------------------------------------------------------------------- */
439 
440  l_1 = l_2 = NULL;
441  l_1 = G_find_key_value("lat_0", proj_info1);
442  l_2 = G_find_key_value("lat_0", proj_info2);
443 
444  if ((l_1 && !l_2) || (!l_1 && l_2))
445  return -10;
446 
447  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
448  return -10;
449 
450  /* -------------------------------------------------------------------- */
451  /* Do they have the same standard parallels? */
452  /* -------------------------------------------------------------------- */
453 
454  l_1 = l_2 = NULL;
455  l_1 = G_find_key_value("lat_1", proj_info1);
456  l_2 = G_find_key_value("lat_1", proj_info2);
457 
458  if ((l_1 && !l_2) || (!l_1 && l_2))
459  return -11;
460 
461  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
462  /* lat_1 differ */
463  /* check for swapped lat_1, lat_2 */
464  l_2 = G_find_key_value("lat_2", proj_info2);
465 
466  if (!l_2)
467  return -11;
468  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
469  return -11;
470  }
471  }
472 
473  l_1 = l_2 = NULL;
474  l_1 = G_find_key_value("lat_2", proj_info1);
475  l_2 = G_find_key_value("lat_2", proj_info2);
476 
477  if ((l_1 && !l_2) || (!l_1 && l_2))
478  return -11;
479 
480  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
481  /* lat_2 differ */
482  /* check for swapped lat_1, lat_2 */
483  l_2 = G_find_key_value("lat_1", proj_info2);
484 
485  if (!l_2)
486  return -11;
487  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
488  return -11;
489  }
490  }
491  }
492 
493  /* towgs84 ? */
494 
495  /* -------------------------------------------------------------------- */
496  /* Add more details in later. */
497  /* -------------------------------------------------------------------- */
498 
499  return 1;
500 }
501 
502 /*!
503  \brief Write WKT definition to file
504 
505  Any WKT string and version recognized by PROJ is supported.
506 
507  \param location_name name of the location to write the WKT definition
508  \param wktstring pointer to WKT string
509 
510  \return 0 success
511  \return -1 error writing
512  */
513 
514 int G_write_projwkt(const char *location_name, const char *wktstring)
515 {
516  FILE *fp;
517  char path[GPATH_MAX];
518  int err, n;
519 
520  if (!wktstring)
521  return 0;
522 
523  if (location_name && *location_name)
524  sprintf(path, "%s/%s/%s/%s", G_gisdbase(), location_name, "PERMANENT", WKT_FILE);
525  else
526  G_file_name(path, "", WKT_FILE, "PERMANENT");
527 
528  fp = fopen(path, "w");
529 
530  if (!fp)
531  G_fatal_error(_("Unable to open output file <%s>: %s"), path, 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, strerror(errno));
546 
547  return err;
548 }
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  sprintf(path, "%s/%s/%s/%s", G_gisdbase(), location_name, "PERMANENT", 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, strerror(errno));
582 
583  err = 0;
584  n = strlen(sridstring);
585  if (sridstring[n - 1] != '\n') {
586  if (n != fprintf(fp, "%s\n", sridstring))
587  err = -1;
588  }
589  else {
590  if (n != fprintf(fp, "%s", sridstring))
591  err = -1;
592  }
593 
594  if (fclose(fp) != 0)
595  G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
596 
597  return err;
598 }
#define TRUE
Definition: gis.h:56
char * G_file_name(char *, const char *, const char *, const char *)
Builds full path names to GIS data files.
Definition: file_name.c:38
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:450
2D/3D raster map header (used also for region)
Definition: gis.h:409
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:235
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:167
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:564
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:498
#define SRID_FILE
Definition: gis.h:121
Definition: path.h:16
int G_mkdir(const char *)
Creates a new directory.
Definition: paths.c:27
#define WKT_FILE
Definition: gis.h:120
#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:514