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