GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72636
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
get_ellipse.c
Go to the documentation of this file.
1 /*!
2  \file lib/gis/get_ellipse.c
3 
4  \brief GIS Library - Getting ellipsoid parameters from the database.
5 
6  This routine returns the ellipsoid parameters from the database.
7  If the PROJECTION_FILE exists in the PERMANENT mapset, read info
8  from that file, otherwise return WGS 84 values.
9 
10  New 05/2000 by al: for datum shift the f parameter is needed too.
11  This all is not a clean design, but it keeps backward-
12  compatibility.
13  Looks up ellipsoid in ellipsoid table and returns the
14  a, e2 and f parameters for the ellipsoid
15 
16  (C) 2001-2009 by the GRASS Development Team
17 
18  This program is free software under the GNU General Public License
19  (>=v2). Read the file COPYING that comes with GRASS for details.
20 
21  \author CERL
22  */
23 
24 #include <unistd.h>
25 #include <ctype.h>
26 #include <string.h>
27 #include <stdlib.h>
28 #include <math.h> /* for sqrt() */
29 #include <grass/gis.h>
30 #include <grass/glocale.h>
31 
32 static const char PERMANENT[] = "PERMANENT";
33 
34 static struct table {
35  struct ellipse
36  {
37  char *name;
38  char *descr;
39  double a;
40  double e2;
41  double f;
42  } *ellipses;
43  int count;
44  int size;
45  int initialized;
46 } table;
47 
48 /* static int get_a_e2 (char *, char *, double *,double *); */
49 static int get_a_e2_f(const char *, const char *, double *, double *, double *);
50 static int compare_ellipse_names(const void *, const void *);
51 static int get_ellipsoid_parameters(struct Key_Value *, double *, double *);
52 
53 /*!
54  * \brief get ellipsoid parameters
55  *
56  * This routine returns the semi-major axis <b>a</b> (in meters) and
57  * the eccentricity squared <b>e2</b> for the ellipsoid associated
58  * with the database. If there is no ellipsoid explicitly associated
59  * with the database, it returns the values for the WGS 84 ellipsoid.
60  *
61  * \param[out] a semi-major axis
62  * \param[out] e2 eccentricity squared
63  *
64  * \return 1 success
65  * \return 0 default values used
66  */
67 int G_get_ellipsoid_parameters(double *a, double *e2)
68 {
69  int stat;
70  char ipath[GPATH_MAX];
71  struct Key_Value *proj_keys;
72 
73  proj_keys = NULL;
74 
76 
77  if (access(ipath, 0) != 0) {
78  *a = 6378137.0;
79  *e2 = .006694385;
80  return 0;
81  }
82 
83  proj_keys = G_read_key_value_file(ipath);
84 
85  stat = get_ellipsoid_parameters(proj_keys, a, e2);
86 
87  G_free_key_value(proj_keys);
88 
89  return stat;
90 }
91 
92 /*!
93  * \brief Get ellipsoid parameters by name
94  *
95  * This routine returns the semi-major axis <i>a</i> (in meters) and
96  * eccentricity squared <i>e2</i> for the named ellipsoid.
97  *
98  * \param name ellipsoid name
99  * \param[out] a semi-major axis
100  * \param[out] e2 eccentricity squared
101  *
102  * \return 1 on success
103  * \return 0 if ellipsoid not found
104  */
105 int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
106 {
107  int i;
108 
110 
111  for (i = 0; i < table.count; i++) {
112  if (G_strcasecmp(name, table.ellipses[i].name) == 0) {
113  *a = table.ellipses[i].a;
114  *e2 = table.ellipses[i].e2;
115  return 1;
116  }
117  }
118  return 0;
119 }
120 
121 /*!
122  * \brief Get ellipsoid name
123  *
124  * This function returns a pointer to the short name for the
125  * <i>n</i><i>th</i> ellipsoid. If <i>n</i> is less than 0 or greater
126  * than the number of known ellipsoids, it returns a NULL pointer.
127  *
128  * \param n ellipsoid identificator
129  *
130  * \return ellipsoid name
131  * \return NULL if no ellipsoid found
132  */
133 const char *G_ellipsoid_name(int n)
134 {
136  return n >= 0 && n < table.count ? table.ellipses[n].name : NULL;
137 }
138 
139 /*!
140  * \brief Get spheroid parameters by name
141  *
142  * This function returns the semi-major axis <i>a</i> (in meters), the
143  * eccentricity squared <i>e2</i> and the inverse flattening <i>f</i>
144  * for the named ellipsoid.
145  *
146  * \param name spheroid name
147  * \param[out] a semi-major axis
148  * \param[out] e2 eccentricity squared
149  * \param[out] f inverse flattening
150  *
151  * \return 1 on success
152  * \return 0 if no spheroid found
153  */
154 int G_get_spheroid_by_name(const char *name, double *a, double *e2, double *f)
155 {
156  int i;
157 
159 
160  for (i = 0; i < table.count; i++) {
161  if (G_strcasecmp(name, table.ellipses[i].name) == 0) {
162  *a = table.ellipses[i].a;
163  *e2 = table.ellipses[i].e2;
164  *f = table.ellipses[i].f;
165  return 1;
166  }
167  }
168  return 0;
169 }
170 
171 
172 /*!
173  * \brief Get description for nth ellipsoid
174  *
175  * This function returns a pointer to the description text for the
176  * <i>n</i>th ellipsoid. If <i>n</i> is less than 0 or greater
177  * than the number of known ellipsoids, it returns a NULL pointer.
178  *
179  * \param n ellipsoid identificator
180  *
181  * \return pointer to ellipsoid description
182  * \return NULL if no ellipsoid found
183  */
184 const char *G_ellipsoid_description(int n)
185 {
187  return n >= 0 && n < table.count ? table.ellipses[n].descr : NULL;
188 }
189 
190 static int get_a_e2_f(const char *s1, const char *s2, double *a, double *e2, double *f)
191 {
192  double b, recipf;
193 
194  if (sscanf(s1, "a=%lf", a) != 1)
195  return 0;
196 
197  if (*a <= 0.0)
198  return 0;
199 
200  if (sscanf(s2, "e=%lf", e2) == 1) {
201  *f = (double)1.0 / -sqrt(((double)1.0 - *e2)) + (double)1.0;
202  return (*e2 >= 0.0);
203  }
204 
205  if (sscanf(s2, "f=1/%lf", f) == 1) {
206  if (*f <= 0.0)
207  return 0;
208  recipf = (double)1.0 / (*f);
209  *e2 = recipf + recipf - recipf * recipf;
210  return (*e2 >= 0.0);
211  }
212 
213  if (sscanf(s2, "b=%lf", &b) == 1) {
214  if (b <= 0.0)
215  return 0;
216  if (b == *a) {
217  *f = 0.0;
218  *e2 = 0.0;
219  }
220  else {
221  recipf = ((*a) - b) / (*a);
222  *f = (double)1.0 / recipf;
223  *e2 = recipf + recipf - recipf * recipf;
224  }
225  return (*e2 >= 0.0);
226  }
227  return 0;
228 }
229 
230 static int compare_ellipse_names(const void *pa, const void *pb)
231 {
232  const struct ellipse *a = pa;
233  const struct ellipse *b = pb;
234 
235  return G_strcasecmp(a->name, b->name);
236 }
237 
238 /*!
239  \brief Read ellipsoid table
240 
241  \param fatal non-zero value for G_fatal_error(), otherwise
242  G_warning() is used
243 
244  \return 1 on success
245  \return 0 on error
246 */
248 {
249  FILE *fd;
250  char file[GPATH_MAX];
251  char buf[1024];
252  char badlines[256];
253  int line;
254  int err;
255 
256  if (G_is_initialized(&table.initialized))
257  return 1;
258 
259  sprintf(file, "%s/etc/proj/ellipse.table", G_gisbase());
260  fd = fopen(file, "r");
261 
262  if (fd == NULL) {
263  (fatal ? G_fatal_error : G_warning)(_("Unable to open ellipsoid table file <%s>"), file);
264  G_initialize_done(&table.initialized);
265  return 0;
266  }
267 
268  err = 0;
269  *badlines = 0;
270  for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
271  char name[100], descr[100], buf1[100], buf2[100];
272  struct ellipse *e;
273 
274  G_strip(buf);
275  if (*buf == 0 || *buf == '#')
276  continue;
277 
278  if (sscanf(buf, "%s \"%99[^\"]\" %s %s", name, descr, buf1, buf2) != 4) {
279  err++;
280  sprintf(buf, " %d", line);
281  if (*badlines)
282  strcat(badlines, ",");
283  strcat(badlines, buf);
284  continue;
285  }
286 
287  if (table.count >= table.size) {
288  table.size += 60;
289  table.ellipses = G_realloc(table.ellipses, table.size * sizeof(struct ellipse));
290  }
291 
292  e = &table.ellipses[table.count];
293 
294  e->name = G_store(name);
295  e->descr = G_store(descr);
296 
297  if (get_a_e2_f(buf1, buf2, &e->a, &e->e2, &e->f) ||
298  get_a_e2_f(buf2, buf1, &e->a, &e->e2, &e->f))
299  table.count++;
300  else {
301  err++;
302  sprintf(buf, " %d", line);
303  if (*badlines)
304  strcat(badlines, ",");
305  strcat(badlines, buf);
306  continue;
307  }
308  }
309 
310  fclose(fd);
311 
312  if (!err) {
313  /* over correct typed version */
314  qsort(table.ellipses, table.count, sizeof(struct ellipse), compare_ellipse_names);
315  G_initialize_done(&table.initialized);
316  return 1;
317  }
318 
319  (fatal ? G_fatal_error : G_warning)(
320  n_(
321  ("Line%s of ellipsoid table file <%s> is invalid"),
322  ("Lines%s of ellipsoid table file <%s> are invalid"),
323  err),
324  badlines, file);
325 
326  G_initialize_done(&table.initialized);
327 
328  return 0;
329 }
330 
331 static int get_ellipsoid_parameters(struct Key_Value *proj_keys, double *a, double *e2)
332 {
333  const char *str, *str1;
334 
335  if (!proj_keys) {
336  return -1;
337  }
338 
339  if ((str = G_find_key_value("ellps", proj_keys)) != NULL) {
340  if (strncmp(str, "sphere", 6) == 0) {
341  str = G_find_key_value("a", proj_keys);
342  if (str != NULL) {
343  if (sscanf(str, "%lf", a) != 1)
344  G_fatal_error(_("Invalid a: field '%s' in file %s in <%s>"),
345  str, PROJECTION_FILE, PERMANENT);
346  }
347  else
348  *a = 6370997.0;
349 
350  *e2 = 0.0;
351 
352  return 0;
353  }
354  else {
355  if (G_get_ellipsoid_by_name(str, a, e2) == 0)
356  G_fatal_error(_("Invalid ellipsoid '%s' in file %s in <%s>"),
357  str, PROJECTION_FILE, PERMANENT);
358  else
359  return 1;
360  }
361  }
362  else {
363  str = G_find_key_value("a", proj_keys);
364  str1 = G_find_key_value("es", proj_keys);
365  if ((str != NULL) && (str1 != NULL)) {
366  if (sscanf(str, "%lf", a) != 1)
367  G_fatal_error(_("Invalid a: field '%s' in file %s in <%s>"),
368  str, PROJECTION_FILE, PERMANENT);
369  if (sscanf(str1, "%lf", e2) != 1)
370  G_fatal_error(_("Invalid es: field '%s' in file %s in <%s>"),
371  str, PROJECTION_FILE, PERMANENT);
372 
373  return 1;
374  }
375  else {
376  str = G_find_key_value("proj", proj_keys);
377  if ((str == NULL) || (strcmp(str, "ll") == 0)) {
378  *a = 6378137.0;
379  *e2 = .006694385;
380  return 0;
381  }
382  else
383  G_fatal_error(_("No ellipsoid info given in file %s in <%s>"),
385  }
386  }
387 
388  return 1;
389 }
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
Definition: strings.c:46
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
#define PERMANENT
Definition: get_projinfo.c:16
void G_strip(char *buf)
Removes all leading and trailing white space from string.
Definition: strings.c:258
#define n_(strs, strp, num)
Definition: glocale.h:14
int count
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition: get_ellipse.c:67
int G_is_initialized(int *p)
Definition: counter.c:59
#define NULL
Definition: ccmath.h:32
void G_initialize_done(int *p)
Definition: counter.c:76
const char * G_ellipsoid_name(int n)
Get ellipsoid name.
Definition: get_ellipse.c:133
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition: key_value1.c:103
fd
Definition: d/range.c:69
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
struct Key_Value * G_read_key_value_file(const char *file)
Read key/values pairs from file.
Definition: key_value3.c:53
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:38
int G_getl2(char *buf, int n, FILE *fd)
Gets a line of text from a file of any pedigree.
Definition: getl.c:64
double b
Definition: r_raster.c:39
fclose(fd)
#define GPATH_MAX
Definition: gis.h:151
const char * G_ellipsoid_description(int n)
Get description for nth ellipsoid.
Definition: get_ellipse.c:184
Definition: gis.h:479
int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
Get ellipsoid parameters by name.
Definition: get_ellipse.c:105
int G_get_spheroid_by_name(const char *name, double *a, double *e2, double *f)
Get spheroid parameters by name.
Definition: get_ellipse.c:154
#define file
#define _(str)
Definition: glocale.h:13
const char * name
Definition: named_colr.c:7
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int G_read_ellipsoid_table(int fatal)
Read ellipsoid table.
Definition: get_ellipse.c:247
#define PROJECTION_FILE
Definition: gis.h:103