GRASS GIS 7 Programmer's Manual  7.4.1(2018)-exported
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 Compare projections including units
97  *
98  * \param proj_info1 projection info to compare
99  * \param proj_units1 projection units to compare
100  * \param proj_info2 projection info to compare
101  * \param proj_units2 projection units to compare
102 
103  * \return -1 if not the same projection
104  * \return -2 if linear unit translation to meters fails
105  * \return -3 if not the same datum,
106  * \return -4 if not the same ellipsoid,
107  * \return -5 if UTM zone differs
108  * \return -6 if UTM hemisphere differs,
109  * \return -7 if false easting differs
110  * \return -8 if false northing differs,
111  * \return -9 if center longitude differs,
112  * \return -10 if center latitude differs,
113  * \return -11 if standard parallels differ,
114  * \return 1 if projections match.
115  */
116 int G_compare_projections(const struct Key_Value *proj_info1,
117  const struct Key_Value *proj_units1,
118  const struct Key_Value *proj_info2,
119  const struct Key_Value *proj_units2)
120 {
121  const char *proj1, *proj2;
122 
123  if (proj_info1 == NULL && proj_info2 == NULL)
124  return TRUE;
125 
126  /* -------------------------------------------------------------------- */
127  /* Are they both in the same projection? */
128  /* -------------------------------------------------------------------- */
129  /* prevent seg fault in G_find_key_value */
130  if (proj_info1 == NULL || proj_info2 == NULL)
131  return -1;
132 
133  proj1 = G_find_key_value("proj", proj_info1);
134  proj2 = G_find_key_value("proj", proj_info2);
135 
136  if (proj1 == NULL || proj2 == NULL || strcmp(proj1, proj2))
137  return -1;
138 
139  /* -------------------------------------------------------------------- */
140  /* Verify that the linear unit translation to meters is OK. */
141  /* -------------------------------------------------------------------- */
142  /* prevent seg fault in G_find_key_value */
143  if (proj_units1 == NULL && proj_units2 == NULL)
144  return 1;
145 
146  if (proj_units1 == NULL || proj_units2 == NULL)
147  return -2;
148 
149  {
150  double a1 = 0, a2 = 0;
151 
152  if (G_find_key_value("meters", proj_units1) != NULL)
153  a1 = atof(G_find_key_value("meters", proj_units1));
154  if (G_find_key_value("meters", proj_units2) != NULL)
155  a2 = atof(G_find_key_value("meters", proj_units2));
156 
157  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
158  return -2;
159  }
160  /* compare unit name only if there is no to meter conversion factor */
161  if (G_find_key_value("meters", proj_units1) == NULL ||
162  G_find_key_value("meters", proj_units2) == NULL) {
163  const char *u_1 = NULL, *u_2 = NULL;
164 
165  u_1 = G_find_key_value("unit", proj_units1);
166  u_2 = G_find_key_value("unit", proj_units2);
167 
168  if ((u_1 && !u_2) || (!u_1 && u_2))
169  return -2;
170 
171  /* the unit name can be arbitrary: the following can be the same
172  * us-ft (proj.4 keyword)
173  * U.S. Surveyor's Foot (proj.4 name)
174  * US survey foot (WKT)
175  * Foot_US (WKT)
176  */
177  if (u_1 && u_2 && G_strcasecmp(u_1, u_2))
178  return -2;
179  }
180 
181  /* -------------------------------------------------------------------- */
182  /* Do they both have the same datum? */
183  /* -------------------------------------------------------------------- */
184  {
185  const char *d_1 = NULL, *d_2 = NULL;
186 
187  d_1 = G_find_key_value("datum", proj_info1);
188  d_2 = G_find_key_value("datum", proj_info2);
189 
190  if ((d_1 && !d_2) || (!d_1 && d_2))
191  return -3;
192 
193  if (d_1 && d_2 && strcmp(d_1, d_2)) {
194  /* different datum short names can mean the same datum,
195  * see lib/gis/datum.table */
196  G_debug(1, "Different datum names");
197  }
198  }
199 
200  /* -------------------------------------------------------------------- */
201  /* Do they both have the same ellipsoid? */
202  /* -------------------------------------------------------------------- */
203  {
204  const char *e_1 = NULL, *e_2 = NULL;
205 
206  e_1 = G_find_key_value("ellps", proj_info1);
207  e_2 = G_find_key_value("ellps", proj_info2);
208 
209  if (e_1 && e_2 && strcmp(e_1, e_2))
210  return -4;
211 
212  if (e_1 == NULL || e_2 == NULL) {
213  double a1 = 0, a2 = 0;
214  double es1 = 0, es2 = 0;
215 
216  /* it may happen that one proj_info has ellps,
217  * while the other has a, es: translate ellps to a, es */
218  if (e_1)
219  G_get_ellipsoid_by_name(e_1, &a1, &es1);
220  else {
221  if (G_find_key_value("a", proj_info1) != NULL)
222  a1 = atof(G_find_key_value("a", proj_info1));
223  if (G_find_key_value("es", proj_info1) != NULL)
224  es1 = atof(G_find_key_value("es", proj_info1));
225  }
226 
227  if (e_2)
228  G_get_ellipsoid_by_name(e_2, &a2, &es2);
229  else {
230  if (G_find_key_value("a", proj_info2) != NULL)
231  a2 = atof(G_find_key_value("a", proj_info2));
232  if (G_find_key_value("es", proj_info2) != NULL)
233  es2 = atof(G_find_key_value("es", proj_info2));
234  }
235 
236  /* it should be an error if a = 0 */
237  if ((a1 == 0 && a2 != 0) || (a1 != 0 && a2 == 0))
238  return -4;
239 
240  if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
241  return -4;
242 
243  if ((es1 == 0 && es2 != 0) || (es1 != 0 && es2 == 0))
244  return -4;
245 
246  if (es1 && es2 && (fabs(es2 - es1) > 0.000001))
247  return -4;
248  }
249  }
250 
251  /* -------------------------------------------------------------------- */
252  /* Zone check specially for UTM */
253  /* -------------------------------------------------------------------- */
254  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
255  && atof(G_find_key_value("zone", proj_info1))
256  != atof(G_find_key_value("zone", proj_info2)))
257  return -5;
258 
259  /* -------------------------------------------------------------------- */
260  /* Hemisphere check specially for UTM */
261  /* -------------------------------------------------------------------- */
262  if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm")
263  && !!G_find_key_value("south", proj_info1)
264  != !!G_find_key_value("south", proj_info2))
265  return -6;
266 
267  /* -------------------------------------------------------------------- */
268  /* Do they both have the same false easting? */
269  /* -------------------------------------------------------------------- */
270 
271  {
272  const char *x_0_1 = NULL, *x_0_2 = NULL;
273 
274  x_0_1 = G_find_key_value("x_0", proj_info1);
275  x_0_2 = G_find_key_value("x_0", proj_info2);
276 
277  if ((x_0_1 && !x_0_2) || (!x_0_1 && x_0_2))
278  return -7;
279 
280  if (x_0_1 && x_0_2 && (fabs(atof(x_0_1) - atof(x_0_2)) > 0.000001))
281  return -7;
282  }
283 
284  /* -------------------------------------------------------------------- */
285  /* Do they both have the same false northing? */
286  /* -------------------------------------------------------------------- */
287 
288  {
289  const char *y_0_1 = NULL, *y_0_2 = NULL;
290 
291  y_0_1 = G_find_key_value("y_0", proj_info1);
292  y_0_2 = G_find_key_value("y_0", proj_info2);
293 
294  if ((y_0_1 && !y_0_2) || (!y_0_1 && y_0_2))
295  return -8;
296 
297  if (y_0_1 && y_0_2 && (fabs(atof(y_0_1) - atof(y_0_2)) > 0.000001))
298  return -8;
299  }
300 
301  /* -------------------------------------------------------------------- */
302  /* Do they have the same center longitude? */
303  /* -------------------------------------------------------------------- */
304 
305  {
306  const char *l_1 = NULL, *l_2 = NULL;
307 
308  l_1 = G_find_key_value("lon_0", proj_info1);
309  l_2 = G_find_key_value("lon_0", proj_info2);
310 
311  if ((l_1 && !l_2) || (!l_1 && l_2))
312  return -9;
313 
314  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
315  return -9;
316 
317  /* -------------------------------------------------------------------- */
318  /* Do they have the same center latitude? */
319  /* -------------------------------------------------------------------- */
320 
321  l_1 = l_2 = NULL;
322  l_1 = G_find_key_value("lat_0", proj_info1);
323  l_2 = G_find_key_value("lat_0", proj_info2);
324 
325  if ((l_1 && !l_2) || (!l_1 && l_2))
326  return -10;
327 
328  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
329  return -10;
330 
331  /* -------------------------------------------------------------------- */
332  /* Do they have the same standard parallels? */
333  /* -------------------------------------------------------------------- */
334 
335  l_1 = l_2 = NULL;
336  l_1 = G_find_key_value("lat_1", proj_info1);
337  l_2 = G_find_key_value("lat_1", proj_info2);
338 
339  if ((l_1 && !l_2) || (!l_1 && l_2))
340  return -11;
341 
342  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
343  /* lat_1 differ */
344  /* check for swapped lat_1, lat_2 */
345  l_2 = G_find_key_value("lat_2", proj_info2);
346 
347  if (!l_2)
348  return -11;
349  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
350  return -11;
351  }
352  }
353 
354  l_1 = l_2 = NULL;
355  l_1 = G_find_key_value("lat_2", proj_info1);
356  l_2 = G_find_key_value("lat_2", proj_info2);
357 
358  if ((l_1 && !l_2) || (!l_1 && l_2))
359  return -11;
360 
361  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
362  /* lat_2 differ */
363  /* check for swapped lat_1, lat_2 */
364  l_2 = G_find_key_value("lat_1", proj_info2);
365 
366  if (!l_2)
367  return -11;
368  if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
369  return -11;
370  }
371  }
372  }
373 
374  /* towgs84 ? */
375 
376  /* -------------------------------------------------------------------- */
377  /* Add more details in later. */
378  /* -------------------------------------------------------------------- */
379 
380  return 1;
381 }
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
#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:116
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
#define TRUE
Definition: dbfopen.c:183
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
int G_put_element_window(const struct Cell_head *window, const char *dir, const char *name)
Write the region.
Definition: put_window.c:74
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