GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
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
32static const char PERMANENT[] = "PERMANENT";
33
34static 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 *); */
49static int get_a_e2_f(const char *, const char *, double *, double *, double *);
50static int compare_ellipse_names(const void *, const void *);
51static 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 */
67int 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
75 G_file_name(ipath, "", PROJECTION_FILE, PERMANENT);
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 */
105int 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 */
133const 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 */
154int 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 */
184const char *G_ellipsoid_description(int n)
185{
187 return n >= 0 && n < table.count ? table.ellipses[n].descr : NULL;
188}
189
190static 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
230static 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
331static 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>"),
384 PROJECTION_FILE, PERMANENT);
385 }
386 }
387
388 return 1;
389}
#define NULL
Definition: ccmath.h:32
void G_initialize_done(int *p)
Definition: counter.c:76
int G_is_initialized(int *p)
Definition: counter.c:59
double b
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
int G_read_ellipsoid_table(int fatal)
Read ellipsoid table.
Definition: get_ellipse.c:247
const char * G_ellipsoid_description(int n)
Get description for nth ellipsoid.
Definition: get_ellipse.c:184
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
const char * G_ellipsoid_name(int n)
Get ellipsoid name.
Definition: get_ellipse.c:133
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition: get_ellipse.c:67
int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
Get ellipsoid parameters by name.
Definition: get_ellipse.c:105
#define PERMANENT
Definition: get_projinfo.c:19
int G_getl2(char *buf, int n, FILE *fd)
Gets a line of text from a file of any pedigree.
Definition: getl.c:64
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
int count
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition: key_value1.c:103
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
struct Key_Value * G_read_key_value_file(const char *file)
Read key/values pairs from file.
Definition: key_value3.c:53
#define file
const char * name
Definition: named_colr.c:7
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
Definition: strings.c:47
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:87
void G_strip(char *buf)
Removes all leading and trailing white space from string.
Definition: strings.c:300
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220