GRASS Programmer's Manual  6.4.2(2012)
ellipse.c
Go to the documentation of this file.
00001 
00016 #include <unistd.h>
00017 #include <ctype.h>
00018 #include <string.h>
00019 #include <stdlib.h>
00020 #include <math.h>               /* for sqrt() */
00021 #include <grass/gis.h>
00022 #include <grass/glocale.h>
00023 #include <grass/gprojects.h>
00024 #include "local_proto.h"
00025 
00026 static int get_a_e2_rf(const char *, const char *, double *, double *,
00027                        double *);
00028 
00037 int GPJ_get_ellipsoid_params(double *a, double *e2, double *rf)
00038 {
00039     int ret;
00040     struct Key_Value *proj_keys = G_get_projinfo();
00041 
00042     if (proj_keys == NULL)
00043         proj_keys = G_create_key_value();
00044 
00045     ret = GPJ__get_ellipsoid_params(proj_keys, a, e2, rf);
00046     G_free_key_value(proj_keys);
00047 
00048     return ret;
00049 }
00050 
00051 int
00052 GPJ__get_ellipsoid_params(struct Key_Value *proj_keys,
00053                           double *a, double *e2, double *rf)
00054 {
00055     struct gpj_ellps estruct;
00056     struct gpj_datum dstruct;
00057     char *str, *str1, *str3;
00058 
00059     str = G_find_key_value("datum", proj_keys);
00060 
00061     if ((str != NULL) && (GPJ_get_datum_by_name(str, &dstruct) > 0)) {
00062         /* If 'datum' key is present, look up correct ellipsoid
00063          * from datum.table */
00064 
00065         str = G_store(dstruct.ellps);
00066         GPJ_free_datum(&dstruct);
00067 
00068     }
00069     else
00070         /* else use ellipsoid defined in PROJ_INFO */
00071         str = G_find_key_value("ellps", proj_keys);
00072 
00073     if (str != NULL) {
00074         if (GPJ_get_ellipsoid_by_name(str, &estruct) < 0) {
00075             G_fatal_error(_("Invalid ellipsoid <%s> in file"), str);
00076         }
00077         else {
00078             *a = estruct.a;
00079             *e2 = estruct.es;
00080             *rf = estruct.rf;
00081             GPJ_free_ellps(&estruct);
00082             return 1;
00083         }
00084     }
00085     else {
00086         str3 = G_find_key_value("a", proj_keys);
00087         if (str3 != NULL) {
00088             G_asprintf(&str, "a=%s", str3);
00089             if ((str3 = G_find_key_value("es", proj_keys)) != NULL)
00090                 G_asprintf(&str1, "e=%s", str3);
00091             else if ((str3 = G_find_key_value("f", proj_keys)) != NULL)
00092                 G_asprintf(&str1, "f=1/%s", str3);
00093             else if ((str3 = G_find_key_value("rf", proj_keys)) != NULL)
00094                 G_asprintf(&str1, "f=1/%s", str3);
00095             else if ((str3 = G_find_key_value("b", proj_keys)) != NULL)
00096                 G_asprintf(&str1, "b=%s", str3);
00097             else
00098                 G_fatal_error(_("No secondary ellipsoid descriptor "
00099                                 "(rf, es or b) in file"));
00100 
00101             if (get_a_e2_rf(str, str1, a, e2, rf) == 0)
00102                 G_fatal_error(_("Invalid ellipsoid descriptors "
00103                                 "(a, rf, es or b) in file"));
00104             return 1;
00105         }
00106         else {
00107             str = G_find_key_value("proj", proj_keys);
00108             if ((str == NULL) || (strcmp(str, "ll") == 0)) {
00109                 *a = 6378137.0;
00110                 *e2 = .006694385;
00111                 *rf = 298.257223563;
00112                 return 0;
00113             }
00114             else {
00115                 G_fatal_error(_("No ellipsoid info given in file"));
00116             }
00117         }
00118     }
00119     return 1;
00120 }
00121 
00122 
00131 int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
00132 {
00133     struct ellps_list *list, *listhead;
00134 
00135     list = listhead = read_ellipsoid_table(0);
00136 
00137     while (list != NULL) {
00138         if (G_strcasecmp(name, list->name) == 0) {
00139             estruct->name = G_store(list->name);
00140             estruct->longname = G_store(list->longname);
00141             estruct->a = list->a;
00142             estruct->es = list->es;
00143             estruct->rf = list->rf;
00144             free_ellps_list(listhead);
00145             return 1;
00146         }
00147         list = list->next;
00148     }
00149     free_ellps_list(listhead);
00150     return -1;
00151 }
00152 
00153 static int
00154 get_a_e2_rf(const char *s1, const char *s2, double *a, double *e2,
00155             double *recipf)
00156 {
00157     double b, f;
00158 
00159     if (sscanf(s1, "a=%lf", a) != 1)
00160         return 0;
00161 
00162     if (*a <= 0.0)
00163         return 0;
00164 
00165     if (sscanf(s2, "e=%lf", e2) == 1) {
00166         f = 1.0 - sqrt(1.0 - *e2);
00167         *recipf = 1.0 / f;
00168         return (*e2 >= 0.0);
00169     }
00170 
00171     if (sscanf(s2, "f=1/%lf", recipf) == 1) {
00172         if (*recipf <= 0.0)
00173             return 0;
00174         f = 1.0 / *recipf;
00175         *e2 = f * (2 - f);
00176         return (*e2 >= 0.0);
00177     }
00178 
00179     if (sscanf(s2, "b=%lf", &b) == 1) {
00180         if (b <= 0.0)
00181             return 0;
00182         if (b == *a) {
00183             f = 0.0;
00184             *e2 = 0.0;
00185         }
00186         else {
00187             f = (*a - b) / *a;
00188             *e2 = f * (2 - f);
00189         }
00190         *recipf = 1.0 / f;
00191         return (*e2 >= 0.0);
00192     }
00193     return 0;
00194 }
00195 
00196 struct ellps_list *read_ellipsoid_table(int fatal)
00197 {
00198     FILE *fd;
00199     char file[GPATH_MAX];
00200     char buf[4096];
00201     char name[100], descr[1024], buf1[1024], buf2[1024];
00202     char badlines[1024];
00203     int line;
00204     int err;
00205     struct ellps_list *current = NULL, *outputlist = NULL;
00206     double a, e2, rf;
00207 
00208     int count = 0;
00209 
00210     sprintf(file, "%s%s", G_gisbase(), ELLIPSOIDTABLE);
00211     fd = fopen(file, "r");
00212 
00213     if (!fd) {
00214         (fatal ? G_fatal_error : G_warning)(
00215             _("Unable to open ellipsoid table file <%s>"), file);
00216         return NULL;
00217     }
00218 
00219     err = 0;
00220     *badlines = 0;
00221     for (line = 1; G_getl2(buf, sizeof buf, fd); line++) {
00222         G_strip(buf);
00223         if (*buf == 0 || *buf == '#')
00224             continue;
00225 
00226         if (sscanf(buf, "%s  \"%1023[^\"]\" %s %s", name, descr, buf1, buf2)
00227             != 4) {
00228             err++;
00229             sprintf(buf, " %d", line);
00230             if (*badlines)
00231                 G_strcat(badlines, ",");
00232             G_strcat(badlines, buf);
00233             continue;
00234         }
00235 
00236 
00237         if (get_a_e2_rf(buf1, buf2, &a, &e2, &rf)
00238             || get_a_e2_rf(buf2, buf1, &a, &e2, &rf)) {
00239             if (current == NULL)
00240                 current = outputlist = G_malloc(sizeof(struct ellps_list));
00241             else
00242                 current = current->next = G_malloc(sizeof(struct ellps_list));
00243             current->name = G_store(name);
00244             current->longname = G_store(descr);
00245             current->a = a;
00246             current->es = e2;
00247             current->rf = rf;
00248             current->next = NULL;
00249             count++;
00250         }
00251         else {
00252             err++;
00253             sprintf(buf, " %d", line);
00254             if (*badlines)
00255                 G_strcat(badlines, ",");
00256             G_strcat(badlines, buf);
00257             continue;
00258         }
00259     }
00260 
00261     fclose(fd);
00262 
00263     if (!err)
00264         return outputlist;
00265 
00266     (fatal ? G_fatal_error : G_warning)(
00267         err == 1
00268         ? _("Line%s of ellipsoid table file <%s> is invalid")
00269         : _("Lines%s of ellipsoid table file <%s> are invalid"),
00270         badlines, file);
00271 
00272     return outputlist;
00273 }
00274 
00275 void GPJ_free_ellps(struct gpj_ellps *estruct)
00276 {
00277     G_free(estruct->name);
00278     G_free(estruct->longname);
00279     return;
00280 }
00281 
00282 void free_ellps_list(struct ellps_list *elist)
00283 {
00284     struct ellps_list *old;
00285 
00286     while (elist != NULL) {
00287         G_free(elist->name);
00288         G_free(elist->longname);
00289         old = elist;
00290         elist = old->next;
00291         G_free(old);
00292     }
00293 
00294     return;
00295 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines