GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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 */
53int 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 */
126int 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 */
185int 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 */
231int 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
511int 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
563int 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
#define TRUE
Definition dbfopen.c:75
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
void G_setenv_nogisrc(const char *name, const char *value)
Set environment name to value (doesn't update .gisrc)
Definition env.c:472
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:61
int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
Get ellipsoid parameters by name.
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
const char * G_gisdbase(void)
Get name of top level database directory.
Definition gisdbase.c:26
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition key_value1.c:85
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
int G_legal_filename(const char *s)
Check for legal database file name.
Definition legal_name.c:34
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
int G_mkdir(const char *path)
Creates a new directory.
Definition paths.c:27
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_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
Definition strings.c:47
Definition path.h:15
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)