GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 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  */
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  */
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 
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 
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 
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 */
309 void GPJ_free_ellps(struct gpj_ellps *estruct)
310 {
311  G_free(estruct->name);
312  G_free(estruct->longname);
313  return;
314 }
315 
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 }
int G_getl2(char *, int, FILE *)
Gets a line of text from a file of any pedigree.
Definition: getl.c:64
#define G_malloc(n)
Definition: defs/gis.h:112
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
Definition: get_projinfo.c:61
#define n_(strs, strp, num)
Definition: glocale.h:11
struct ellps_list * read_ellipsoid_table(int fatal)
Definition: ellipse.c:224
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int count
double a
Definition: gprojects.h:110
void G_strip(char *)
Removes all leading and trailing white space from string.
Definition: strings.c:300
#define NULL
Definition: ccmath.h:32
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:73
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
int GPJ_get_datum_by_name(const char *, struct gpj_datum *)
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:37
int int G_strcasecmp(const char *, const char *)
String compare ignoring case (upper or lower)
Definition: strings.c:47
double b
Definition: r_raster.c:39
struct list * list
Definition: read_list.c:24
void GPJ_free_ellps(struct gpj_ellps *estruct)
Free ellipsoid data structure.
Definition: ellipse.c:309
char * ellps
Definition: gprojects.h:86
char * name
Definition: gprojects.h:109
#define GPATH_MAX
Definition: gis.h:170
void GPJ_free_datum(struct gpj_datum *)
Free the memory used for the strings in a gpj_datum struct.
Definition: proj/datum.c:402
Definition: gis.h:501
char * longname
Definition: gprojects.h:109
double es
Definition: gprojects.h:110
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
void G_warning(const char *,...) __attribute__((format(printf
#define file
int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
Get the ellipsoid parameters from the database.
Definition: ellipse.c:43
#define _(str)
Definition: glocale.h:10
char * G_store(const char *)
Copy string to allocated memory.
Definition: strings.c:87
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition: key_value1.c:23
int G_asprintf(char **, const char *,...) __attribute__((format(printf
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
const char * name
Definition: named_colr.c:7
#define ELLIPSOIDTABLE
Definition: gprojects.h:62
double rf
Definition: gprojects.h:110
void G_free_key_value(struct Key_Value *)
Free allocated Key_Value structure.
Definition: key_value1.c:103
const char * G_find_key_value(const char *, const struct Key_Value *)
Find given key (case sensitive)
Definition: key_value1.c:84
void free_ellps_list(struct ellps_list *elist)
Definition: ellipse.c:316