GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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 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 *); */
48static int get_a_e2_f(const char *, const char *, double *, double *, double *);
49static int compare_ellipse_names(const void *, const void *);
50static 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 */
66int 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
74 G_file_name(ipath, "", PROJECTION_FILE, PERMANENT);
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 */
104int 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 */
132const 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 */
153int 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 */
182const char *G_ellipsoid_description(int n)
183{
185 return n >= 0 && n < table.count ? table.ellipses[n].descr : NULL;
186}
187
188static 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
229static 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
332static 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
void G_initialize_done(int *p)
Definition counter.c:77
int G_is_initialized(int *p)
Definition counter.c:60
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.
const char * G_ellipsoid_description(int n)
Get description for nth ellipsoid.
int G_get_spheroid_by_name(const char *name, double *a, double *e2, double *f)
Get spheroid parameters by name.
const char * G_ellipsoid_name(int n)
Get ellipsoid name.
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.
#define PERMANENT
int G_getl2(char *buf, int n, FILE *fd)
Gets a line of text from a file of any pedigree.
Definition getl.c:65
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition gisbase.c:39
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition key_value1.c:104
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition key_value1.c:85
struct Key_Value * G_read_key_value_file(const char *file)
Read key/values pairs from file.
Definition key_value3.c:55
#define file
const char * name
Definition named_colr.c:6
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)