GRASS Programmer's Manual  6.4.2(2012)
proj/datum.c
Go to the documentation of this file.
00001 
00016 #include <unistd.h>
00017 #include <string.h>
00018 #include <ctype.h>
00019 #include <stdlib.h>
00020 
00021 #include <grass/gis.h>
00022 #include <grass/glocale.h>
00023 #include <grass/gprojects.h>
00024 #include "local_proto.h"
00025 
00037 int GPJ_get_datum_by_name(const char *name, struct gpj_datum *dstruct)
00038 {
00039     struct datum_list *list, *listhead;
00040 
00041     list = listhead = read_datum_table();
00042 
00043     while (list != NULL) {
00044         if (G_strcasecmp(name, list->name) == 0) {
00045             dstruct->name = G_store(list->name);
00046             dstruct->longname = G_store(list->longname);
00047             dstruct->ellps = G_store(list->ellps);
00048             dstruct->dx = list->dx;
00049             dstruct->dy = list->dy;
00050             dstruct->dz = list->dz;
00051             free_datum_list(listhead);
00052             return 1;
00053         }
00054         list = list->next;
00055     }
00056     free_datum_list(listhead);
00057     return -1;
00058 }
00059 
00085 int GPJ_get_default_datum_params_by_name(const char *name, char **params)
00086 {
00087     struct gpj_datum_transform_list *list, *old;
00088     int count = 1;
00089 
00090     list = GPJ_get_datum_transform_by_name(name);
00091 
00092     if (list == NULL) {
00093         *params = NULL;
00094         return -1;
00095     }
00096 
00097     /* Take the first parameter set in the list as the default
00098      * (will normally be a 3-parameter transformation)        */
00099     *params = G_store(list->params);
00100 
00101     while (list->next != NULL) {
00102         count++;
00103         old = list;
00104         list = list->next;
00105         G_free(old);
00106     }
00107 
00108     G_free(list);
00109     return count;
00110 
00111 }
00112 
00136 int GPJ_get_datum_params(char **name, char **params)
00137 {
00138     int ret;
00139     struct Key_Value *proj_keys = G_get_projinfo();
00140 
00141     ret = GPJ__get_datum_params(proj_keys, name, params);
00142     G_free_key_value(proj_keys);
00143 
00144     return ret;
00145 }
00146 
00174 int GPJ__get_datum_params(struct Key_Value *projinfo,
00175                           char **datumname, char **params)
00176 {
00177     int returnval = -1;
00178 
00179     if (NULL != G_find_key_value("datum", projinfo)) {
00180         *datumname = G_store(G_find_key_value("datum", projinfo));
00181         returnval = 1;
00182     }
00183     else
00184         *datumname = NULL;
00185 
00186     if (G_find_key_value("datumparams", projinfo) != NULL) {
00187         *params = G_store(G_find_key_value("datumparams", projinfo));
00188         returnval = 2;
00189     }
00190     else if (G_find_key_value("nadgrids", projinfo) != NULL) {
00191         const char *gisbase = G_gisbase();
00192 
00193         G_asprintf(params, "nadgrids=%s%s/%s", gisbase, GRIDDIR,
00194                    G_find_key_value("nadgrids", projinfo));
00195         returnval = 2;
00196     }
00197     else if (G_find_key_value("towgs84", projinfo) != NULL) {
00198         G_asprintf(params, "towgs84=%s",
00199                    G_find_key_value("towgs84", projinfo));
00200         returnval = 2;
00201     }
00202     else if (G_find_key_value("dx", projinfo) != NULL
00203              && G_find_key_value("dy", projinfo) != NULL
00204              && G_find_key_value("dz", projinfo) != NULL) {
00205         G_asprintf(params, "towgs84=%s,%s,%s",
00206                    G_find_key_value("dx", projinfo),
00207                    G_find_key_value("dy", projinfo),
00208                    G_find_key_value("dz", projinfo));
00209         returnval = 2;
00210     }
00211     else
00212         *params = NULL;
00213 
00214     return returnval;
00215 
00216 }
00217 
00240 int GPJ_ask_datum_params(const char *datumname, char **params)
00241 {
00242     char buff[1024], answer[100];
00243     char *Tmp_file;
00244     FILE *Tmp_fd = NULL;
00245     struct gpj_datum_transform_list *list, *listhead, *old;
00246     int transformcount, currenttransform;
00247 
00248     if (G_strcasecmp(datumname, "custom") != 0) {
00249         Tmp_file = G_tempfile();
00250         if (NULL == (Tmp_fd = fopen(Tmp_file, "w"))) {
00251             G_warning(_("Unable to open temporary file"));
00252         }
00253 
00254         fprintf(Tmp_fd, "Number\tDetails\t\n---\n");
00255         listhead = GPJ_get_datum_transform_by_name(datumname);
00256         list = listhead;
00257         transformcount = 0;
00258         while (list != NULL) {
00259             /* Count how many sets of transformation paramters have been 
00260              * defined for this datum and print them to a temporary file 
00261              * in case the user asks for them to be displayed */
00262             fprintf(Tmp_fd,
00263                     "%d\tUsed in %s\n\t(PROJ.4 Params %s)\n\t%s\n---\n",
00264                     list->count, list->where_used, list->params,
00265                     list->comment);
00266             list = list->next;
00267             transformcount++;
00268         }
00269         fclose(Tmp_fd);
00270 
00271         for (;;) {
00272             do {
00273                 fprintf(stderr,
00274                         ("\nNow select Datum Transformation Parameters\n"));
00275                 fprintf(stderr,
00276                         ("Please think carefully about the area covered by your data\n"
00277                          "and the accuracy you require before making your selection.\n"));
00278                 fprintf(stderr,
00279                         ("\nEnter 'list' to see the list of available Parameter sets\n"));
00280                 fprintf(stderr,
00281                         ("Enter the corresponding number, or <RETURN> to cancel request\n"));
00282                 fprintf(stderr, ">");
00283             } while (!G_gets(answer));
00284             G_strip(answer);
00285             if (strlen(answer) == 0) {
00286                 remove(Tmp_file);
00287                 G_free(Tmp_file);
00288                 return -1;
00289             }
00290             if (strcmp(answer, "list") == 0) {
00291                 char *pager;
00292 
00293                 pager = getenv("GRASS_PAGER");
00294                 if (!pager || strlen(pager) == 0)
00295                     pager = "cat";
00296 
00297                 /* Always print interactive output to stderr */
00298                 sprintf(buff, "%s \"%s\" 1>&2", pager,
00299                         G_convert_dirseps_to_host(Tmp_file));
00300                 G_system(buff);
00301             }
00302             else {
00303                 if ((sscanf(answer, "%d", &currenttransform) != 1) ||
00304                     currenttransform > transformcount ||
00305                     currenttransform < 1) {
00306 
00307                     /* If a number was not typed, was less than 0 or greater
00308                      * than the number of sets of parameters, ask again */
00309                     fprintf(stderr, ("\ninvalid transformation number\n"));
00310                 }
00311                 else
00312                     break;
00313             }
00314 
00315         }
00316         remove(Tmp_file);
00317         G_free(Tmp_file);
00318 
00319         list = listhead;
00320         while (list != NULL) {
00321             /* Search through the linked list to find the parameter string
00322              * that corresponds to the number entered */
00323             if (list->count == currenttransform)
00324                 G_asprintf(params, list->params);
00325 
00326             /* Continue to end of list even after we find it, to free all
00327              * the memory used */
00328             old = list;
00329             list = old->next;
00330             G_free(old);
00331         }
00332     }
00333     else {
00334         /* Here we ask the user to enter customised parameters */
00335         for (;;) {
00336             do {
00337                 fprintf(stderr,
00338                         ("\nPlease specify datum transformation parameters in PROJ.4 syntax. Examples:\n"));
00339                 fprintf(stderr,
00340                         ("\ttowgs84=dx,dy,dz\t(3-parameter transformation)\n"));
00341                 fprintf(stderr,
00342                         ("\ttowgs84=dx,dy,dz,rx,ry,rz,m\t(7-parameter transformation)\n"));
00343                 fprintf(stderr,
00344                         ("\tnadgrids=alaska\t(Tables-based grid-shifting transformation)\n"));
00345                 fprintf(stderr, _("Hit RETURN to cancel request\n"));
00346                 fprintf(stderr, ">");
00347             } while (!G_gets(answer));
00348             G_strip(answer);
00349             if (strlen(answer) == 0)
00350                 return -1;
00351             G_asprintf(params, answer);
00352             sprintf(buff,
00353                     "Parameters to be used are:\n\"%s\"\nIs this correct?",
00354                     *params);
00355             if (G_yes(buff, 1))
00356                 break;
00357 
00358         }
00359 
00360     }
00361 
00362     return 1;
00363 
00364 }
00365 
00378 struct gpj_datum_transform_list *GPJ_get_datum_transform_by_name(const char
00379                                                                  *inputname)
00380 {
00381     FILE *fd;
00382     char file[GPATH_MAX];
00383     char buf[1024];
00384     int line;
00385     struct gpj_datum_transform_list *current = NULL, *outputlist = NULL;
00386     struct gpj_datum dstruct;
00387     int count = 0;
00388 
00389     GPJ_get_datum_by_name(inputname, &dstruct);
00390     if (dstruct.dx < 99999 && dstruct.dy < 99999 && dstruct.dz < 99999) {
00391         /* Include the old-style dx dy dz parameters from datum.table at the 
00392          * start of the list, unless these have been set to all 99999 to 
00393          * indicate only entries in datumtransform.table should be used */
00394         if (current == NULL)
00395             current = outputlist =
00396                 G_malloc(sizeof(struct gpj_datum_transform_list));
00397         else
00398             current = current->next =
00399                 G_malloc(sizeof(struct gpj_datum_transform_list));
00400         G_asprintf(&(current->params), "towgs84=%.3f,%.3f,%.3f", dstruct.dx,
00401                    dstruct.dy, dstruct.dz);
00402         G_asprintf(&(current->where_used), "whole %s region", inputname);
00403         G_asprintf(&(current->comment),
00404                    "Default 3-Parameter Transformation (May not be optimum for "
00405                    "older datums; use this only if no more appropriate options "
00406                    "are available.)");
00407         count++;
00408         current->count = count;
00409         current->next = NULL;
00410     }
00411     GPJ_free_datum(&dstruct);
00412 
00413     /* Now check for additional parameters in datumtransform.table */
00414 
00415     sprintf(file, "%s%s", G_gisbase(), DATUMTRANSFORMTABLE);
00416 
00417     fd = fopen(file, "r");
00418     if (!fd) {
00419         G_warning(_("Unable to open datum table file <%s>"), file);
00420         return outputlist;
00421     }
00422 
00423     for (line = 1; G_getl2(buf, sizeof(buf), fd); line++) {
00424         char name[100], params[1024], where_used[1024], comment[1024];
00425 
00426         G_strip(buf);
00427         if (*buf == '\0' || *buf == '#')
00428             continue;
00429 
00430         if (sscanf(buf, "%99s \"%1023[^\"]\" \"%1023[^\"]\" \"%1023[^\"]\"",
00431                    name, params, where_used, comment) != 4) {
00432             G_warning(_("Error in datum table file <%s>, line %d"), file,
00433                       line);
00434             continue;
00435         }
00436 
00437         if (G_strcasecmp(inputname, name) == 0) {
00438             /* If the datum name in this line matches the one we are 
00439              * looking for, add an entry to the linked list */
00440             if (current == NULL)
00441                 current = outputlist =
00442                     G_malloc(sizeof(struct gpj_datum_transform_list));
00443             else
00444                 current = current->next =
00445                     G_malloc(sizeof(struct gpj_datum_transform_list));
00446             current->params = G_store(params);
00447             current->where_used = G_store(where_used);
00448             current->comment = G_store(comment);
00449             count++;
00450             current->count = count;
00451             current->next = NULL;
00452         }
00453     }
00454 
00455     fclose(fd);
00456 
00457     return outputlist;
00458 
00459 }
00460 
00471 struct datum_list *read_datum_table(void)
00472 {
00473     FILE *fd;
00474     char file[GPATH_MAX];
00475     char buf[4096];
00476     int line;
00477     struct datum_list *current = NULL, *outputlist = NULL;
00478     int count = 0;
00479 
00480     sprintf(file, "%s%s", G_gisbase(), DATUMTABLE);
00481 
00482     fd = fopen(file, "r");
00483     if (!fd) {
00484         G_warning(_("Unable to open datum table file <%s>"), file);
00485         return NULL;
00486     }
00487 
00488     for (line = 1; G_getl2(buf, sizeof(buf), fd); line++) {
00489         char name[100], descr[1024], ellps[100];
00490         double dx, dy, dz;
00491 
00492         G_strip(buf);
00493         if (*buf == '\0' || *buf == '#')
00494             continue;
00495 
00496         if (sscanf(buf, "%s \"%1023[^\"]\" %s dx=%lf dy=%lf dz=%lf",
00497                    name, descr, ellps, &dx, &dy, &dz) != 6) {
00498             G_warning(_("Error in datum table file <%s>, line %d"), file,
00499                       line);
00500             continue;
00501         }
00502 
00503         if (current == NULL)
00504             current = outputlist = G_malloc(sizeof(struct datum_list));
00505         else
00506             current = current->next = G_malloc(sizeof(struct datum_list));
00507         current->name = G_store(name);
00508         current->longname = G_store(descr);
00509         current->ellps = G_store(ellps);
00510         current->dx = dx;
00511         current->dy = dy;
00512         current->dz = dz;
00513         current->next = NULL;
00514 
00515         count++;
00516     }
00517 
00518     fclose(fd);
00519 
00520     return outputlist;
00521 }
00522 
00529 void GPJ_free_datum(struct gpj_datum *dstruct)
00530 {
00531     G_free(dstruct->name);
00532     G_free(dstruct->longname);
00533     G_free(dstruct->ellps);
00534     return;
00535 }
00536 
00543 void free_datum_list(struct datum_list *dstruct)
00544 {
00545     struct datum_list *old;
00546 
00547     while (dstruct != NULL) {
00548         G_free(dstruct->name);
00549         G_free(dstruct->longname);
00550         G_free(dstruct->ellps);
00551         old = dstruct;
00552         dstruct = old->next;
00553         G_free(old);
00554     }
00555 
00556     return;
00557 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines