Loading...
Searching...
No Matches
get_proj.c
Go to the documentation of this file.
/**
\file get_proj.c
\brief GProj library - Functions for re-projecting point data
\author Original Author unknown, probably Soil Conservation Service,
Eric Miller, Paul Kelly, Markus Metz
(C) 2003-2008, 2018 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.
**/
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>
#include <string.h>
#include <grass/gis.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
/* Finder function for datum transformation grids */
#define FINDERFUNC set_proj_share
#define PERMANENT "PERMANENT"
#define MAX_PARGS 100
static void alloc_options(char *);
static char *opt_in[MAX_PARGS];
static int nopt;
/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
/**
* \brief Create a pj_info struct Co-ordinate System definition from a set of
* PROJ_INFO / PROJ_UNITS-style key-value pairs
*
* This function takes a GRASS-style co-ordinate system definition as stored
* in the PROJ_INFO and PROJ_UNITS files and processes it to create a pj_info
* representation for use in re-projecting with pj_do_proj(). In addition to
* the parameters passed to it it may also make reference to the system
* ellipse.table and datum.table files if necessary.
*
* \param info Pointer to a pj_info struct (which must already exist) into
* which the co-ordinate system definition will be placed
* \param in_proj_keys PROJ_INFO-style key-value pairs
* \param in_units_keys PROJ_UNITS-style key-value pairs
*
* \return -1 on error (unable to initialise PROJ.4)
* 2 if "default" 3-parameter datum shift values from datum.table
* were used
* 3 if an unrecognised datum name was passed on to PROJ.4 (and
* initialization was successful)
* 1 otherwise
**/
int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
const struct Key_Value *in_units_keys)
{
const char *str;
int i;
double a, es, rf;
int returnval = 1;
char buffa[300], factbuff[50];
int deflen;
char proj_in[250], *datum, *params;
PJ *pj;
PJ_CONTEXT *pjc;
proj_in[0] = '\0';
info->zone = 0;
info->meters = 1.0;
info->proj[0] = '\0';
info->def = NULL;
info->pj = NULL;
info->srid = NULL;
info->wkt = NULL;
str = G_find_key_value("meters", in_units_keys);
if (str != NULL) {
strcpy(factbuff, str);
if (strlen(factbuff) > 0)
sscanf(factbuff, "%lf", &(info->meters));
}
str = G_find_key_value("name", in_proj_keys);
if (str != NULL) {
snprintf(proj_in, sizeof(proj_in), "%s", str);
}
str = G_find_key_value("proj", in_proj_keys);
if (str != NULL) {
snprintf(info->proj, sizeof(info->proj), "%s", str);
}
if (strlen(info->proj) <= 0)
snprintf(info->proj, sizeof(info->proj), "ll");
str = G_find_key_value("init", in_proj_keys);
if (str != NULL) {
info->srid = G_store(str);
}
nopt = 0;
for (i = 0; i < in_proj_keys->nitems; i++) {
/* the name parameter is just for grasses use */
if (strcmp(in_proj_keys->key[i], "name") == 0) {
continue;
/* init is here ignored */
}
else if (strcmp(in_proj_keys->key[i], "init") == 0) {
continue;
/* zone handled separately at end of loop */
}
else if (strcmp(in_proj_keys->key[i], "zone") == 0) {
continue;
/* Datum and ellipsoid-related parameters will be handled
* separately after end of this loop PK */
}
else if (strcmp(in_proj_keys->key[i], "datum") == 0 ||
strcmp(in_proj_keys->key[i], "dx") == 0 ||
strcmp(in_proj_keys->key[i], "dy") == 0 ||
strcmp(in_proj_keys->key[i], "dz") == 0 ||
strcmp(in_proj_keys->key[i], "datumparams") == 0 ||
strcmp(in_proj_keys->key[i], "nadgrids") == 0 ||
strcmp(in_proj_keys->key[i], "towgs84") == 0 ||
strcmp(in_proj_keys->key[i], "ellps") == 0 ||
strcmp(in_proj_keys->key[i], "a") == 0 ||
strcmp(in_proj_keys->key[i], "b") == 0 ||
strcmp(in_proj_keys->key[i], "es") == 0 ||
strcmp(in_proj_keys->key[i], "f") == 0 ||
strcmp(in_proj_keys->key[i], "rf") == 0) {
continue;
/* PROJ.4 uses longlat instead of ll as 'projection name' */
}
else if (strcmp(in_proj_keys->key[i], "proj") == 0) {
if (strcmp(in_proj_keys->value[i], "ll") == 0)
snprintf(buffa, sizeof(buffa), "proj=longlat");
else
snprintf(buffa, sizeof(buffa), "proj=%s",
in_proj_keys->value[i]);
/* 'One-sided' PROJ.4 flags will have the value in
* the key-value pair set to 'defined' and only the
* key needs to be passed on. */
}
else if (strcmp(in_proj_keys->value[i], "defined") == 0)
snprintf(buffa, sizeof(buffa), "%s", in_proj_keys->key[i]);
else
snprintf(buffa, sizeof(buffa), "%s=%s", in_proj_keys->key[i],
in_proj_keys->value[i]);
alloc_options(buffa);
}
str = G_find_key_value("zone", in_proj_keys);
if (str != NULL) {
if (sscanf(str, "%d", &(info->zone)) != 1) {
G_fatal_error(_("Invalid zone %s specified"), str);
}
if (info->zone < 0) {
/* if zone is negative, write abs(zone) and define south */
info->zone = -info->zone;
if (G_find_key_value("south", in_proj_keys) == NULL) {
snprintf(buffa, sizeof(buffa), "south");
alloc_options(buffa);
}
}
snprintf(buffa, sizeof(buffa), "zone=%d", info->zone);
alloc_options(buffa);
}
if ((GPJ__get_ellipsoid_params(in_proj_keys, &a, &es, &rf) == 0) &&
(str = G_find_key_value("ellps", in_proj_keys)) != NULL) {
/* Default values were returned but an ellipsoid name not recognised
* by GRASS is present---perhaps it will be recognised by
* PROJ.4 even though it wasn't by GRASS */
snprintf(buffa, sizeof(buffa), "ellps=%s", str);
alloc_options(buffa);
}
else {
snprintf(buffa, sizeof(buffa), "a=%.16g", a);
alloc_options(buffa);
/* Cannot use es directly because the OSRImportFromProj4()
* function in OGR only accepts b or rf as the 2nd parameter */
if (es == 0)
snprintf(buffa, sizeof(buffa), "b=%.16g", a);
else
snprintf(buffa, sizeof(buffa), "rf=%.16g", rf);
alloc_options(buffa);
}
/* Workaround to stop PROJ reading values from defaults file when
* rf (and sometimes ellps) is not specified */
if (G_find_key_value("no_defs", in_proj_keys) == NULL) {
snprintf(buffa, sizeof(buffa), "no_defs");
alloc_options(buffa);
}
/* If datum parameters are present in the PROJ_INFO keys, pass them on */
if (GPJ__get_datum_params(in_proj_keys, &datum, &params) == 2) {
snprintf(buffa, sizeof(buffa), "%s", params);
alloc_options(buffa);
G_free(params);
/* else if a datum name is present take it and look up the parameters
* from the datum.table file */
}
else if (datum != NULL) {
if (GPJ_get_default_datum_params_by_name(datum, &params) > 0) {
snprintf(buffa, sizeof(buffa), "%s", params);
alloc_options(buffa);
returnval = 2;
G_free(params);
/* else just pass the datum name on and hope it is recognised by
* PROJ.4 even though it isn't recognised by GRASS */
}
else {
snprintf(buffa, sizeof(buffa), "datum=%s", datum);
alloc_options(buffa);
returnval = 3;
}
/* else there'll be no datum transformation taking place here... */
}
else {
returnval = 4;
}
G_free(datum);
/* without type=crs, PROJ6 does not recognize what this is,
* a crs or some kind of coordinate operation, falling through to
* PJ_TYPE_OTHER_COORDINATE_OPERATION */
alloc_options("type=crs");
pjc = proj_context_create();
if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
strcpy(
buffa,
_("Unable to initialise PROJ with the following parameter list:"));
for (i = 0; i < nopt; i++) {
char err[50];
snprintf(err, sizeof(err), " +%s", opt_in[i]);
strcat(buffa, err);
}
G_warning("%s", buffa);
return -1;
}
int perr = proj_errno(pj);
if (perr)
G_fatal_error("PROJ error %d", perr);
if (proj_get_type(pj) == PJ_TYPE_BOUND_CRS) {
PJ *source_crs = proj_get_source_crs(pjc, pj);
if (source_crs) {
proj_destroy(pj);
pj = source_crs;
}
}
info->pj = pj;
deflen = 0;
for (i = 0; i < nopt; i++)
deflen += strlen(opt_in[i]) + 2;
info->def = G_malloc(deflen + 1);
snprintf(buffa, sizeof(buffa), "+%s ", opt_in[0]);
strcpy(info->def, buffa);
G_free(opt_in[0]);
for (i = 1; i < nopt; i++) {
snprintf(buffa, sizeof(buffa), "+%s ", opt_in[i]);
strcat(info->def, buffa);
G_free(opt_in[i]);
}
return returnval;
}
static void alloc_options(char *buffa)
{
size_t nsize;
nsize = strlen(buffa) + 1;
opt_in[nopt++] = (char *)G_malloc(nsize);
snprintf(opt_in[nopt - 1], nsize, "%s", buffa);
return;
}
/**
* \brief Create a pj_info struct Co-ordinate System definition from a
* string with a sequence of key=value pairs
*
* This function takes a GRASS- or PROJ style co-ordinate system definition
* and processes it to create a pj_info representation for use in
* re-projecting with pj_do_proj(). In addition to the parameters passed
* to it it may also make reference to the system ellipse.table and
* datum.table files if necessary.
*
* \param info Pointer to a pj_info struct (which must already exist) into
* which the co-ordinate system definition will be placed
* \param str input string with projection definition
*
* \return -1 on error (unable to initialise PROJ.4)
* 1 on success
**/
int pj_get_string(struct pj_info *info, char *str)
{
char *s;
int i, nsize;
char zonebuff[50], buffa[300];
int deflen;
PJ *pj;
PJ_CONTEXT *pjc;
info->zone = 0;
info->proj[0] = '\0';
info->meters = 1.0;
info->def = NULL;
info->srid = NULL;
info->pj = NULL;
nopt = 0;
if ((str == NULL) || (str[0] == '\0')) {
/* Null Pointer or empty string is supplied for parameters,
* implying latlong projection; just need to set proj
* parameter and call pj_init PK */
snprintf(info->proj, sizeof(info->proj), "ll");
snprintf(buffa, sizeof(buffa), "proj=latlong ellps=WGS84");
alloc_options(buffa);
}
else {
/* Parameters have been provided; parse through them but don't
* bother with most of the checks in pj_get_kv; assume the
* programmer knows what he / she is doing when using this
* function rather than reading a PROJ_INFO file PK */
s = str;
while (s = strtok(s, " \t\n"), s) {
if (strncmp(s, "+unfact=", 8) == 0) {
s = s + 8;
info->meters = atof(s);
}
else {
if (strncmp(s, "+", 1) == 0)
++s;
if (nsize = strlen(s), nsize) {
if (nopt >= MAX_PARGS) {
fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
G_fatal_error(
_("Option input overflowed option table"));
}
if (strncmp("zone=", s, 5) == 0) {
snprintf(zonebuff, sizeof(zonebuff), "%s", s + 5);
sscanf(zonebuff, "%d", &(info->zone));
}
if (strncmp(s, "init=", 5) == 0) {
info->srid = G_store(s + 6);
}
if (strncmp("proj=", s, 5) == 0) {
snprintf(info->proj, sizeof(info->proj), "%s", s + 5);
if (strcmp(info->proj, "ll") == 0)
snprintf(buffa, sizeof(buffa), "proj=latlong");
else
snprintf(buffa, sizeof(buffa), "%s", s);
}
else {
snprintf(buffa, sizeof(buffa), "%s", s);
}
alloc_options(buffa);
}
}
s = 0;
}
}
/* without type=crs, PROJ6 does not recognize what this is,
* a crs or some kind of coordinate operation, falling through to
* PJ_TYPE_OTHER_COORDINATE_OPERATION */
alloc_options("type=crs");
pjc = proj_context_create();
if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
G_warning(_("Unable to initialize pj cause: %s"),
proj_errno_string(proj_context_errno(pjc)));
return -1;
}
if (proj_get_type(pj) == PJ_TYPE_BOUND_CRS) {
PJ *source_crs = proj_get_source_crs(pjc, pj);
if (source_crs) {
proj_destroy(pj);
pj = source_crs;
}
}
info->pj = pj;
deflen = 0;
for (i = 0; i < nopt; i++)
deflen += strlen(opt_in[i]) + 2;
info->def = G_malloc(deflen + 1);
snprintf(buffa, sizeof(buffa), "+%s ", opt_in[0]);
strcpy(info->def, buffa);
G_free(opt_in[0]);
for (i = 1; i < nopt; i++) {
snprintf(buffa, sizeof(buffa), "+%s ", opt_in[i]);
strcat(info->def, buffa);
G_free(opt_in[i]);
}
return 1;
}
/* set_proj_share()
* 'finder function' for use with PROJ.4 pj_set_finder() function
* this is used to find grids, usually in /usr/share/proj
* GRASS no longer provides copies of proj grids in GRIDDIR
* -> do not use gisbase/GRIDDIR */
const char *set_proj_share(const char *name)
{
static char *buf = NULL;
const char *projshare;
static size_t buf_len = 0;
size_t len;
projshare = getenv("GRASS_PROJSHARE");
if (!projshare)
return NULL;
len = strlen(projshare) + strlen(name) + 2;
if (buf_len < len) {
if (buf != NULL)
G_free(buf);
buf_len = len + 20;
buf = G_malloc(buf_len);
}
snprintf(buf, buf_len, "%s/%s", projshare, name);
return buf;
}
/**
* \brief Print projection parameters as used by PROJ.4 for input and
* output co-ordinate systems
*
* \param iproj 'Input' co-ordinate system
* \param oproj 'Output' co-ordinate system
*
* \return 1 on success, -1 on error (i.e. if the PROJ-style definition
* is NULL for either co-ordinate system)
**/
int pj_print_proj_params(const struct pj_info *iproj,
const struct pj_info *oproj)
{
char *str;
if (iproj) {
str = iproj->def;
if (str != NULL) {
fprintf(stderr, "%s: %s\n", _("Input Projection Parameters"), str);
fprintf(stderr, "%s: %.16g\n", _("Input Unit Factor"),
iproj->meters);
}
else
return -1;
}
if (oproj) {
str = oproj->def;
if (str != NULL) {
fprintf(stderr, "%s: %s\n", _("Output Projection Parameters"), str);
fprintf(stderr, "%s: %.16g\n", _("Output Unit Factor"),
oproj->meters);
}
else
return -1;
}
return 1;
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
NULL
#define NULL
Definition ccmath.h:32
GPJ__get_ellipsoid_params
int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a, double *e2, double *rf)
Get the ellipsoid parameters from proj keys structure.
Definition ellipse.c:74
pj_get_kv
int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys, const struct Key_Value *in_units_keys)
Create a pj_info struct Co-ordinate System definition from a set of PROJ_INFO / PROJ_UNITS-style key-...
Definition get_proj.c:60
MAX_PARGS
#define MAX_PARGS
Definition get_proj.c:28
pj_print_proj_params
int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *oproj)
Print projection parameters as used by PROJ.4 for input and output co-ordinate systems.
Definition get_proj.c:474
set_proj_share
const char * set_proj_share(const char *name)
Definition get_proj.c:438
pj_get_string
int pj_get_string(struct pj_info *info, char *str)
Create a pj_info struct Co-ordinate System definition from a string with a sequence of key=value pair...
Definition get_proj.c:319
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
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
name
const char * name
Definition named_colr.c:6
strcpy
#define strcpy
Definition parson.c:66
GPJ_get_default_datum_params_by_name
int GPJ_get_default_datum_params_by_name(const char *name, char **params)
"Last resort" function to retrieve a "default" set of datum parameters for a datum (N....
Definition proj/datum.c:86
GPJ__get_datum_params
int GPJ__get_datum_params(const struct Key_Value *projinfo, char **datumname, char **params)
Extract the datum transformation-related parameters from a set of general PROJ_INFO parameters.
Definition proj/datum.c:173
G_store
char * G_store(const char *s)
Copy string to allocated memory.
Definition strings.c:87
err
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition symbol/read.c:216