1 /*!
2 \file lib/proj/ellipse.c
3
4 \brief GProj library - Functions for reading datum parameters from the location database
5
6 \author Paul Kelly <paul-grass stjohnspoint.co.uk>
7
8 (C) 2003-2008 by the GRASS Development Team
9
10 This program is free software under the GNU General Public
11 License (>=v2). Read the file COPYING that comes with GRASS
12 for details.
13 */
14
15 #include <unistd.h>
16 #include <ctype.h>
17 #include <string.h>
18 #include <stdlib.h>
19 #include <math.h> /* for sqrt() */
20 #include <grass/gis.h>
21 #include <grass/glocale.h>
22 #include <grass/gprojects.h>
23 #include "local_proto.h"
24
25 static int get_a_e2_rf(const char *, const char *, double *, double *,
26 double *);
27
28 /*!
29 * \brief Get the ellipsoid parameters from the database.
30 *
31 * If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
32 * that file, otherwise return WGS 84 values.
33 *
34 * Dies with diagnostic if there is an error.
35 *
36 * \param[out] a semi-major axis
37 * \param[out] e2 first eccentricity squared
38 * \param[out] rf reciprocal of the ellipsoid flattening term
39 *
40 * \return 1 on success
41 * \return 0 default values used.
42 */
GPJ_get_ellipsoid_params(double * a,double * e2,double * rf)43 int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
44 {
45 int ret;
46 struct Key_Value *proj_keys = G_get_projinfo();
47
48 if (proj_keys == NULL)
49 proj_keys = G_create_key_value();
50
51 ret = GPJ__get_ellipsoid_params(proj_keys, a, e2, rf);
52 G_free_key_value(proj_keys);
53
54 return ret;
55 }
56
57 /*!
58 * \brief Get the ellipsoid parameters from proj keys structure.
59 *
60 * If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
61 * that file, otherwise return WGS 84 values.
62 *
63 * Dies with diagnostic if there is an error.
64 *
65 * \param proj_keys proj definition
66 * \param[out] a semi-major axis
67 * \param[out] e2 first eccentricity squared
68 * \param[out] rf reciprocal of the ellipsoid flattening term
69 *
70 * \return 1 on success
71 * \return 0 default values used.
72 */
GPJ__get_ellipsoid_params(const struct Key_Value * proj_keys,double * a,double * e2,double * rf)73 int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys,
74 double *a, double *e2, double *rf)
75 {
76 struct gpj_ellps estruct;
77 struct gpj_datum dstruct;
78 const char *str, *str3;
79 char *str1, *ellps;
80
81 str = G_find_key_value("datum", proj_keys);
82
83 if ((str != NULL) && (GPJ_get_datum_by_name(str, &dstruct) > 0)) {
84 /* If 'datum' key is present, look up correct ellipsoid
85 * from datum.table */
86
87 ellps = G_store(dstruct.ellps);
88 GPJ_free_datum(&dstruct);
89
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 G_asprintf(&str4, "a=%s", str3);
115 if ((str3 = G_find_key_value("es", proj_keys)) != NULL)
116 G_asprintf(&str1, "e=%s", str3);
117 else if ((str3 = G_find_key_value("f", proj_keys)) != NULL)
118 G_asprintf(&str1, "f=1/%s", str3);
119 else if ((str3 = G_find_key_value("rf", proj_keys)) != NULL)
120 G_asprintf(&str1, "f=1/%s", str3);
121 else if ((str3 = G_find_key_value("b", proj_keys)) != NULL)
122 G_asprintf(&str1, "b=%s", str3);
123 else
124 G_fatal_error(_("No secondary ellipsoid descriptor "
125 "(rf, es or b) in file"));
126
127 if (get_a_e2_rf(str4, str1, a, e2, rf) == 0)
128 G_fatal_error(_("Invalid ellipsoid descriptors "
129 "(a, rf, es or b) in file"));
130 return 1;
131 }
132 else {
133 str = G_find_key_value("proj", proj_keys);
134 if ((str == NULL) || (strcmp(str, "ll") == 0)) {
135 *a = 6378137.0;
136 *e2 = .006694385;
137 *rf = 298.257223563;
138 return 0;
139 }
140 else {
141 G_fatal_error(_("No ellipsoid info given in file"));
142 }
143 }
144 }
145 return 1;
146 }
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
GPJ_get_ellipsoid_by_name(const char * name,struct gpj_ellps * estruct)160 int 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
get_a_e2_rf(const char * s1,const char * s2,double * a,double * e2,double * recipf)182 int 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
read_ellipsoid_table(int fatal)224 struct 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 int count = 0;
237
238 sprintf(file, "%s%s", G_gisbase(), ELLIPSOIDTABLE);
239 fd = fopen(file, "r");
240
241 if (!fd) {
242 (fatal ? G_fatal_error : G_warning)(
243 _("Unable to open ellipsoid table file <%s>"), file);
244 return NULL;
245 }
246
247 err = 0;
248 *badlines = 0;
249 for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
250 G_strip(buf);
251 if (*buf == 0 || *buf == '#')
252 continue;
253
254 if (sscanf(buf, "%s \"%1023[^\"]\" %s %s", name, descr, buf1, buf2)
255 != 4) {
256 err++;
257 sprintf(buf, " %d", line);
258 if (*badlines)
259 strcat(badlines, ",");
260 strcat(badlines, buf);
261 continue;
262 }
263
264
265 if (get_a_e2_rf(buf1, buf2, &a, &e2, &rf)
266 || get_a_e2_rf(buf2, buf1, &a, &e2, &rf)) {
267 if (current == NULL)
268 current = outputlist = G_malloc(sizeof(struct ellps_list));
269 else
270 current = current->next = G_malloc(sizeof(struct ellps_list));
271 current->name = G_store(name);
272 current->longname = G_store(descr);
273 current->a = a;
274 current->es = e2;
275 current->rf = rf;
276 current->next = NULL;
277 count++;
278 }
279 else {
280 err++;
281 sprintf(buf, " %d", line);
282 if (*badlines)
283 strcat(badlines, ",");
284 strcat(badlines, buf);
285 continue;
286 }
287 }
288
289 fclose(fd);
290
291 if (!err)
292 return outputlist;
293
294 (fatal ? G_fatal_error : G_warning)(
295 n_(
296 ("Line%s of ellipsoid table file <%s> is invalid"),
297 ("Lines%s of ellipsoid table file <%s> are invalid"),
298 err),
299 badlines, file);
300
301 return outputlist;
302 }
303
304 /*!
305 \brief Free ellipsoid data structure.
306
307 \param estruct data structure to be freed
308 */
GPJ_free_ellps(struct gpj_ellps * estruct)309 void GPJ_free_ellps(struct gpj_ellps *estruct)
310 {
311 G_free(estruct->name);
312 G_free(estruct->longname);
313 return;
314 }
315
free_ellps_list(struct ellps_list * elist)316 void free_ellps_list(struct ellps_list *elist)
317 {
318 struct ellps_list *old;
319
320 while (elist != NULL) {
321 G_free(elist->name);
322 G_free(elist->longname);
323 old = elist;
324 elist = old->next;
325 G_free(old);
326 }
327
328 return;
329 }
330