GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
ellipse.c
Go to the documentation of this file.
1/*!
2 \file lib/proj/ellipse.c
3
4 \brief GProj library - Functions for reading datum parameters from the
5 location database
6
7 \author Paul Kelly <paul-grass stjohnspoint.co.uk>
8
9 (C) 2003-2008 by the GRASS Development Team
10
11 This program is free software under the GNU General Public
12 License (>=v2). Read the file COPYING that comes with GRASS
13 for details.
14 */
15
16#include <unistd.h>
17#include <ctype.h>
18#include <string.h>
19#include <stdlib.h>
20#include <math.h> /* for sqrt() */
21#include <grass/gis.h>
22#include <grass/glocale.h>
23#include <grass/gprojects.h>
24#include "local_proto.h"
25
26static int get_a_e2_rf(const char *, const char *, double *, double *,
27 double *);
28
29/*!
30 * \brief Get the ellipsoid parameters from the database.
31 *
32 * If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
33 * that file, otherwise return WGS 84 values.
34 *
35 * Dies with diagnostic if there is an error.
36 *
37 * \param[out] a semi-major axis
38 * \param[out] e2 first eccentricity squared
39 * \param[out] rf reciprocal of the ellipsoid flattening term
40 *
41 * \return 1 on success
42 * \return 0 default values used.
43 */
44int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
45{
46 int ret;
47 struct Key_Value *proj_keys = G_get_projinfo();
48
49 if (proj_keys == NULL)
50 proj_keys = G_create_key_value();
51
52 ret = GPJ__get_ellipsoid_params(proj_keys, a, e2, rf);
53 G_free_key_value(proj_keys);
54
55 return ret;
56}
57
58/*!
59 * \brief Get the ellipsoid parameters from proj keys structure.
60 *
61 * If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
62 * that file, otherwise return WGS 84 values.
63 *
64 * Dies with diagnostic if there is an error.
65 *
66 * \param proj_keys proj definition
67 * \param[out] a semi-major axis
68 * \param[out] e2 first eccentricity squared
69 * \param[out] rf reciprocal of the ellipsoid flattening term
70 *
71 * \return 1 on success
72 * \return 0 default values used.
73 */
74int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a,
75 double *e2, double *rf)
76{
77 struct gpj_ellps estruct;
78 struct gpj_datum dstruct;
79 const char *str, *str3;
80 char *str1, *ellps;
81
82 str = G_find_key_value("datum", proj_keys);
83
84 if ((str != NULL) && (GPJ_get_datum_by_name(str, &dstruct) > 0)) {
85 /* If 'datum' key is present, look up correct ellipsoid
86 * from datum.table */
87
88 ellps = G_store(dstruct.ellps);
89 GPJ_free_datum(&dstruct);
90 }
91 else
92 /* else use ellipsoid defined in PROJ_INFO */
93 ellps = G_store(G_find_key_value("ellps", proj_keys));
94
95 if (ellps != NULL && *ellps) {
96 if (GPJ_get_ellipsoid_by_name(ellps, &estruct) < 0)
97 G_fatal_error(_("Invalid ellipsoid <%s> in file"), ellps);
98
99 *a = estruct.a;
100 *e2 = estruct.es;
101 *rf = estruct.rf;
102 GPJ_free_ellps(&estruct);
103 G_free(ellps);
104
105 return 1;
106 }
107 else {
108 if (ellps) /* *ellps = '\0' */
109 G_free(ellps);
110
111 str3 = G_find_key_value("a", proj_keys);
112 if (str3 != NULL) {
113 char *str4;
114
115 G_asprintf(&str4, "a=%s", str3);
116 if ((str3 = G_find_key_value("es", proj_keys)) != NULL)
117 G_asprintf(&str1, "e=%s", str3);
118 else if ((str3 = G_find_key_value("f", proj_keys)) != NULL)
119 G_asprintf(&str1, "f=1/%s", str3);
120 else if ((str3 = G_find_key_value("rf", proj_keys)) != NULL)
121 G_asprintf(&str1, "f=1/%s", str3);
122 else if ((str3 = G_find_key_value("b", proj_keys)) != NULL)
123 G_asprintf(&str1, "b=%s", str3);
124 else
125 G_fatal_error(_("No secondary ellipsoid descriptor "
126 "(rf, es or b) in file"));
127
128 if (get_a_e2_rf(str4, str1, a, e2, rf) == 0)
129 G_fatal_error(_("Invalid ellipsoid descriptors "
130 "(a, rf, es or b) in file"));
131 return 1;
132 }
133 else {
134 str = G_find_key_value("proj", proj_keys);
135 if ((str == NULL) || (strcmp(str, "ll") == 0)) {
136 *a = 6378137.0;
137 *e2 = .006694385;
138 *rf = 298.257223563;
139 return 0;
140 }
141 else {
142 G_fatal_error(_("No ellipsoid info given in file"));
143 }
144 }
145 }
146 return 1;
147}
148
149/*!
150 * \brief Looks up ellipsoid in ellipsoid table and returns the a, e2
151 * parameters for the ellipsoid.
152 *
153 * \param name ellipsoid name
154 * \param[out] estruct ellipsoid
155 *
156 * \return 1 on success
157 * \return -1 if not found in table
158 */
159
160int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
161{
162 struct ellps_list *list, *listhead;
163
164 list = listhead = read_ellipsoid_table(0);
165
166 while (list != NULL) {
167 if (G_strcasecmp(name, list->name) == 0) {
168 estruct->name = G_store(list->name);
169 estruct->longname = G_store(list->longname);
170 estruct->a = list->a;
171 estruct->es = list->es;
172 estruct->rf = list->rf;
173 free_ellps_list(listhead);
174 return 1;
175 }
176 list = list->next;
177 }
178 free_ellps_list(listhead);
179 return -1;
180}
181
182int get_a_e2_rf(const char *s1, const char *s2, double *a, double *e2,
183 double *recipf)
184{
185 double b, f;
186
187 if (sscanf(s1, "a=%lf", a) != 1)
188 return 0;
189
190 if (*a <= 0.0)
191 return 0;
192
193 if (sscanf(s2, "e=%lf", e2) == 1) {
194 f = 1.0 - sqrt(1.0 - *e2);
195 *recipf = 1.0 / f;
196 return (*e2 >= 0.0);
197 }
198
199 if (sscanf(s2, "f=1/%lf", recipf) == 1) {
200 if (*recipf <= 0.0)
201 return 0;
202 f = 1.0 / *recipf;
203 *e2 = f * (2 - f);
204 return (*e2 >= 0.0);
205 }
206
207 if (sscanf(s2, "b=%lf", &b) == 1) {
208 if (b <= 0.0)
209 return 0;
210 if (b == *a) {
211 f = 0.0;
212 *e2 = 0.0;
213 }
214 else {
215 f = (*a - b) / *a;
216 *e2 = f * (2 - f);
217 }
218 *recipf = 1.0 / f;
219 return (*e2 >= 0.0);
220 }
221 return 0;
222}
223
224struct ellps_list *read_ellipsoid_table(int fatal)
225{
226 FILE *fd;
227 char file[GPATH_MAX];
228 char buf[4096];
229 char name[100], descr[1024], buf1[1024], buf2[1024];
230 char badlines[1024];
231 int line;
232 int err;
233 struct ellps_list *current = NULL, *outputlist = NULL;
234 double a, e2, rf;
235
236 sprintf(file, "%s%s", G_gisbase(), ELLIPSOIDTABLE);
237 fd = fopen(file, "r");
238
239 if (!fd) {
240 (fatal ? G_fatal_error : G_warning)(
241 _("Unable to open ellipsoid table file <%s>"), file);
242 return NULL;
243 }
244
245 err = 0;
246 *badlines = 0;
247 for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
248 G_strip(buf);
249 if (*buf == 0 || *buf == '#')
250 continue;
251
252 if (sscanf(buf, "%s \"%1023[^\"]\" %s %s", name, descr, buf1, buf2) !=
253 4) {
254 err++;
255 sprintf(buf, " %d", line);
256 if (*badlines)
257 strcat(badlines, ",");
258 strcat(badlines, buf);
259 continue;
260 }
261
262 if (get_a_e2_rf(buf1, buf2, &a, &e2, &rf) ||
263 get_a_e2_rf(buf2, buf1, &a, &e2, &rf)) {
264 if (current == NULL)
265 current = outputlist = G_malloc(sizeof(struct ellps_list));
266 else
267 current = current->next = G_malloc(sizeof(struct ellps_list));
268 current->name = G_store(name);
269 current->longname = G_store(descr);
270 current->a = a;
271 current->es = e2;
272 current->rf = rf;
273 current->next = NULL;
274 }
275 else {
276 err++;
277 sprintf(buf, " %d", line);
278 if (*badlines)
279 strcat(badlines, ",");
280 strcat(badlines, buf);
281 continue;
282 }
283 }
284
285 fclose(fd);
286
287 if (!err)
288 return outputlist;
289
290 (fatal ? G_fatal_error : G_warning)(
291 n_(("Line%s of ellipsoid table file <%s> is invalid"),
292 ("Lines%s of ellipsoid table file <%s> are invalid"), err),
293 badlines, file);
294
295 return outputlist;
296}
297
298/*!
299 \brief Free ellipsoid data structure.
300
301 \param estruct data structure to be freed
302 */
303void GPJ_free_ellps(struct gpj_ellps *estruct)
304{
305 G_free(estruct->name);
306 G_free(estruct->longname);
307 return;
308}
309
310void free_ellps_list(struct ellps_list *elist)
311{
312 struct ellps_list *old;
313
314 while (elist != NULL) {
315 G_free(elist->name);
316 G_free(elist->longname);
317 old = elist;
318 elist = old->next;
319 G_free(old);
320 }
321
322 return;
323}
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
int G_asprintf(char **out, const char *fmt,...)
Definition asprintf.c:69
#define NULL
Definition ccmath.h:32
double b
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:160
struct ellps_list * read_ellipsoid_table(int fatal)
Definition ellipse.c:224
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
void free_ellps_list(struct ellps_list *elist)
Definition ellipse.c:310
void GPJ_free_ellps(struct gpj_ellps *estruct)
Free ellipsoid data structure.
Definition ellipse.c:303
int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
Get the ellipsoid parameters from the database.
Definition ellipse.c:44
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
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_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition key_value1.c:23
#define file
const char * name
Definition named_colr.c:6
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
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
struct list * list
Definition read_list.c:24
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)