GRASS Programmer's Manual
6.4.2(2012)
|
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 }