datum.c

Go to the documentation of this file.
00001 /*
00002  * $Id: datum.c,v 2.0 2004/11/09 12:16:06 bernhard Exp $
00003  *
00004  ****************************************************************************
00005  *
00006  * MODULE:       gis library
00007  * AUTHOR(S):    Andreas Lange - andreas.lange@rhein-main.de
00008  *               Paul Kelly - paul-grass@stjohnspoint.co.uk
00009  * PURPOSE:      provide functions for reading datum parameters from the
00010  *               location database.     
00011  * COPYRIGHT:    (C) 2000, 2003 by the GRASS Development Team
00012  *
00013  *               This program is free software under the GNU General Public
00014  *               License (>=v2). Read the file COPYING that comes with GRASS
00015  *               for details.
00016  *
00017  *****************************************************************************/
00018 
00019 #define DATUMTABLE "/etc/datum.table"
00020 #define DATUMTRANSFORMTABLE "/etc/datumtransform.table"
00021 
00022 #include <unistd.h>
00023 #include <string.h>
00024 #include <ctype.h>
00025 #include <stdlib.h>
00026 
00027 #include "gis.h"
00028 #include "glocale.h"
00029 
00030 static struct table
00031 {
00032     char *name;   /* Short Name / acronym of map datum */
00033     char *descr;  /* Long Name for map datum */
00034     char *ellps;  /* acronym for ellipsoid used with this datum */
00035     double dx;    /* delta x */
00036     double dy;    /* delta y */
00037     double dz;    /* delta z */
00038 } *table;
00039 
00040 struct datum_transform_list
00041 {
00042     int count;        /* Transform Number (ordered list) */
00043     char *params;     /* PROJ.4-style datum transform parameters */
00044     char *where_used; /* Comment text describing where (geographically)
00045                          the transform is valid */
00046     char *comment;    /* Additional Comments */
00047     struct datum_transform_list *next;  /* Pointer to next set of 
00048                                       transform parameters in linked list */
00049 };
00050 
00051 static int size;
00052 static int count = -1;
00053 
00054 static int compare_table_names(const void *, const void *);
00055 static void read_datum_table(void);
00056 static struct datum_transform_list * 
00057                           get_datum_transform_by_name(const char *inputname);
00058 
00059 int 
00060 G_get_datum_by_name(const char *name)
00061 {
00062     int i;
00063 
00064     read_datum_table();
00065 
00066     for (i = 0; i < count; i++)
00067         if (G_strcasecmp(name, table[i].name) == 0)
00068             return i;
00069 
00070     return -1;
00071 }
00072 
00073 /* this sets the datum shift parameters for datum pointed to by name */
00074 int 
00075 G_datum_shift(int n, double *dx, double *dy, double *dz)
00076 {
00077     read_datum_table();
00078 
00079     if (n < 0 || n >= count) 
00080         return 0; 
00081 
00082     *dx = table[n].dx;
00083     *dy = table[n].dy;
00084     *dz = table[n].dz;
00085     return 1;
00086 }
00087 
00088 /* set the ellipsoid name and parameters for datum */
00089 int 
00090 G_datum_parameters(int n, char *ellps, double *dx, double *dy, double *dz)
00091 {
00092     read_datum_table();
00093 
00094     if (n < 0 || n >= count) 
00095         return 0; 
00096 
00097     G_strcpy(ellps, table[n].ellps);
00098     *dx = table[n].dx;
00099     *dy = table[n].dy;
00100     *dz = table[n].dz;
00101     return 1;
00102 }
00103 
00104 char *
00105 G_datum_name(int n)
00106 {
00107     read_datum_table();
00108 
00109     if (n < 0 || n >= count) 
00110         return NULL; 
00111 
00112     return table[n].name;
00113 }
00114 
00115 char *
00116 G_datum_description(int n)
00117 { 
00118     read_datum_table();
00119   
00120     if (n < 0 || n >= count)
00121         return NULL;
00122 
00123     return table[n].descr;
00124 }
00125 
00126 char *
00127 G_datum_ellipsoid(int n)
00128 {
00129     read_datum_table();
00130 
00131     if (n < 0 || n >= count)
00132         return NULL;
00133 
00134     return table[n].ellps;
00135 }
00136 
00137 /***********************************************************
00138  *  G_ask_datum_params(datumname, params)
00139  *     char *datumname   String containing datum name that
00140  *                       parameters are to be found for. Must
00141  *                       exist in datum.table or be "custom"
00142  *     char *params      Pointer into which a string containing
00143  *                       the datum parameters chosen by the user
00144  *                       will be placed.
00145  *
00146  *  Interactively ask for datum parameters for a particular datum
00147  *  
00148  *  returns: 1 ok, -1 error
00149  ************************************************************/
00150 
00151 int G_ask_datum_params(char *datumname, char *params)
00152 {
00153     char buff[1024], answer[100];
00154     char *Tmp_file;
00155     FILE  *Tmp_fd = NULL;
00156     struct datum_transform_list *list, *listhead, *old;
00157     int transformcount, currenttransform, i;
00158 
00159     if( G_strcasecmp(datumname, "custom") != 0)
00160     {
00161         Tmp_file = G_tempfile ();
00162         if (NULL == (Tmp_fd = fopen (Tmp_file, "w"))) {
00163             G_warning(_("Cannot open temp file") );
00164         }
00165 
00166         fprintf(Tmp_fd,"Number\tDetails\t\n---\n");
00167         listhead = get_datum_transform_by_name( datumname );
00168         list = listhead;
00169         transformcount = 0;
00170         while( list != NULL)
00171         {          
00172             /* Count how many sets of transformation paramters have been 
00173              * defined for this datum and print them to a temporary file 
00174              * in case the user asks for them to be displayed */
00175             fprintf(Tmp_fd,"%d\tUsed in %s\n\t(PROJ.4 Params %s)\n\t%s\n---\n",
00176                    list->count, list->where_used, list->params, list->comment);
00177             list = list->next;
00178             transformcount++;
00179         }      
00180         fclose(Tmp_fd);
00181 
00182         for(;;) {
00183             do {
00184                 fprintf(stderr,("\nNow select Datum Transformation Parameters\n"));
00185                 fprintf(stderr,("Enter 'list' to see the list of available Parameter sets\n"));
00186                 fprintf (stderr,("Enter the corresponding number, or <RETURN> to cancel request\n"));
00187                 fprintf(stderr,">");
00188             } while(!G_gets(answer));
00189             G_strip(answer); 
00190             if(strlen(answer)==0) 
00191             {
00192                 remove( Tmp_file );
00193                 G_free( Tmp_file );
00194                 return -1;
00195             }
00196             if (strcmp(answer,"list") == 0) {
00197                 if (isatty(1))
00198                     sprintf(buff,"$GRASS_PAGER %s",Tmp_file);
00199                 else
00200                     sprintf(buff,"cat %s",Tmp_file);
00201                 G_system(buff);
00202             }
00203             else {
00204                 if ( (sscanf(answer, "%d", &currenttransform) != 1) ||
00205                     currenttransform > transformcount || currenttransform < 1) {
00206 
00207                     /* If a number was not typed, was less than 0 or greater
00208                      * than the number of sets of parameters, ask again */
00209                     fprintf(stderr,("\ninvalid transformation number\n"));
00210                 }
00211                 else break;
00212             }
00213 
00214         }
00215         remove ( Tmp_file );
00216         G_free ( Tmp_file );
00217    
00218         list = listhead;
00219         while (list != NULL)
00220         {
00221             /* Search through the linked list to find the parameter string
00222              * that corresponds to the number entered */
00223             if( list->count == currenttransform )
00224                 sprintf(params, list->params);
00225            
00226             /* Continue to end of list even after we find it, to free all
00227              * the memory used */
00228             old = list;
00229             list = old->next;
00230             G_free( old );
00231         }
00232     }
00233     else
00234     {
00235         /* Here we ask the user to enter customised parameters */
00236         for(;;) {
00237             do {
00238                 fprintf(stderr,("\nPlease specify datum transformation parameters in PROJ.4 syntax. Examples:\n"));
00239                 fprintf(stderr,("\ttowgs84=dx,dy,dz\t(3-parameter transformation)\n"));
00240                 fprintf(stderr,("\ttowgs84=dx,dy,dz,rx,ry,rz,m\t(7-parameter transformation)\n"));
00241                 fprintf(stderr,("\tnadgrids=alaska\t(Tables-based grid-shifting transformation)\n"));
00242                 fprintf (stderr,_("Hit RETURN to cancel request\n"));
00243                 fprintf(stderr,">");
00244             } while(!G_gets(answer));
00245             G_strip(answer); 
00246             if(strlen(answer)==0)
00247                 return -1;
00248             sprintf(params, answer);
00249             sprintf(buff, "Parameters to be used are:\n\"%s\"\nIs this correct?", params);
00250             if (G_yes(buff, 1))
00251                 break;
00252 
00253         }
00254 
00255     }
00256    
00257     return 1;
00258 
00259 }
00260 
00261 /***********************************************************
00262  *  G_get_datumparams_from_projinfo(projinfo, datumname, params)
00263  *     struct Key_Value *projinfo Set of key_value pairs containing
00264  *                       projection information in PROJ_INFO file
00265  *                       format
00266  *     char *datumname   Pointer into which a string containing
00267  *                       the datum name (if present) will be
00268  *                       placed.
00269  *     char *params      Pointer into which a string containing
00270  *                       the datum parameters (if present) will
00271  *                       be placed.
00272  *
00273  *  Extract the datum transformation-related parameters from a 
00274  *  set of general PROJ_INFO parameters.
00275  *  This function can be used to test if a location set-up 
00276  *  supports datum transformation.
00277  *  
00278  *  returns: -1 error or no datum information found, 
00279  *           1 only datum name found, 2 params found
00280  ************************************************************/
00281 
00282 int G_get_datumparams_from_projinfo(struct Key_Value *projinfo, 
00283                                     char *datumname, char *params)
00284 {
00285     int returnval = -1;
00286    
00287     if( NULL != G_find_key_value("datum", projinfo) )
00288     {
00289         sprintf(datumname, G_find_key_value("datum", projinfo));
00290         returnval = 1;
00291     }
00292           
00293     if( G_find_key_value("datumparams", projinfo) != NULL )
00294     {
00295         sprintf(params, G_find_key_value("datumparams", projinfo));
00296         returnval = 2;
00297     }
00298     else if( G_find_key_value("nadgrids", projinfo) != NULL )
00299     {
00300         sprintf(params, "nadgrids=%s", G_find_key_value("nadgrids", projinfo));
00301         returnval = 2;
00302     }
00303     else if( G_find_key_value("towgs84", projinfo) != NULL )
00304     {
00305         sprintf(params, "towgs84=%s", G_find_key_value("towgs84", projinfo));
00306         returnval = 2;
00307     }
00308     else if( G_find_key_value("dx", projinfo) != NULL
00309           && G_find_key_value("dy", projinfo) != NULL
00310           && G_find_key_value("dz", projinfo) != NULL ) 
00311     {
00312         sprintf(params, "towgs84=%s,%s,%s",
00313                 G_find_key_value("dx", projinfo),
00314                 G_find_key_value("dy", projinfo),
00315                 G_find_key_value("dz", projinfo) );
00316         returnval = 2;
00317     }
00318 
00319     return returnval;
00320    
00321 }
00322 
00323 /*
00324  ***********************************************************
00325  *  get_datum_transform_by_name(inputname)
00326  *     char *inputname   String containing the datum name we
00327  *                       are going to look up parameters for
00328  * 
00329  *  Internal function to find all possible sets of 
00330  *  transformation parameters for a particular datum
00331  *  
00332  *  returns: Pointer to struct datum_transform_list (a linked
00333  *           list containing transformation parameters),
00334  *           or NULL if no suitable parameters were found.
00335  ************************************************************/
00336 
00337 static struct datum_transform_list * 
00338                           get_datum_transform_by_name(const char *inputname)
00339 {
00340     FILE *fd;
00341     char file[1024];
00342     char buf[1024];
00343     int line;
00344     struct datum_transform_list *current=NULL, *outputlist=NULL;
00345     double dx, dy, dz;
00346     int count = 0;
00347    
00348     sprintf(file, "%s%s", G_gisbase(), DATUMTRANSFORMTABLE);
00349 
00350     fd = fopen(file, "r");
00351     if (!fd)
00352     {
00353         G_warning(_("unable to open datum table file: %s"), file);
00354         return NULL;
00355     }
00356 
00357     for (line = 1; G_getl(buf, sizeof(buf), fd); line++)
00358     {
00359         char name[100], params[256], where_used[256], comment[256];
00360 
00361         G_strip(buf);
00362         if (*buf == '\0' || *buf == '#')
00363             continue;
00364 
00365         if (sscanf(buf, "%99s \"%255[^\"]\" \"%255[^\"]\" \"%255[^\"]\"",
00366                    name, params, where_used, comment) != 4)
00367         {
00368             G_warning(_("error in datum table file, line %d"), line);
00369             continue;
00370         }
00371 
00372         if ( G_strcasecmp(inputname, name) == 0 )
00373         {
00374             /* If the datum name in this line matches the one we are 
00375              * looking for, add an entry to the linked list */
00376             if(current == NULL)
00377                 current = outputlist = G_malloc( sizeof(struct datum_transform_list) );
00378             else
00379                 current = current->next = G_malloc( sizeof(struct datum_transform_list) );
00380             current->params = G_store(params);
00381             current->where_used = G_store(where_used);
00382             current->comment = G_store(comment);           
00383             count++;
00384             current->count = count;
00385             current->next = NULL;
00386         }           
00387     } 
00388    
00389     G_datum_shift( G_get_datum_by_name( inputname ), &dx, &dy, &dz);
00390     if( dx < 99999 && dy < 99999 && dz < 99999 )
00391     {
00392         /* Include the old-style dx dy dz parameters from datum.table at the 
00393          * end of the list, unless these have been set to all 99999 to 
00394          * indicate only entries in datumtransform.table should be used */
00395         if(current == NULL)
00396             current = outputlist = G_malloc( sizeof(struct datum_transform_list) );
00397         else
00398             current = current->next = G_malloc( sizeof(struct datum_transform_list) );
00399         sprintf(buf, "towgs84=%.3f,%.3f,%.3f", dx, dy, dz);
00400         current->params = G_store(buf);
00401         sprintf(buf, "Default %s region", inputname);
00402         current->where_used = G_store(buf);
00403         sprintf(buf, "Default 3-Parameter Transformation");
00404         current->comment = G_store(buf);
00405         count++;
00406         current->count = count;
00407         current->next = NULL;
00408     }   
00409    
00410     
00411     return outputlist;
00412 
00413 }
00414 
00415 static void
00416 read_datum_table(void) 
00417 {
00418     FILE *fd;
00419     char file[1024];
00420     char buf[1024];
00421     int line;
00422 
00423     if (count >= 0)
00424         return;
00425 
00426     count = 0;
00427 
00428     sprintf(file, "%s%s", G_gisbase(), DATUMTABLE);
00429 
00430     fd = fopen(file, "r");
00431     if (!fd)
00432     {
00433         G_warning(_("unable to open datum table file: %s"), file);
00434         return;
00435     }
00436 
00437     for (line = 1; G_getl(buf, sizeof(buf), fd); line++)
00438     {
00439         char name[100], descr[100], ellps[100];
00440         double dx, dy, dz;
00441         struct table *t;
00442 
00443         G_strip(buf);
00444         if (*buf == '\0' || *buf == '#')
00445             continue;
00446 
00447         if (count >= size)
00448         {
00449             size += 50;
00450             table = G_realloc(table, size * sizeof(struct table));
00451         }
00452 
00453         t = &table[count];
00454 
00455         if (sscanf(buf, "%s \"%99[^\"]\" %s dx=%lf dy=%lf dz=%lf",
00456                    name, descr, ellps, &t->dx, &t->dy, &t->dz) != 6)
00457         {
00458             G_warning(_("error in datum table file, line %d"), line);
00459             continue;
00460         }
00461 
00462         t->name  = G_store (name);
00463         t->descr = G_store (descr);
00464         t->ellps = G_store (ellps);
00465 
00466         count++;
00467     }
00468  
00469     qsort(table, count, sizeof(struct table), compare_table_names);
00470 }
00471 
00472 static int
00473 compare_table_names(const void *aa, const void *bb)
00474 {
00475     const struct table *a = aa;
00476     const struct table *b = bb;
00477 
00478     return G_strcasecmp(a->name, b->name);
00479 }

Generated on Wed Aug 23 17:49:22 2006 for GRASS by  doxygen 1.4.7