GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72636
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 <unistd.h>
22 #include <sys/stat.h>
23 #include <math.h>
24 
25 /*!
26  * \brief Create a new location
27  *
28  * This function creates a new location in the current database,
29  * initializes the projection, default window and current window.
30  *
31  * \param location_name Name of the new location. Should not include
32  * the full path, the location will be created within
33  * the current database.
34  * \param wind default window setting for the new location.
35  * All fields should be set in this
36  * structure, and care should be taken to ensure that
37  * the proj/zone fields match the definition in the
38  * proj_info parameter(see G_set_cellhd_from_projinfo()).
39  *
40  * \param proj_info projection definition suitable to write to the
41  * PROJ_INFO file, or NULL for PROJECTION_XY.
42  *
43  * \param proj_units projection units suitable to write to the PROJ_UNITS
44  * file, or NULL.
45  *
46  * \return 0 on success
47  * \return -1 to indicate a system error (check errno).
48  * \return -2 failed to create projection file (currently not used)
49  * \return -3 illegal name
50  */
51 int G_make_location(const char *location_name,
52  struct Cell_head *wind,
53  const struct Key_Value *proj_info,
54  const struct Key_Value *proj_units)
55 {
56  char path[GPATH_MAX];
57 
58  /* check if location name is legal */
59  if (G_legal_filename(location_name) != 1)
60  return -3;
61 
62  /* Try to create the location directory, under the gisdbase. */
63  sprintf(path, "%s/%s", G_gisdbase(), location_name);
64  if (G_mkdir(path) != 0)
65  return -1;
66 
67  /* Make the PERMANENT mapset. */
68  sprintf(path, "%s/%s/%s", G_gisdbase(), location_name, "PERMANENT");
69  if (G_mkdir(path) != 0) {
70  return -1;
71  }
72 
73  /* make these the new current location and mapset */
74  G_setenv_nogisrc("LOCATION_NAME", location_name);
75  G_setenv_nogisrc("MAPSET", "PERMANENT");
76 
77  /* Create the default, and current window files */
78  G_put_element_window(wind, "", "DEFAULT_WIND");
79  G_put_element_window(wind, "", "WIND");
80 
81  /* Write out the PROJ_INFO, and PROJ_UNITS if available. */
82  if (proj_info != NULL) {
83  G_file_name(path, "", "PROJ_INFO", "PERMANENT");
84  G_write_key_value_file(path, proj_info);
85  }
86 
87  if (proj_units != NULL) {
88  G_file_name(path, "", "PROJ_UNITS", "PERMANENT");
89  G_write_key_value_file(path, proj_units);
90  }
91 
92  return 0;
93 }
94 
95 /*!
96  * \brief Create a new location
97  *
98  * This function creates a new location in the current database,
99  * initializes the projection, default window and current window,
100  * and sets the EPSG code if present
101  *
102  * \param location_name Name of the new location. Should not include
103  * the full path, the location will be created within
104  * the current database.
105  * \param wind default window setting for the new location.
106  * All fields should be set in this
107  * structure, and care should be taken to ensure that
108  * the proj/zone fields match the definition in the
109  * proj_info parameter(see G_set_cellhd_from_projinfo()).
110  *
111  * \param proj_info projection definition suitable to write to the
112  * PROJ_INFO file, or NULL for PROJECTION_XY.
113  *
114  * \param proj_units projection units suitable to write to the PROJ_UNITS
115  * file, or NULL.
116  *
117  * \param proj_epsg EPSG code suitable to write to the PROJ_EPSG
118  * file, or NULL.
119  *
120  * \return 0 on success
121  * \return -1 to indicate a system error (check errno).
122  * \return -2 failed to create projection file (currently not used)
123  * \return -3 illegal name
124  */
125 int G_make_location_epsg(const char *location_name,
126  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  char path[GPATH_MAX];
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  G_file_name(path, "", "PROJ_EPSG", "PERMANENT");
142  G_write_key_value_file(path, proj_epsg);
143  }
144 
145  return 0;
146 }
147 
148 /*!
149  * \brief Compare projections including units
150  *
151  * \param proj_info1 projection info to compare
152  * \param proj_units1 projection units to compare
153  * \param proj_info2 projection info to compare
154  * \param proj_units2 projection units to compare
155 
156  * \return -1 if not the same projection
157  * \return -2 if linear unit translation to meters fails
158  * \return -3 if not the same datum,
159  * \return -4 if not the same ellipsoid,
160  * \return -5 if UTM zone differs
161  * \return -6 if UTM hemisphere differs,
162  * \return -7 if false easting differs
163  * \return -8 if false northing differs,
164  * \return -9 if center longitude differs,
165  * \return -10 if center latitude differs,
166  * \return -11 if standard parallels differ,
167  * \return 1 if projections match.
168  */
169 int G_compare_projections(const struct Key_Value *proj_info1,
170  const struct Key_Value *proj_units1,
171  const struct Key_Value *proj_info2,
172  const struct Key_Value *proj_units2)
173 {
174  const char *proj1, *proj2;
175 
176  if (proj_info1 == NULL && proj_info2 == NULL)
177  return TRUE;
178 
179  /* -------------------------------------------------------------------- */
180  /* Are they both in the same projection? */
181  /* -------------------------------------------------------------------- */
182  /* prevent seg fault in G_find_key_value */
183  if (proj_info1 == NULL || proj_info2 == NULL)
184  return -1;
185 
186  proj1 = G_find_key_value("proj", proj_info1);
187  proj2 = G_find_key_value("proj", proj_info2);
188 
189  if (proj1 == NULL || proj2 == NULL || strcmp(proj1, proj2))
190  return -1;
191 
192  /* -------------------------------------------------------------------- */
193  /* Verify that the linear unit translation to meters is OK. */
194  /* -------------------------------------------------------------------- */
195  /* prevent seg fault in G_find_key_value */
196  if (proj_units1 == NULL && proj_units2 == NULL)
197  return 1;
198 
199  if (proj_units1 == NULL || proj_units2 == NULL)
200  return -2;
201 
202  {
203  double a1 = 0, a2 = 0;
204 
205  if (G_find_key_value("meters", proj_units1) != NULL)
206  a1 = atof(G_find_key_value("meters", proj_units1));
207  if (G_find_key_value("meters", proj_units2) != NULL)
208  a2 = atof(G_find_key_value("meters", proj_units2));
209 
210  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
211  return -2;
212  }
213  /* compare unit name only if there is no to meter conversion factor */
214  if (G_find_key_value("meters", proj_units1) == NULL ||
215  G_find_key_value("meters", proj_units2) == NULL) {
216  const char *u_1 = NULL, *u_2 = NULL;
217 
218  u_1 = G_find_key_value("unit", proj_units1);
219  u_2 = G_find_key_value("unit", proj_units2);
220 
221  if ((u_1 && !u_2) || (!u_1 && u_2))
222  return -2;
223 
224  /* the unit name can be arbitrary: the following can be the same
225  * us-ft (proj.4 keyword)
226  * U.S. Surveyor's Foot (proj.4 name)
227  * US survey foot (WKT)
228  * Foot_US (WKT)
229  */
230  if (u_1 && u_2 && G_strcasecmp(u_1, u_2))
231  return -2;
232  }
233 
234  /* -------------------------------------------------------------------- */
235  /* Do they both have the same datum? */
236  /* -------------------------------------------------------------------- */
237  {
238  const char *d_1 = NULL, *d_2 = NULL;
239 
240  d_1 = G_find_key_value("datum", proj_info1);
241  d_2 = G_find_key_value("datum", proj_info2);
242 
243  if ((d_1 && !d_2) || (!d_1 && d_2))
244  return -3;
245 
246  if (d_1 && d_2 && strcmp(d_1, d_2)) {
247  /* different datum short names can mean the same datum,
248  * see lib/gis/datum.table */
249  G_debug(1, "Different datum names");
250  }
251  }
252 
253  /* -------------------------------------------------------------------- */
254  /* Do they both have the same ellipsoid? */
255  /* -------------------------------------------------------------------- */
256  {
257  const char *e_1 = NULL, *e_2 = NULL;
258 
259  e_1 = G_find_key_value("ellps", proj_info1);
260  e_2 = G_find_key_value("ellps", proj_info2);
261 
262  if (e_1 && e_2 && strcmp(e_1, e_2))
263  return -4;
264 
265  if (e_1 == NULL || e_2 == NULL) {
266  double a1 = 0, a2 = 0;
267  double es1 = 0, es2 = 0;
268 
269  /* it may happen that one proj_info has ellps,
270  * while the other has a, es: translate ellps to a, es */
271  if (e_1)
272  G_get_ellipsoid_by_name(e_1, &a1, &es1);
273  else {
274  if (G_find_key_value("a", proj_info1) != NULL)
275  a1 = atof(G_find_key_value("a", proj_info1));
276  if (G_find_key_value("es", proj_info1) != NULL)
277  es1 = atof(G_find_key_value("es", proj_info1));
278  }
279 
280  if (e_2)
281  G_get_ellipsoid_by_name(e_2, &a2, &es2);
282  else {
283  if (G_find_key_value("a", proj_info2) != NULL)
284  a2 = atof(G_find_key_value("a", proj_info2));
285  if (G_find_key_value("es", proj_info2) != NULL)
286  es2 = atof(G_find_key_value("es", proj_info2));
287  }
288 
289  /* it should be an error if a = 0 */
290  if ((a1 == 0 && a2 != 0) || (a1 != 0 && a2 == 0))
291  return -4;
292 
293  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
294  return -4;
295 
296  if ((es1 == 0 && es2 != 0) || (es1 != 0 && es2 == 0))
297  return -4;
298 
299  if (es1 && es2 && (fabs(es2 - es1) > 0.000001))
300  return -4;
301  }
302  }
303 
304  /* -------------------------------------------------------------------- */
305  /* Zone check specially for UTM */
306  /* -------------------------------------------------------------------- */
307  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
308  && atof(G_find_key_value("zone", proj_info1))
309  != atof(G_find_key_value("zone", proj_info2)))
310  return -5;
311 
312  /* -------------------------------------------------------------------- */
313  /* Hemisphere check specially for UTM */
314  /* -------------------------------------------------------------------- */
315  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
316  && !!G_find_key_value("south", proj_info1)
317  != !!G_find_key_value("south", proj_info2))
318  return -6;
319 
320  /* -------------------------------------------------------------------- */
321  /* Do they both have the same false easting? */
322  /* -------------------------------------------------------------------- */
323 
324  {
325  const char *x_0_1 = NULL, *x_0_2 = NULL;
326 
327  x_0_1 = G_find_key_value("x_0", proj_info1);
328  x_0_2 = G_find_key_value("x_0", proj_info2);
329 
330  if ((x_0_1 && !x_0_2) || (!x_0_1 && x_0_2))
331  return -7;
332 
333  if (x_0_1 && x_0_2 && (fabs(atof(x_0_1) - atof(x_0_2)) > 0.000001))
334  return -7;
335  }
336 
337  /* -------------------------------------------------------------------- */
338  /* Do they both have the same false northing? */
339  /* -------------------------------------------------------------------- */
340 
341  {
342  const char *y_0_1 = NULL, *y_0_2 = NULL;
343 
344  y_0_1 = G_find_key_value("y_0", proj_info1);
345  y_0_2 = G_find_key_value("y_0", proj_info2);
346 
347  if ((y_0_1 && !y_0_2) || (!y_0_1 && y_0_2))
348  return -8;
349 
350  if (y_0_1 && y_0_2 && (fabs(atof(y_0_1) - atof(y_0_2)) > 0.000001))
351  return -8;
352  }
353 
354  /* -------------------------------------------------------------------- */
355  /* Do they have the same center longitude? */
356  /* -------------------------------------------------------------------- */
357 
358  {
359  const char *l_1 = NULL, *l_2 = NULL;
360 
361  l_1 = G_find_key_value("lon_0", proj_info1);
362  l_2 = G_find_key_value("lon_0", proj_info2);
363 
364  if ((l_1 && !l_2) || (!l_1 && l_2))
365  return -9;
366 
367  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
368  return -9;
369 
370  /* -------------------------------------------------------------------- */
371  /* Do they have the same center latitude? */
372  /* -------------------------------------------------------------------- */
373 
374  l_1 = l_2 = NULL;
375  l_1 = G_find_key_value("lat_0", proj_info1);
376  l_2 = G_find_key_value("lat_0", proj_info2);
377 
378  if ((l_1 && !l_2) || (!l_1 && l_2))
379  return -10;
380 
381  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
382  return -10;
383 
384  /* -------------------------------------------------------------------- */
385  /* Do they have the same standard parallels? */
386  /* -------------------------------------------------------------------- */
387 
388  l_1 = l_2 = NULL;
389  l_1 = G_find_key_value("lat_1", proj_info1);
390  l_2 = G_find_key_value("lat_1", proj_info2);
391 
392  if ((l_1 && !l_2) || (!l_1 && l_2))
393  return -11;
394 
395  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
396  /* lat_1 differ */
397  /* check for swapped lat_1, lat_2 */
398  l_2 = G_find_key_value("lat_2", proj_info2);
399 
400  if (!l_2)
401  return -11;
402  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
403  return -11;
404  }
405  }
406 
407  l_1 = l_2 = NULL;
408  l_1 = G_find_key_value("lat_2", proj_info1);
409  l_2 = G_find_key_value("lat_2", proj_info2);
410 
411  if ((l_1 && !l_2) || (!l_1 && l_2))
412  return -11;
413 
414  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
415  /* lat_2 differ */
416  /* check for swapped lat_1, lat_2 */
417  l_2 = G_find_key_value("lat_1", proj_info2);
418 
419  if (!l_2)
420  return -11;
421  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
422  return -11;
423  }
424  }
425  }
426 
427  /* towgs84 ? */
428 
429  /* -------------------------------------------------------------------- */
430  /* Add more details in later. */
431  /* -------------------------------------------------------------------- */
432 
433  return 1;
434 }
#define TRUE
Definition: gis.h:49
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
Definition: strings.c:46
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
int G_mkdir(const char *path)
Creates a new directory.
Definition: paths.c:27
2D/3D raster map header (used also for region)
Definition: gis.h:390
#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:169
void G_write_key_value_file(const char *file, const struct Key_Value *kv)
Write key/value pairs to file.
Definition: key_value3.c:28
char * G_file_name(char *path, const char *element, const char *name, const char *mapset)
Builds full path names to GIS data files.
Definition: file_name.c:38
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
void G_setenv_nogisrc(const char *name, const char *value)
Set environment name to value (doesn&#39;t update .gisrc)
Definition: env.c:448
#define GPATH_MAX
Definition: gis.h:151
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:125
int G_put_element_window(const struct Cell_head *window, const char *dir, const char *name)
Write the region.
Definition: put_window.c:74
Definition: gis.h:479
int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
Get ellipsoid parameters by name.
Definition: get_ellipse.c:105
const char * G_gisdbase(void)
Get name of top level database directory.
Definition: gisdbase.c:26
Definition: path.h:16
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:51