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