Loading...
Searching...
No Matches
make_loc.c
Go to the documentation of this file.
/*!
* \file lib/gis/make_loc.c
*
* \brief GIS Library - Functions to create a new location
*
* Creates a new location automatically given a "Cell_head", PROJ_INFO
* and PROJ_UNITS information.
*
* (C) 2000-2013 by the GRASS Development Team
*
* This program is free software under the GNU General Public License
* (>=v2). Read the file COPYING that comes with GRASS for details.
*
* \author Frank Warmerdam
*/
#include <grass/gis.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <unistd.h>
#include <sys/stat.h>
#include <math.h>
#include <grass/glocale.h>
/*!
* \brief Create a new location
*
* This function creates a new location in the current database,
* initializes the projection, default window and current window.
*
* \param location_name Name of the new location. Should not include
* the full path, the location will be created within
* the current database.
* \param wind default window setting for the new location.
* All fields should be set in this
* structure, and care should be taken to ensure that
* the proj/zone fields match the definition in the
* proj_info parameter(see G_set_cellhd_from_projinfo()).
*
* \param proj_info projection definition suitable to write to the
* PROJ_INFO file, or NULL for PROJECTION_XY.
*
* \param proj_units projection units suitable to write to the PROJ_UNITS
* file, or NULL.
*
* \return 0 on success
* \return -1 to indicate a system error (check errno).
* \return -2 failed to create projection file (currently not used)
* \return -3 illegal name
*/
int G_make_location(const char *location_name, struct Cell_head *wind,
const struct Key_Value *proj_info,
const struct Key_Value *proj_units)
{
char path[GPATH_MAX];
/* check if location name is legal */
if (G_legal_filename(location_name) != 1)
return -3;
/* Try to create the location directory, under the gisdbase. */
snprintf(path, sizeof(path), "%s/%s", G_gisdbase(), location_name);
if (G_mkdir(path) != 0)
return -1;
/* Make the PERMANENT mapset. */
snprintf(path, sizeof(path), "%s/%s/%s", G_gisdbase(), location_name,
"PERMANENT");
if (G_mkdir(path) != 0) {
return -1;
}
/* make these the new current location and mapset */
G_setenv_nogisrc("LOCATION_NAME", location_name);
G_setenv_nogisrc("MAPSET", "PERMANENT");
/* Create the default, and current window files */
G_put_element_window(wind, "", "DEFAULT_WIND");
G_put_element_window(wind, "", "WIND");
/* Write out the PROJ_INFO, and PROJ_UNITS if available. */
if (proj_info != NULL) {
G_file_name(path, "", "PROJ_INFO", "PERMANENT");
G_write_key_value_file(path, proj_info);
}
if (proj_units != NULL) {
G_file_name(path, "", "PROJ_UNITS", "PERMANENT");
G_write_key_value_file(path, proj_units);
}
return 0;
}
/*!
* \brief Create a new location
*
* This function creates a new location in the current database,
* initializes the projection, default window and current window,
* and sets the EPSG code if present
*
* \param location_name Name of the new location. Should not include
* the full path, the location will be created within
* the current database.
* \param wind default window setting for the new location.
* All fields should be set in this
* structure, and care should be taken to ensure that
* the proj/zone fields match the definition in the
* proj_info parameter(see G_set_cellhd_from_projinfo()).
*
* \param proj_info projection definition suitable to write to the
* PROJ_INFO file, or NULL for PROJECTION_XY.
*
* \param proj_units projection units suitable to write to the PROJ_UNITS
* file, or NULL.
*
* \param proj_epsg EPSG code suitable to write to the PROJ_EPSG
* file, or NULL.
*
* \return 0 on success
* \return -1 to indicate a system error (check errno).
* \return -2 failed to create projection file (currently not used)
* \return -3 illegal name
*/
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)
{
int ret;
ret = G_make_location(location_name, wind, proj_info, proj_units);
if (ret != 0)
return ret;
/* Write out the PROJ_EPSG if available. */
if (proj_epsg != NULL) {
char path[GPATH_MAX];
G_file_name(path, "", "PROJ_EPSG", "PERMANENT");
G_write_key_value_file(path, proj_epsg);
}
return 0;
}
/*!
* \brief Create a new location
*
* This function creates a new location in the current database,
* initializes the projection, default window and current window,
* and sets WKT, srid, and EPSG code if present
*
* \param location_name Name of the new location. Should not include
* the full path, the location will be created within
* the current database.
* \param wind default window setting for the new location.
* All fields should be set in this
* structure, and care should be taken to ensure that
* the proj/zone fields match the definition in the
* proj_info parameter(see G_set_cellhd_from_projinfo()).
*
* \param proj_info projection definition suitable to write to the
* PROJ_INFO file, or NULL for PROJECTION_XY.
*
* \param proj_units projection units suitable to write to the PROJ_UNITS
* file, or NULL.
*
* \param proj_srid Spatial reference ID suitable to write to the PROJ_SRID
* file, or NULL.
*
* \param proj_wkt WKT definition suitable to write to the PROJ_WKT
* file, or NULL.
*
* \return 0 on success
* \return -1 to indicate a system error (check errno).
* \return -2 failed to create projection file (currently not used)
* \return -3 illegal name
*/
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)
{
int ret;
ret = G_make_location(location_name, wind, proj_info, proj_units);
if (ret != 0)
return ret;
/* Write out PROJ_SRID if srid is available. */
if (proj_srid != NULL) {
G_write_projsrid(location_name, proj_srid);
}
/* Write out PROJ_WKT if WKT is available. */
if (proj_wkt != NULL) {
G_write_projwkt(location_name, proj_wkt);
}
return 0;
}
/*!
* \brief Compare projections including units
*
* \param proj_info1 projection info to compare
* \param proj_units1 projection units to compare
* \param proj_info2 projection info to compare
* \param proj_units2 projection units to compare
* \return -1 if not the same projection
* \return -2 if linear unit translation to meters fails
* \return -3 if not the same datum,
* \return -4 if not the same ellipsoid,
* \return -5 if UTM zone differs
* \return -6 if UTM hemisphere differs,
* \return -7 if false easting differs
* \return -8 if false northing differs,
* \return -9 if center longitude differs,
* \return -10 if center latitude differs,
* \return -11 if standard parallels differ,
* \return 1 if projections match.
*/
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)
{
const char *proj1, *proj2;
if (proj_info1 == NULL && proj_info2 == NULL)
return TRUE;
/* -------------------------------------------------------------------- */
/* Are they both in the same projection? */
/* -------------------------------------------------------------------- */
/* prevent seg fault in G_find_key_value */
if (proj_info1 == NULL || proj_info2 == NULL)
return -1;
proj1 = G_find_key_value("proj", proj_info1);
proj2 = G_find_key_value("proj", proj_info2);
if (proj1 == NULL || proj2 == NULL || strcmp(proj1, proj2))
return -1;
/* -------------------------------------------------------------------- */
/* Verify that the linear unit translation to meters is OK. */
/* -------------------------------------------------------------------- */
/* prevent seg fault in G_find_key_value */
if (proj_units1 == NULL && proj_units2 == NULL)
return 1;
if (proj_units1 == NULL || proj_units2 == NULL)
return -2;
{
double a1 = 0, a2 = 0;
if (G_find_key_value("meters", proj_units1) != NULL)
a1 = atof(G_find_key_value("meters", proj_units1));
if (G_find_key_value("meters", proj_units2) != NULL)
a2 = atof(G_find_key_value("meters", proj_units2));
if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
return -2;
}
/* compare unit name only if there is no to meter conversion factor */
if (G_find_key_value("meters", proj_units1) == NULL ||
G_find_key_value("meters", proj_units2) == NULL) {
const char *u_1 = NULL, *u_2 = NULL;
u_1 = G_find_key_value("unit", proj_units1);
u_2 = G_find_key_value("unit", proj_units2);
if ((u_1 && !u_2) || (!u_1 && u_2))
return -2;
/* the unit name can be arbitrary: the following can be the same
* us-ft (proj.4 keyword)
* U.S. Surveyor's Foot (proj.4 name)
* US survey foot (WKT)
* Foot_US (WKT)
*/
if (u_1 && u_2 && G_strcasecmp(u_1, u_2))
return -2;
}
/* -------------------------------------------------------------------- */
/* Do they both have the same datum? */
/* -------------------------------------------------------------------- */
{
const char *d_1 = NULL, *d_2 = NULL;
d_1 = G_find_key_value("datum", proj_info1);
d_2 = G_find_key_value("datum", proj_info2);
if ((d_1 && !d_2) || (!d_1 && d_2))
return -3;
if (d_1 && d_2 && strcmp(d_1, d_2)) {
/* different datum short names can mean the same datum,
* see lib/gis/datum.table */
G_debug(1, "Different datum names");
}
}
/* -------------------------------------------------------------------- */
/* Do they both have the same ellipsoid? */
/* -------------------------------------------------------------------- */
{
const char *e_1 = NULL, *e_2 = NULL;
e_1 = G_find_key_value("ellps", proj_info1);
e_2 = G_find_key_value("ellps", proj_info2);
if (e_1 && e_2 && strcmp(e_1, e_2))
return -4;
if (e_1 == NULL || e_2 == NULL) {
double a1 = 0, a2 = 0;
double es1 = 0, es2 = 0;
/* it may happen that one proj_info has ellps,
* while the other has a, es: translate ellps to a, es */
if (e_1)
G_get_ellipsoid_by_name(e_1, &a1, &es1);
else {
if (G_find_key_value("a", proj_info1) != NULL)
a1 = atof(G_find_key_value("a", proj_info1));
if (G_find_key_value("es", proj_info1) != NULL)
es1 = atof(G_find_key_value("es", proj_info1));
}
if (e_2)
G_get_ellipsoid_by_name(e_2, &a2, &es2);
else {
if (G_find_key_value("a", proj_info2) != NULL)
a2 = atof(G_find_key_value("a", proj_info2));
if (G_find_key_value("es", proj_info2) != NULL)
es2 = atof(G_find_key_value("es", proj_info2));
}
/* it should be an error if a = 0 */
if ((a1 == 0 && a2 != 0) || (a1 != 0 && a2 == 0))
return -4;
if (a1 && a2 && (fabs(a2 - a1) > 0.000001))
return -4;
if ((es1 == 0 && es2 != 0) || (es1 != 0 && es2 == 0))
return -4;
if (es1 && es2 && (fabs(es2 - es1) > 0.000001))
return -4;
}
}
/* -------------------------------------------------------------------- */
/* Zone check specially for UTM */
/* -------------------------------------------------------------------- */
if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm") &&
atof(G_find_key_value("zone", proj_info1)) !=
atof(G_find_key_value("zone", proj_info2)))
return -5;
/* -------------------------------------------------------------------- */
/* Hemisphere check specially for UTM */
/* -------------------------------------------------------------------- */
if (!strcmp(proj1, "utm") && !strcmp(proj2, "utm") &&
!!G_find_key_value("south", proj_info1) !=
!!G_find_key_value("south", proj_info2))
return -6;
/* -------------------------------------------------------------------- */
/* Do they both have the same false easting? */
/* -------------------------------------------------------------------- */
{
const char *x_0_1 = NULL, *x_0_2 = NULL;
x_0_1 = G_find_key_value("x_0", proj_info1);
x_0_2 = G_find_key_value("x_0", proj_info2);
if ((x_0_1 && !x_0_2) || (!x_0_1 && x_0_2))
return -7;
if (x_0_1 && x_0_2 && (fabs(atof(x_0_1) - atof(x_0_2)) > 0.000001))
return -7;
}
/* -------------------------------------------------------------------- */
/* Do they both have the same false northing? */
/* -------------------------------------------------------------------- */
{
const char *y_0_1 = NULL, *y_0_2 = NULL;
y_0_1 = G_find_key_value("y_0", proj_info1);
y_0_2 = G_find_key_value("y_0", proj_info2);
if ((y_0_1 && !y_0_2) || (!y_0_1 && y_0_2))
return -8;
if (y_0_1 && y_0_2 && (fabs(atof(y_0_1) - atof(y_0_2)) > 0.000001))
return -8;
}
/* -------------------------------------------------------------------- */
/* Do they have the same center longitude? */
/* -------------------------------------------------------------------- */
{
const char *l_1 = NULL, *l_2 = NULL;
l_1 = G_find_key_value("lon_0", proj_info1);
l_2 = G_find_key_value("lon_0", proj_info2);
if ((l_1 && !l_2) || (!l_1 && l_2))
return -9;
if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
return -9;
/* --------------------------------------------------------------------
*/
/* Do they have the same center latitude? */
/* --------------------------------------------------------------------
*/
l_1 = G_find_key_value("lat_0", proj_info1);
l_2 = G_find_key_value("lat_0", proj_info2);
if ((l_1 && !l_2) || (!l_1 && l_2))
return -10;
if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001))
return -10;
/* --------------------------------------------------------------------
*/
/* Do they have the same standard parallels? */
/* --------------------------------------------------------------------
*/
l_1 = G_find_key_value("lat_1", proj_info1);
l_2 = G_find_key_value("lat_1", proj_info2);
if ((l_1 && !l_2) || (!l_1 && l_2))
return -11;
if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
/* lat_1 differ */
/* check for swapped lat_1, lat_2 */
l_2 = G_find_key_value("lat_2", proj_info2);
if (!l_2)
return -11;
if (fabs(atof(l_1) - atof(l_2)) > 0.000001) {
return -11;
}
}
l_1 = G_find_key_value("lat_2", proj_info1);
l_2 = G_find_key_value("lat_2", proj_info2);
if ((l_1 && !l_2) || (!l_1 && l_2))
return -11;
if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
/* lat_2 differ */
/* check for swapped lat_1, lat_2 */
l_2 = G_find_key_value("lat_1", proj_info2);
if (!l_2)
return -11;
if (fabs(atof(l_1) - atof(l_2)) > 0.000001) {
return -11;
}
}
}
/* towgs84 ? */
/* -------------------------------------------------------------------- */
/* Add more details in later. */
/* -------------------------------------------------------------------- */
return 1;
}
/*!
\brief Write WKT definition to file
Any WKT string and version recognized by PROJ is supported.
\param location_name name of the location to write the WKT definition
\param wktstring pointer to WKT string
\return 0 success
\return -1 error writing
*/
int G_write_projwkt(const char *location_name, const char *wktstring)
{
FILE *fp;
char path[GPATH_MAX];
int err, n;
if (!wktstring)
return 0;
if (location_name && *location_name)
snprintf(path, sizeof(path), "%s/%s/%s/%s", G_gisdbase(), location_name,
"PERMANENT", WKT_FILE);
else
G_file_name(path, "", WKT_FILE, "PERMANENT");
fp = fopen(path, "w");
if (!fp)
G_fatal_error(_("Unable to open output file <%s>: %s"), path,
strerror(errno));
err = 0;
n = strlen(wktstring);
if (wktstring[n - 1] != '\n') {
if (n != fprintf(fp, "%s\n", wktstring))
err = -1;
}
else {
if (n != fprintf(fp, "%s", wktstring))
err = -1;
}
if (fclose(fp) != 0)
G_fatal_error(_("Error closing output file <%s>: %s"), path,
strerror(errno));
return err;
}
/*!
\brief Write srid (spatial reference id) to file
A srid consists of an authority name and code and must be known to
PROJ.
\param location_name name of the location to write the srid
\param sridstring pointer to srid string
\return 0 success
\return -1 error writing
*/
int G_write_projsrid(const char *location_name, const char *sridstring)
{
FILE *fp;
char path[GPATH_MAX];
int err, n;
if (!sridstring)
return 0;
if (location_name && *location_name)
snprintf(path, sizeof(path), "%s/%s/%s/%s", G_gisdbase(), location_name,
"PERMANENT", SRID_FILE);
else
G_file_name(path, "", SRID_FILE, "PERMANENT");
fp = fopen(path, "w");
if (!fp)
G_fatal_error(_("Unable to open output file <%s>: %s"), path,
strerror(errno));
err = 0;
n = strlen(sridstring);
if (sridstring[n - 1] != '\n') {
if (n != fprintf(fp, "%s\n", sridstring))
err = -1;
}
else {
if (n != fprintf(fp, "%s", sridstring))
err = -1;
}
if (fclose(fp) != 0)
G_fatal_error(_("Error closing output file <%s>: %s"), path,
strerror(errno));
return err;
}
NULL
#define NULL
Definition ccmath.h:32
TRUE
#define TRUE
Definition dbfopen.c:56
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
G_setenv_nogisrc
void G_setenv_nogisrc(const char *name, const char *value)
Set environment name to value (doesn't update .gisrc).
Definition env.c:472
G_file_name
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
G_get_ellipsoid_by_name
int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
Get ellipsoid parameters by name.
Definition get_ellipse.c:104
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_gisdbase
const char * G_gisdbase(void)
Get name of top level database directory.
Definition gisdbase.c:26
G_find_key_value
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive).
Definition key_value1.c:85
G_write_key_value_file
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
G_legal_filename
int G_legal_filename(const char *s)
Check for legal database file name.
Definition legal_name.c:34
G_make_location_crs
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:183
G_write_projsrid
int G_write_projsrid(const char *location_name, const char *sridstring)
Write srid (spatial reference id) to file.
Definition make_loc.c:559
G_write_projwkt
int G_write_projwkt(const char *location_name, const char *wktstring)
Write WKT definition to file.
Definition make_loc.c:508
G_compare_projections
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:229
G_make_location
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
G_make_location_epsg
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
G_mkdir
int G_mkdir(const char *path)
Creates a new directory.
Definition paths.c:27
G_put_element_window
int G_put_element_window(const struct Cell_head *window, const char *dir, const char *name)
Write the region.
Definition put_window.c:74
G_strcasecmp
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower).
Definition strings.c:47
path
Definition path.h:15
err
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition symbol/read.c:216