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 */
G_make_location(const char * location_name,struct Cell_head * wind,const struct Key_Value * proj_info,const struct Key_Value * proj_units)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 */
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)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 defintion 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 */
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)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 */
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)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
G_write_projwkt(const char * location_name,const char * wktstring)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
G_write_projsrid(const char * location_name,const char * sridstring)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 }
599