00001
00023 #include <unistd.h>
00024 #include <ctype.h>
00025 #include <string.h>
00026 #include <stdlib.h>
00027 #include <math.h>
00028 #include <grass/gis.h>
00029 #include <grass/glocale.h>
00030
00031 static struct table
00032 {
00033 char *name;
00034 char *descr;
00035 double a;
00036 double e2;
00037 double f;
00038 } *table = NULL;
00039
00040 static int count = -1;
00041 static char *PERMANENT = "PERMANENT";
00042
00043
00044 static int get_a_e2_f(const char *, const char *, double *, double *,
00045 double *);
00046 void ellipsoid_table_file(char *);
00047 static int compare_table_names(const void *, const void *);
00048 static int read_ellipsoid_table(int);
00049 static int get_ellipsoid_parameters(struct Key_Value *, double *, double *);
00050
00051
00066 int G_get_ellipsoid_parameters(double *a, double *e2)
00067 {
00068 int in_stat, stat;
00069 char ipath[GPATH_MAX];
00070 struct Key_Value *proj_keys;
00071
00072 proj_keys = NULL;
00073
00074 G__file_name(ipath, "", PROJECTION_FILE, PERMANENT);
00075
00076 if (access(ipath, 0) != 0) {
00077 *a = 6378137.0;
00078 *e2 = .006694385;
00079 return 0;
00080 }
00081
00082 proj_keys = G_read_key_value_file(ipath, &in_stat);
00083
00084 if (in_stat != 0) {
00085 G_fatal_error(_("Unable to open file %s in <%s>"),
00086 PROJECTION_FILE, PERMANENT);
00087 }
00088
00089 stat = get_ellipsoid_parameters(proj_keys, a, e2);
00090
00091 G_free_key_value(proj_keys);
00092
00093 return stat;
00094 }
00095
00110 int G_get_ellipsoid_by_name(const char *name, double *a, double *e2)
00111 {
00112 int i;
00113
00114 (void)read_ellipsoid_table(0);
00115
00116 for (i = 0; i < count; i++) {
00117 if (G_strcasecmp(name, table[i].name) == 0) {
00118 *a = table[i].a;
00119 *e2 = table[i].e2;
00120 return 1;
00121 }
00122 }
00123 return 0;
00124 }
00125
00138 char *G_ellipsoid_name(int n)
00139 {
00140 (void)read_ellipsoid_table(0);
00141 return n >= 0 && n < count ? table[n].name : NULL;
00142 }
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00171 int G_get_spheroid_by_name(const char *name, double *a, double *e2, double *f)
00172 {
00173 int i;
00174
00175 (void)read_ellipsoid_table(0);
00176
00177 for (i = 0; i < count; i++) {
00178 if (G_strcasecmp(name, table[i].name) == 0) {
00179 *a = table[i].a;
00180 *e2 = table[i].e2;
00181 *f = table[i].f;
00182 return 1;
00183 }
00184 }
00185 return 0;
00186 }
00187
00188
00202 char *G_ellipsoid_description(int n)
00203 {
00204 (void)read_ellipsoid_table(0);
00205 return n >= 0 && n < count ? table[n].descr : NULL;
00206 }
00207
00208 static int
00209 get_a_e2_f(const char *s1, const char *s2, double *a, double *e2, double *f)
00210 {
00211 double b, recipf;
00212
00213 if (sscanf(s1, "a=%lf", a) != 1)
00214 return 0;
00215
00216 if (*a <= 0.0)
00217 return 0;
00218
00219 if (sscanf(s2, "e=%lf", e2) == 1) {
00220 *f = (double)1.0 / -sqrt(((double)1.0 - *e2)) + (double)1.0;
00221 return (*e2 >= 0.0);
00222 }
00223
00224 if (sscanf(s2, "f=1/%lf", f) == 1) {
00225 if (*f <= 0.0)
00226 return 0;
00227 recipf = (double)1.0 / (*f);
00228 *e2 = recipf + recipf - recipf * recipf;
00229 return (*e2 >= 0.0);
00230 }
00231
00232 if (sscanf(s2, "b=%lf", &b) == 1) {
00233 if (b <= 0.0)
00234 return 0;
00235 if (b == *a) {
00236 *f = 0.0;
00237 *e2 = 0.0;
00238 }
00239 else {
00240 recipf = ((*a) - b) / (*a);
00241 *f = (double)1.0 / recipf;
00242 *e2 = recipf + recipf - recipf * recipf;
00243 }
00244 return (*e2 >= 0.0);
00245 }
00246 return 0;
00247 }
00248
00249 void ellipsoid_table_file(char *file)
00250 {
00251 sprintf(file, "%s/etc/ellipse.table", G_gisbase());
00252 return;
00253 }
00254
00255 static int compare_table_names(const void *pa, const void *pb)
00256 {
00257 const struct table *a = pa, *b = pb;
00258
00259
00260 return G_strcasecmp(a->name, b->name);
00261 }
00262
00263 static int read_ellipsoid_table(int fatal)
00264 {
00265 FILE *fd;
00266 char file[GPATH_MAX];
00267 char buf[1024];
00268 char name[100], descr[100], buf1[100], buf2[100];
00269 char badlines[256];
00270 int line;
00271 int err;
00272
00273 if (count >= 0)
00274 return 1;
00275 count = 0;
00276 table = NULL;
00277
00278 (void)ellipsoid_table_file(file);
00279 fd = fopen(file, "r");
00280
00281 if (fd == NULL) {
00282 perror(file);
00283 sprintf(buf, _("Unable to open ellipsoid table file <%s>"), file);
00284 fatal ? G_fatal_error(buf) : G_warning(buf);
00285 return 0;
00286 }
00287
00288 err = 0;
00289 *badlines = 0;
00290 for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
00291 G_strip(buf);
00292 if (*buf == 0 || *buf == '#')
00293 continue;
00294
00295 if (sscanf(buf, "%s \"%99[^\"]\" %s %s", name, descr, buf1, buf2) !=
00296 4) {
00297 err++;
00298 sprintf(buf, " %d", line);
00299 if (*badlines)
00300 G_strcat(badlines, ",");
00301 G_strcat(badlines, buf);
00302 continue;
00303 }
00304
00305 table =
00306 (struct table *)G_realloc((char *)table,
00307 (count + 1) * sizeof(*table));
00308 table[count].name = G_store(name);
00309 table[count].descr = G_store(descr);
00310
00311 if (get_a_e2_f
00312 (buf1, buf2, &table[count].a, &table[count].e2, &table[count].f)
00313 || get_a_e2_f(buf2, buf1, &table[count].a, &table[count].e2,
00314 &table[count].f))
00315 count++;
00316 else {
00317 err++;
00318 sprintf(buf, " %d", line);
00319 if (*badlines)
00320 G_strcat(badlines, ",");
00321 G_strcat(badlines, buf);
00322 continue;
00323 }
00324 }
00325
00326 fclose(fd);
00327
00328 if (!err) {
00329
00330 qsort(table, count, sizeof(*table), compare_table_names);
00331 return 1;
00332 }
00333
00334 (fatal ? G_fatal_error : G_warning) ((err > 1)
00335 ?
00336 _("Lines%s of ellipsoid table file <%s> are invalid")
00337 :
00338 _("Line%s of ellipsoid table file <%s> is invalid"),
00339 badlines, file);
00340
00341 return 0;
00342 }
00343
00344 static int get_ellipsoid_parameters(struct Key_Value *proj_keys, double *a,
00345 double *e2)
00346 {
00347 char *str, *str1;
00348
00349 if (!proj_keys) {
00350 return -1;
00351 }
00352
00353 if ((str = G_find_key_value("ellps", proj_keys)) != NULL) {
00354 if (strncmp(str, "sphere", 6) == 0) {
00355 str = G_find_key_value("a", proj_keys);
00356 if (str != NULL) {
00357 if (sscanf(str, "%lf", a) != 1) {
00358 G_fatal_error(_("Invalid a: field '%s' in file %s in <%s>"),
00359 str, PROJECTION_FILE, PERMANENT);
00360 }
00361 }
00362 else {
00363 *a = 6370997.0;
00364 }
00365 *e2 = 0.0;
00366
00367 return 0;
00368 }
00369 else {
00370 if (G_get_ellipsoid_by_name(str, a, e2) == 0) {
00371 G_fatal_error(_("Invalid ellipsoid '%s' in file %s in <%s>"),
00372 str, PROJECTION_FILE, PERMANENT);
00373 }
00374 else {
00375 return 1;
00376 }
00377 }
00378 }
00379 else {
00380 str = G_find_key_value("a", proj_keys);
00381 str1 = G_find_key_value("es", proj_keys);
00382 if ((str != NULL) && (str1 != NULL)) {
00383 if (sscanf(str, "%lf", a) != 1) {
00384 G_fatal_error(_("Invalid a: field '%s' in file %s in <%s>"),
00385 str, PROJECTION_FILE, PERMANENT);
00386 }
00387 if (sscanf(str1, "%lf", e2) != 1) {
00388 G_fatal_error(_("Invalid es: field '%s' in file %s in <%s>"),
00389 str, PROJECTION_FILE, PERMANENT);
00390 }
00391
00392 return 1;
00393 }
00394 else {
00395 str = G_find_key_value("proj", proj_keys);
00396 if ((str == NULL) || (strcmp(str, "ll") == 0)) {
00397 *a = 6378137.0;
00398 *e2 = .006694385;
00399 return 0;
00400 }
00401 else {
00402 G_fatal_error(_("No ellipsoid info given in file %s in <%s>"),
00403 PROJECTION_FILE, PERMANENT);
00404 }
00405 }
00406 }
00407
00408 return 1;
00409 }