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