Loading...
Searching...
No Matches
ellipse.c
Go to the documentation of this file.
/*!
\file lib/proj/ellipse.c
\brief GProj library - Functions for reading datum parameters from the
location database
\author Paul Kelly <paul-grass stjohnspoint.co.uk>
(C) 2003-2008 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 <unistd.h>
#include <ctype.h>
#include <string.h>
#include <stdlib.h>
#include <math.h> /* for sqrt() */
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/gprojects.h>
#include "local_proto.h"
static int get_a_e2_rf(const char *, const char *, double *, double *,
double *);
/*!
* \brief Get the ellipsoid parameters from the database.
*
* If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
* that file, otherwise return WGS 84 values.
*
* Dies with diagnostic if there is an error.
*
* \param[out] a semi-major axis
* \param[out] e2 first eccentricity squared
* \param[out] rf reciprocal of the ellipsoid flattening term
*
* \return 1 on success
* \return 0 default values used.
*/
int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
{
int ret;
struct Key_Value *proj_keys = G_get_projinfo();
if (proj_keys == NULL)
proj_keys = G_create_key_value();
ret = GPJ__get_ellipsoid_params(proj_keys, a, e2, rf);
G_free_key_value(proj_keys);
return ret;
}
/*!
* \brief Get the ellipsoid parameters from proj keys structure.
*
* If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
* that file, otherwise return WGS 84 values.
*
* Dies with diagnostic if there is an error.
*
* \param proj_keys proj definition
* \param[out] a semi-major axis
* \param[out] e2 first eccentricity squared
* \param[out] rf reciprocal of the ellipsoid flattening term
*
* \return 1 on success
* \return 0 default values used.
*/
int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a,
double *e2, double *rf)
{
struct gpj_ellps estruct;
struct gpj_datum dstruct;
const char *str, *str3;
char *str1, *ellps;
str = G_find_key_value("datum", proj_keys);
if ((str != NULL) && (GPJ_get_datum_by_name(str, &dstruct) > 0)) {
/* If 'datum' key is present, look up correct ellipsoid
* from datum.table */
ellps = G_store(dstruct.ellps);
GPJ_free_datum(&dstruct);
}
else
/* else use ellipsoid defined in PROJ_INFO */
ellps = G_store(G_find_key_value("ellps", proj_keys));
if (ellps != NULL && *ellps) {
if (GPJ_get_ellipsoid_by_name(ellps, &estruct) < 0)
G_fatal_error(_("Invalid ellipsoid <%s> in file"), ellps);
*a = estruct.a;
*e2 = estruct.es;
*rf = estruct.rf;
GPJ_free_ellps(&estruct);
G_free(ellps);
return 1;
}
else {
if (ellps) /* *ellps = '\0' */
G_free(ellps);
str3 = G_find_key_value("a", proj_keys);
if (str3 != NULL) {
char *str4;
G_asprintf(&str4, "a=%s", str3);
if ((str3 = G_find_key_value("es", proj_keys)) != NULL)
G_asprintf(&str1, "e=%s", str3);
else if ((str3 = G_find_key_value("f", proj_keys)) != NULL)
G_asprintf(&str1, "f=1/%s", str3);
else if ((str3 = G_find_key_value("rf", proj_keys)) != NULL)
G_asprintf(&str1, "f=1/%s", str3);
else if ((str3 = G_find_key_value("b", proj_keys)) != NULL)
G_asprintf(&str1, "b=%s", str3);
else
G_fatal_error(_("No secondary ellipsoid descriptor "
"(rf, es or b) in file"));
if (get_a_e2_rf(str4, str1, a, e2, rf) == 0)
G_fatal_error(_("Invalid ellipsoid descriptors "
"(a, rf, es or b) in file"));
return 1;
}
else {
str = G_find_key_value("proj", proj_keys);
if ((str == NULL) || (strcmp(str, "ll") == 0)) {
*a = 6378137.0;
*e2 = .006694385;
*rf = 298.257223563;
return 0;
}
else {
G_fatal_error(_("No ellipsoid info given in file"));
}
}
}
return 1;
}
/*!
* \brief Looks up ellipsoid in ellipsoid table and returns the a, e2
* parameters for the ellipsoid.
*
* \param name ellipsoid name
* \param[out] estruct ellipsoid
*
* \return 1 on success
* \return -1 if not found in table
*/
int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
{
struct ellps_list *list, *listhead;
list = listhead = read_ellipsoid_table(0);
while (list != NULL) {
if (G_strcasecmp(name, list->name) == 0) {
estruct->name = G_store(list->name);
estruct->longname = G_store(list->longname);
estruct->a = list->a;
estruct->es = list->es;
estruct->rf = list->rf;
free_ellps_list(listhead);
return 1;
}
list = list->next;
}
free_ellps_list(listhead);
return -1;
}
int get_a_e2_rf(const char *s1, const char *s2, double *a, double *e2,
double *recipf)
{
double b, f;
if (sscanf(s1, "a=%lf", a) != 1)
return 0;
if (*a <= 0.0)
return 0;
if (sscanf(s2, "e=%lf", e2) == 1) {
f = 1.0 - sqrt(1.0 - *e2);
*recipf = 1.0 / f;
return (*e2 >= 0.0);
}
if (sscanf(s2, "f=1/%lf", recipf) == 1) {
if (*recipf <= 0.0)
return 0;
f = 1.0 / *recipf;
*e2 = f * (2 - f);
return (*e2 >= 0.0);
}
if (sscanf(s2, "b=%lf", &b) == 1) {
if (b <= 0.0)
return 0;
if (b == *a) {
f = 0.0;
*e2 = 0.0;
}
else {
f = (*a - b) / *a;
*e2 = f * (2 - f);
}
*recipf = 1.0 / f;
return (*e2 >= 0.0);
}
return 0;
}
struct ellps_list *read_ellipsoid_table(int fatal)
{
FILE *fd;
char file[GPATH_MAX];
char buf[4096];
char name[100], descr[1024], buf1[1024], buf2[1024];
char badlines[1024];
int line;
int err;
struct ellps_list *current = NULL, *outputlist = NULL;
double a, e2, rf;
snprintf(file, sizeof(file), "%s%s", G_gisbase(), ELLIPSOIDTABLE);
fd = fopen(file, "r");
if (!fd) {
(fatal ? G_fatal_error : G_warning)(
_("Unable to open ellipsoid table file <%s>"), file);
return NULL;
}
err = 0;
*badlines = 0;
for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
G_strip(buf);
if (*buf == 0 || *buf == '#')
continue;
if (sscanf(buf, "%s \"%1023[^\"]\" %s %s", name, descr, buf1, buf2) !=
4) {
err++;
snprintf(buf, sizeof(buf), " %d", line);
if (*badlines)
strcat(badlines, ",");
strcat(badlines, buf);
continue;
}
if (get_a_e2_rf(buf1, buf2, &a, &e2, &rf) ||
get_a_e2_rf(buf2, buf1, &a, &e2, &rf)) {
if (current == NULL)
current = outputlist = G_malloc(sizeof(struct ellps_list));
else
current = current->next = G_malloc(sizeof(struct ellps_list));
current->name = G_store(name);
current->longname = G_store(descr);
current->a = a;
current->es = e2;
current->rf = rf;
current->next = NULL;
}
else {
err++;
snprintf(buf, sizeof(buf), " %d", line);
if (*badlines)
strcat(badlines, ",");
strcat(badlines, buf);
continue;
}
}
fclose(fd);
if (!err)
return outputlist;
(fatal ? G_fatal_error : G_warning)(
n_(("Line%s of ellipsoid table file <%s> is invalid"),
("Lines%s of ellipsoid table file <%s> are invalid"), err),
badlines, file);
return outputlist;
}
/*!
\brief Free ellipsoid data structure.
\param estruct data structure to be freed
*/
void GPJ_free_ellps(struct gpj_ellps *estruct)
{
G_free(estruct->name);
G_free(estruct->longname);
return;
}
void free_ellps_list(struct ellps_list *elist)
{
struct ellps_list *old;
while (elist != NULL) {
G_free(elist->name);
G_free(elist->longname);
old = elist;
elist = old->next;
G_free(old);
}
return;
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
G_asprintf
int G_asprintf(char **out, const char *fmt,...)
Definition asprintf.c:68
NULL
#define NULL
Definition ccmath.h:32
b
double b
Definition driver/set_window.c:5
GPJ_get_ellipsoid_by_name
int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
Looks up ellipsoid in ellipsoid table and returns the a, e2 parameters for the ellipsoid.
Definition ellipse.c:159
read_ellipsoid_table
struct ellps_list * read_ellipsoid_table(int fatal)
Definition ellipse.c:223
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
free_ellps_list
void free_ellps_list(struct ellps_list *elist)
Definition ellipse.c:309
GPJ_free_ellps
void GPJ_free_ellps(struct gpj_ellps *estruct)
Free ellipsoid data structure.
Definition ellipse.c:302
GPJ_get_ellipsoid_params
int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
Get the ellipsoid parameters from the database.
Definition ellipse.c:44
G_get_projinfo
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
Definition get_projinfo.c:61
G_getl2
int G_getl2(char *buf, int n, FILE *fd)
Gets a line of text from a file of any pedigree.
Definition getl.c:60
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_gisbase
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition gisbase.c:39
G_free_key_value
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition key_value1.c:104
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_create_key_value
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition key_value1.c:23
file
#define file
name
const char * name
Definition named_colr.c:6
GPJ_free_datum
void GPJ_free_datum(struct gpj_datum *dstruct)
Free the memory used for the strings in a gpj_datum struct.
Definition proj/datum.c:396
GPJ_get_datum_by_name
int GPJ_get_datum_by_name(const char *name, struct gpj_datum *dstruct)
Look up a string in datum.table file to see if it is a valid datum name and if so place its informati...
Definition proj/datum.c:38
list
struct list * list
Definition read_list.c:24
G_strcasecmp
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower).
Definition strings.c:47
G_store
char * G_store(const char *s)
Copy string to allocated memory.
Definition strings.c:87
G_strip
void G_strip(char *buf)
Removes all leading and trailing white space from string.
Definition strings.c:300
err
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition symbol/read.c:216