GRASS Programmer's Manual  6.4.2(2012)
get_proj.c
Go to the documentation of this file.
00001 
00017 #include <stdio.h>
00018 #include <stdlib.h>
00019 #include <ctype.h>
00020 #include <math.h>
00021 #include <string.h>
00022 #include <grass/gis.h>
00023 #include <grass/gprojects.h>
00024 #include <grass/glocale.h>
00025 
00026 #define MAIN
00027 
00028 /* Finder function for datum conversion lookup tables */
00029 #define FINDERFUNC set_proj_lib
00030 #define PERMANENT "PERMANENT"
00031 #define MAX_PARGS 100
00032 
00033 static void alloc_options(char *);
00034 
00035 static char *opt_in[MAX_PARGS];
00036 static int nopt1;
00037 
00061 int pj_get_kv(struct pj_info *info, struct Key_Value *in_proj_keys,
00062               struct Key_Value *in_units_keys)
00063 {
00064     char *str;
00065     int i;
00066     double a, es, rf;
00067     int returnval = 1;
00068     char buffa[300], factbuff[50];
00069     char proj_in[50], *datum, *params;
00070     projPJ *pj;
00071 
00072     proj_in[0] = '\0';
00073     info->zone = 0;
00074     info->meters = 1.0;
00075     info->proj[0] = '\0';
00076 
00077     str = G_find_key_value("meters", in_units_keys);
00078     if (str != NULL) {
00079         strcpy(factbuff, str);
00080         if (strlen(factbuff) > 0)
00081             sscanf(factbuff, "%lf", &(info->meters));
00082     }
00083     str = G_find_key_value("name", in_proj_keys);
00084     if (str != NULL) {
00085         sprintf(proj_in, "%s", str);
00086     }
00087     str = G_find_key_value("proj", in_proj_keys);
00088     if (str != NULL) {
00089         sprintf(info->proj, "%s", str);
00090     }
00091     if (strlen(info->proj) <= 0)
00092         sprintf(info->proj, "ll");
00093 
00094     nopt1 = 0;
00095     for (i = 0; i < in_proj_keys->nitems; i++) {
00096         /* the name parameter is just for grasses use */
00097         if (strcmp(in_proj_keys->key[i], "name") == 0) {
00098             continue;
00099 
00100             /* zone handled separately at end of loop */
00101         }
00102         else if (strcmp(in_proj_keys->key[i], "zone") == 0) {
00103             continue;
00104 
00105             /* Datum and ellipsoid-related parameters will be handled 
00106              * separately after end of this loop PK */
00107 
00108         }
00109         else if (strcmp(in_proj_keys->key[i], "datum") == 0
00110                  || strcmp(in_proj_keys->key[i], "dx") == 0
00111                  || strcmp(in_proj_keys->key[i], "dy") == 0
00112                  || strcmp(in_proj_keys->key[i], "dz") == 0
00113                  || strcmp(in_proj_keys->key[i], "datumparams") == 0
00114                  || strcmp(in_proj_keys->key[i], "nadgrids") == 0
00115                  || strcmp(in_proj_keys->key[i], "towgs84") == 0
00116                  || strcmp(in_proj_keys->key[i], "ellps") == 0
00117                  || strcmp(in_proj_keys->key[i], "a") == 0
00118                  || strcmp(in_proj_keys->key[i], "b") == 0
00119                  || strcmp(in_proj_keys->key[i], "es") == 0
00120                  || strcmp(in_proj_keys->key[i], "f") == 0
00121                  || strcmp(in_proj_keys->key[i], "rf") == 0) {
00122             continue;
00123 
00124             /* PROJ.4 uses longlat instead of ll as 'projection name' */
00125 
00126         }
00127         else if (strcmp(in_proj_keys->key[i], "proj") == 0) {
00128             if (strcmp(in_proj_keys->value[i], "ll") == 0)
00129                 sprintf(buffa, "proj=longlat");
00130             else
00131                 sprintf(buffa, "proj=%s", in_proj_keys->value[i]);
00132 
00133             /* 'One-sided' PROJ.4 flags will have the value in
00134              * the key-value pair set to 'defined' and only the
00135              * key needs to be passed on. */
00136         }
00137         else if (strcmp(in_proj_keys->value[i], "defined") == 0)
00138             sprintf(buffa, in_proj_keys->key[i]);
00139 
00140         else
00141             sprintf(buffa, "%s=%s",
00142                     in_proj_keys->key[i], in_proj_keys->value[i]);
00143 
00144         alloc_options(buffa);
00145     }
00146 
00147     str = G_find_key_value("zone", in_proj_keys);
00148     if (str != NULL) {
00149         if (sscanf(str, "%d", &(info->zone)) != 1) {
00150             G_fatal_error(_("Invalid zone %s specified"), str);
00151         }
00152         if (info->zone < 0) {
00153 
00154             /* if zone is negative, write abs(zone) and define south */
00155             info->zone = -info->zone;
00156 
00157             if (G_find_key_value("south", in_proj_keys) == NULL) {
00158                 sprintf(buffa, "south");
00159                 alloc_options(buffa);
00160             }
00161         }
00162         sprintf(buffa, "zone=%d", info->zone);
00163         alloc_options(buffa);
00164     }
00165 
00166     if ((GPJ__get_ellipsoid_params(in_proj_keys, &a, &es, &rf) == 0)
00167         && (str = G_find_key_value("ellps", in_proj_keys)) != NULL) {
00168         /* Default values were returned but an ellipsoid name not recognised
00169          * by GRASS is present---perhaps it will be recognised by
00170          * PROJ.4 even though it wasn't by GRASS */
00171         sprintf(buffa, "ellps=%s", str);
00172         alloc_options(buffa);
00173     }
00174     else {
00175         sprintf(buffa, "a=%.16g", a);
00176         alloc_options(buffa);
00177         /* Cannot use es directly because the OSRImportFromProj4() 
00178          * function in OGR only accepts b or rf as the 2nd parameter */
00179         if (es == 0)
00180             sprintf(buffa, "b=%.16g", a);
00181         else
00182             sprintf(buffa, "rf=%.16g", rf);
00183         alloc_options(buffa);
00184 
00185     }
00186     /* Workaround to stop PROJ reading values from defaults file when
00187      * rf (and sometimes ellps) is not specified */
00188     if (G_find_key_value("no_defs", in_proj_keys) == NULL) {
00189         sprintf(buffa, "no_defs");
00190         alloc_options(buffa);
00191     }
00192 
00193     /* If datum parameters are present in the PROJ_INFO keys, pass them on */
00194     if (GPJ__get_datum_params(in_proj_keys, &datum, &params) == 2) {
00195         sprintf(buffa, params);
00196         alloc_options(buffa);
00197         G_free(params);
00198 
00199         /* else if a datum name is present take it and look up the parameters 
00200          * from the datum.table file */
00201     }
00202     else if (datum != NULL) {
00203 
00204         if (GPJ_get_default_datum_params_by_name(datum, &params) > 0) {
00205             sprintf(buffa, params);
00206             alloc_options(buffa);
00207             returnval = 2;
00208             G_free(params);
00209 
00210             /* else just pass the datum name on and hope it is recognised by 
00211              * PROJ.4 even though it isn't recognised by GRASS */
00212         }
00213         else {
00214             sprintf(buffa, "datum=%s", datum);
00215             alloc_options(buffa);
00216             returnval = 3;
00217         }
00218         G_free(datum);
00219         /* else there'll be no datum transformation taking place here... */
00220     }
00221     else {
00222         returnval = 4;
00223     }
00224 
00225     /* Set finder function for locating datum conversion tables PK */
00226     pj_set_finder(FINDERFUNC);
00227 
00228     if (!(pj = pj_init(nopt1, opt_in))) {
00229         strcpy(buffa,
00230                _("Unable to initialise PROJ.4 with the following parameter list:"));
00231         for (i = 0; i < nopt1; i++) {
00232             char err[50];
00233 
00234             sprintf(err, " +%s", opt_in[i]);
00235             strcat(buffa, err);
00236         }
00237         G_warning(buffa);
00238         G_warning(_("The error message: %s"), pj_strerrno(pj_errno));
00239         return -1;
00240     }
00241     info->pj = pj;
00242 
00243     return returnval;
00244 }
00245 
00246 static void alloc_options(char *buffa)
00247 {
00248     int nsize;
00249 
00250     nsize = strlen(buffa);
00251     opt_in[nopt1++] = (char *)G_malloc(nsize + 1);
00252     sprintf(opt_in[nopt1 - 1], buffa);
00253     return;
00254 }
00255 
00256 int pj_get_string(struct pj_info *info, char *str)
00257 {
00258     char *opt_in[MAX_PARGS];
00259     char *s;
00260     int nopt = 0;
00261     int nsize;
00262     char zonebuff[50], buffa[300];
00263     projPJ *pj;
00264 
00265     info->zone = 0;
00266     info->proj[0] = '\0';
00267     info->meters = 1.0;
00268 
00269     if ((str == NULL) || (str[0] == '\0')) {
00270         /* Null Pointer or empty string is supplied for parameters, 
00271          * implying latlong projection; just need to set proj 
00272          * parameter and call pj_init PK */
00273         sprintf(info->proj, "ll");
00274         sprintf(buffa, "proj=latlong ellps=WGS84");
00275         nsize = strlen(buffa);
00276         opt_in[nopt] = (char *)G_malloc(nsize + 1);
00277         sprintf(opt_in[nopt++], buffa);
00278     }
00279     else {
00280         /* Parameters have been provided; parse through them but don't
00281          * bother with most of the checks in pj_get_kv; assume the
00282          * programmer knows what he / she is doing when using this 
00283          * function rather than reading a PROJ_INFO file       PK */
00284         s = str;
00285         while (s = strtok(s, " \t\n"), s) {
00286             if (strncmp(s, "+unfact=", 8) == 0) {
00287                 s = s + 8;
00288                 info->meters = atof(s);
00289             }
00290             else {
00291                 if (strncmp(s, "+", 1) == 0)
00292                     ++s;
00293                 if (nsize = strlen(s), nsize) {
00294                     if (nopt >= MAX_PARGS) {
00295                         fprintf(stderr, "nopt = %d, s=%s\n", nopt, str);
00296                         G_fatal_error(_("Option input overflowed option table"));
00297                     }
00298 
00299                     if (strncmp("zone=", s, 5) == 0) {
00300                         sprintf(zonebuff, "%s", s + 5);
00301                         sscanf(zonebuff, "%d", &(info->zone));
00302                     }
00303 
00304                     if (strncmp("proj=", s, 5) == 0) {
00305                         sprintf(info->proj, "%s", s + 5);
00306                         if (strcmp(info->proj, "ll") == 0)
00307                             sprintf(buffa, "proj=latlong");
00308                         else
00309                             sprintf(buffa, s);
00310                     }
00311                     else {
00312                         sprintf(buffa, s);
00313                     }
00314                     nsize = strlen(buffa);
00315                     opt_in[nopt] = (char *)G_malloc(nsize + 1);
00316                     sprintf(opt_in[nopt++], buffa);
00317                 }
00318             }
00319             s = 0;
00320         }
00321     }
00322 
00323     /* Set finder function for locating datum conversion tables PK */
00324     pj_set_finder(FINDERFUNC);
00325 
00326     if (!(pj = pj_init(nopt, opt_in))) {
00327         G_warning(_("Unable to initialize pj cause: %s"),
00328                   pj_strerrno(pj_errno));
00329         return -1;
00330     }
00331     info->pj = pj;
00332 
00333     return 1;
00334 }
00335 
00352 int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
00353 {
00354     pjnew->meters = 1.;
00355     pjnew->zone = 0;
00356     sprintf(pjnew->proj, "ll");
00357     if ((pjnew->pj = pj_latlong_from_proj(pjold->pj)) == NULL)
00358         return -1;
00359     else
00360         return 1;
00361 }
00362 
00363 /* set_proj_lib()
00364  * 'finder function' for use with PROJ.4 pj_set_finder() function */
00365 
00366 const char *set_proj_lib(const char *name)
00367 {
00368     const char *gisbase = G_gisbase();
00369     static char *buf = NULL;
00370     static size_t buf_len;
00371     size_t len = strlen(gisbase) + sizeof(GRIDDIR) + strlen(name) + 1;
00372 
00373     if (buf_len < len) {
00374         if (buf != NULL)
00375             G_free(buf);
00376         buf_len = len + 20;
00377         buf = G_malloc(buf_len);
00378     }
00379 
00380     sprintf(buf, "%s%s/%s", gisbase, GRIDDIR, name);
00381 
00382     return buf;
00383 }
00384 
00396 int pj_print_proj_params(struct pj_info *iproj, struct pj_info *oproj)
00397 {
00398     char *str;
00399 
00400     if (iproj) {
00401         str = pj_get_def(iproj->pj, 1);
00402         if (str != NULL) {
00403             fprintf(stderr, "%s: %s\n", _("Input Projection Parameters"),
00404                     str);
00405             pj_dalloc(str);
00406             fprintf(stderr, "%s: %.16g\n", _("Input Unit Factor"),
00407                     iproj->meters);
00408         }
00409         else
00410             return -1;
00411     }
00412 
00413     if (oproj) {
00414         str = pj_get_def(oproj->pj, 1);
00415         if (str != NULL) {
00416             fprintf(stderr, "%s: %s\n", _("Output Projection Parameters"),
00417                     str);
00418             pj_dalloc(str);
00419             fprintf(stderr, "%s: %.16g\n", _("Output Unit Factor"),
00420                     oproj->meters);
00421         }
00422         else
00423             return -1;
00424     }
00425 
00426     return 1;
00427 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines