GRASS Programmer's Manual
6.4.1(2011)
|
00001 00023 #include <unistd.h> 00024 #include <ctype.h> 00025 #include <string.h> 00026 #include <stdlib.h> 00027 #include <math.h> /* for sqrt() */ 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 /* static int get_a_e2 (char *, char *, double *,double *); */ 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 * new 05/2000 by al: for datum shift the f parameter is needed too. 00146 * this all is not a clean design, but it keeps backward- 00147 * compatibility. 00148 * looks up ellipsoid in ellipsoid table and returns the 00149 * a, e2 and f parameters for the ellipsoid 00150 * 00151 * returns 1 if ok, 00152 * 0 if not found in table 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 /* return strcmp(a->name,b->name); */ 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 /* over correct typed version */ 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 }