GRASS GIS 7 Programmer's Manual  7.0.3(2016)-r00000
ellipse.c
Go to the documentation of this file.
1 
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 
26 static int get_a_e2_rf(const char *, const char *, double *, double *,
27  double *);
28 
37 int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
38 {
39  int ret;
40  struct Key_Value *proj_keys = G_get_projinfo();
41 
42  if (proj_keys == NULL)
43  proj_keys = G_create_key_value();
44 
45  ret = GPJ__get_ellipsoid_params(proj_keys, a, e2, rf);
46  G_free_key_value(proj_keys);
47 
48  return ret;
49 }
50 
51 int
52 GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys,
53  double *a, double *e2, double *rf)
54 {
55  struct gpj_ellps estruct;
56  struct gpj_datum dstruct;
57  const char *str, *str3;
58  char *str1, *ellps;
59 
60  str = G_find_key_value("datum", proj_keys);
61 
62  if ((str != NULL) && (GPJ_get_datum_by_name(str, &dstruct) > 0)) {
63  /* If 'datum' key is present, look up correct ellipsoid
64  * from datum.table */
65 
66  ellps = G_store(dstruct.ellps);
67  GPJ_free_datum(&dstruct);
68 
69  }
70  else
71  /* else use ellipsoid defined in PROJ_INFO */
72  ellps = G_store(G_find_key_value("ellps", proj_keys));
73 
74  if (ellps != NULL && *ellps) {
75  if (GPJ_get_ellipsoid_by_name(ellps, &estruct) < 0)
76  G_fatal_error(_("Invalid ellipsoid <%s> in file"), ellps);
77 
78  *a = estruct.a;
79  *e2 = estruct.es;
80  *rf = estruct.rf;
81  GPJ_free_ellps(&estruct);
82  G_free(ellps);
83 
84  return 1;
85  }
86  else {
87  if (ellps) /* *ellps = '\0' */
88  G_free(ellps);
89 
90  str3 = G_find_key_value("a", proj_keys);
91  if (str3 != NULL) {
92  char *str4;
93  G_asprintf(&str4, "a=%s", str3);
94  if ((str3 = G_find_key_value("es", proj_keys)) != NULL)
95  G_asprintf(&str1, "e=%s", str3);
96  else if ((str3 = G_find_key_value("f", proj_keys)) != NULL)
97  G_asprintf(&str1, "f=1/%s", str3);
98  else if ((str3 = G_find_key_value("rf", proj_keys)) != NULL)
99  G_asprintf(&str1, "f=1/%s", str3);
100  else if ((str3 = G_find_key_value("b", proj_keys)) != NULL)
101  G_asprintf(&str1, "b=%s", str3);
102  else
103  G_fatal_error(_("No secondary ellipsoid descriptor "
104  "(rf, es or b) in file"));
105 
106  if (get_a_e2_rf(str4, str1, a, e2, rf) == 0)
107  G_fatal_error(_("Invalid ellipsoid descriptors "
108  "(a, rf, es or b) in file"));
109  return 1;
110  }
111  else {
112  str = G_find_key_value("proj", proj_keys);
113  if ((str == NULL) || (strcmp(str, "ll") == 0)) {
114  *a = 6378137.0;
115  *e2 = .006694385;
116  *rf = 298.257223563;
117  return 0;
118  }
119  else {
120  G_fatal_error(_("No ellipsoid info given in file"));
121  }
122  }
123  }
124  return 1;
125 }
126 
127 
136 int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
137 {
138  struct ellps_list *list, *listhead;
139 
140  list = listhead = read_ellipsoid_table(0);
141 
142  while (list != NULL) {
143  if (G_strcasecmp(name, list->name) == 0) {
144  estruct->name = G_store(list->name);
145  estruct->longname = G_store(list->longname);
146  estruct->a = list->a;
147  estruct->es = list->es;
148  estruct->rf = list->rf;
149  free_ellps_list(listhead);
150  return 1;
151  }
152  list = list->next;
153  }
154  free_ellps_list(listhead);
155  return -1;
156 }
157 
158 static int
159 get_a_e2_rf(const char *s1, const char *s2, double *a, double *e2,
160  double *recipf)
161 {
162  double b, f;
163 
164  if (sscanf(s1, "a=%lf", a) != 1)
165  return 0;
166 
167  if (*a <= 0.0)
168  return 0;
169 
170  if (sscanf(s2, "e=%lf", e2) == 1) {
171  f = 1.0 - sqrt(1.0 - *e2);
172  *recipf = 1.0 / f;
173  return (*e2 >= 0.0);
174  }
175 
176  if (sscanf(s2, "f=1/%lf", recipf) == 1) {
177  if (*recipf <= 0.0)
178  return 0;
179  f = 1.0 / *recipf;
180  *e2 = f * (2 - f);
181  return (*e2 >= 0.0);
182  }
183 
184  if (sscanf(s2, "b=%lf", &b) == 1) {
185  if (b <= 0.0)
186  return 0;
187  if (b == *a) {
188  f = 0.0;
189  *e2 = 0.0;
190  }
191  else {
192  f = (*a - b) / *a;
193  *e2 = f * (2 - f);
194  }
195  *recipf = 1.0 / f;
196  return (*e2 >= 0.0);
197  }
198  return 0;
199 }
200 
201 struct ellps_list *read_ellipsoid_table(int fatal)
202 {
203  FILE *fd;
204  char file[GPATH_MAX];
205  char buf[4096];
206  char name[100], descr[1024], buf1[1024], buf2[1024];
207  char badlines[1024];
208  int line;
209  int err;
210  struct ellps_list *current = NULL, *outputlist = NULL;
211  double a, e2, rf;
212 
213  int count = 0;
214 
215  sprintf(file, "%s%s", G_gisbase(), ELLIPSOIDTABLE);
216  fd = fopen(file, "r");
217 
218  if (!fd) {
219  (fatal ? G_fatal_error : G_warning)(
220  _("Unable to open ellipsoid table file <%s>"), file);
221  return NULL;
222  }
223 
224  err = 0;
225  *badlines = 0;
226  for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
227  G_strip(buf);
228  if (*buf == 0 || *buf == '#')
229  continue;
230 
231  if (sscanf(buf, "%s \"%1023[^\"]\" %s %s", name, descr, buf1, buf2)
232  != 4) {
233  err++;
234  sprintf(buf, " %d", line);
235  if (*badlines)
236  strcat(badlines, ",");
237  strcat(badlines, buf);
238  continue;
239  }
240 
241 
242  if (get_a_e2_rf(buf1, buf2, &a, &e2, &rf)
243  || get_a_e2_rf(buf2, buf1, &a, &e2, &rf)) {
244  if (current == NULL)
245  current = outputlist = G_malloc(sizeof(struct ellps_list));
246  else
247  current = current->next = G_malloc(sizeof(struct ellps_list));
248  current->name = G_store(name);
249  current->longname = G_store(descr);
250  current->a = a;
251  current->es = e2;
252  current->rf = rf;
253  current->next = NULL;
254  count++;
255  }
256  else {
257  err++;
258  sprintf(buf, " %d", line);
259  if (*badlines)
260  strcat(badlines, ",");
261  strcat(badlines, buf);
262  continue;
263  }
264  }
265 
266  fclose(fd);
267 
268  if (!err)
269  return outputlist;
270 
271  (fatal ? G_fatal_error : G_warning)(
272  n_(
273  ("Line%s of ellipsoid table file <%s> is invalid"),
274  ("Lines%s of ellipsoid table file <%s> are invalid"),
275  err),
276  badlines, file);
277 
278  return outputlist;
279 }
280 
281 void GPJ_free_ellps(struct gpj_ellps *estruct)
282 {
283  G_free(estruct->name);
284  G_free(estruct->longname);
285  return;
286 }
287 
288 void free_ellps_list(struct ellps_list *elist)
289 {
290  struct ellps_list *old;
291 
292  while (elist != NULL) {
293  G_free(elist->name);
294  G_free(elist->longname);
295  old = elist;
296  elist = old->next;
297  G_free(old);
298  }
299 
300  return;
301 }
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
Definition: strings.c:46
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
void GPJ_free_datum(struct gpj_datum *dstruct)
Free the memory used for the strings in a gpj_datum struct.
Definition: proj/datum.c:395
void G_strip(char *buf)
Removes all leading and trailing white space from string.
Definition: strings.c:258
struct ellps_list * read_ellipsoid_table(int fatal)
Definition: ellipse.c:201
int count
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
int G_asprintf(char **out, const char *fmt,...)
Definition: asprintf.c:70
#define NULL
Definition: ccmath.h:32
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:37
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition: key_value1.c:103
int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a, double *e2, double *rf)
Definition: ellipse.c:52
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
double b
int G_getl2(char *buf, int n, FILE *fd)
Gets a line of text from a file of any pedigree.
Definition: getl.c:64
struct list * list
Definition: read_list.c:24
void GPJ_free_ellps(struct gpj_ellps *estruct)
Definition: ellipse.c:281
struct Key_Value * G_get_projinfo(void)
Gets projection information for location.
Definition: get_projinfo.c:58
#define file
int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
Definition: ellipse.c:37
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:136
const char * name
Definition: named_colr.c:7
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition: key_value1.c:23
void free_ellps_list(struct ellps_list *elist)
Definition: ellipse.c:288