/Proj4/pj_init.c
C | 434 lines | 244 code | 55 blank | 135 comment | 113 complexity | 0b27f5927aa5290611f5c6665c0e0134 MD5 | raw file
1/****************************************************************************** 2 * $Id: pj_init.c,v 1.19 2007/11/26 00:21:59 fwarmerdam Exp $ 3 * 4 * Project: PROJ.4 5 * Purpose: Initialize projection object from string definition. Includes 6 * pj_init(), pj_init_plus() and pj_free() function. 7 * Author: Gerald Evenden, Frank Warmerdam <warmerdam@pobox.com> 8 * 9 ****************************************************************************** 10 * Copyright (c) 1995, Gerald Evenden 11 * Copyright (c) 2002, Frank Warmerdam <warmerdam@pobox.com> 12 * 13 * Permission is hereby granted, free of charge, to any person obtaining a 14 * copy of this software and associated documentation files (the "Software"), 15 * to deal in the Software without restriction, including without limitation 16 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 17 * and/or sell copies of the Software, and to permit persons to whom the 18 * Software is furnished to do so, subject to the following conditions: 19 * 20 * The above copyright notice and this permission notice shall be included 21 * in all copies or substantial portions of the Software. 22 * 23 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 24 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 25 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 26 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 27 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 28 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 29 * DEALINGS IN THE SOFTWARE. 30 ****************************************************************************** 31 * 32 * $Log: pj_init.c,v $ 33 * Revision 1.19 2007/11/26 00:21:59 fwarmerdam 34 * Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before 35 * adjustment for spherical projections. 36 * Modified pj_datum_transform() to use the original ellipsoid parameters, 37 * not the ones adjusted for spherical projections. 38 * Modified pj_datum_transform() to not attempt any datum shift via 39 * geocentric coordinates if the source *or* destination are raw ellipsoids 40 * (ie. PJD_UNKNOWN). All per PROJ bug #1602, GDAL bug #2025. 41 * 42 * Revision 1.18 2006/10/12 21:04:39 fwarmerdam 43 * Added experimental +lon_wrap argument to set a "center point" for 44 * longitude wrapping of longitude values coming out of pj_transform(). 45 * 46 * Revision 1.17 2006/09/22 23:06:24 fwarmerdam 47 * remote static start variable in pj_init (bug 1283) 48 * 49 * Revision 1.16 2004/09/08 15:23:37 warmerda 50 * added new error for unknown prime meridians 51 * 52 * Revision 1.15 2004/05/05 01:45:41 warmerda 53 * Made sword even longer. 54 * 55 * Revision 1.14 2004/05/05 01:45:00 warmerda 56 * Make sword buffer larger so long +towgs84 parameters don't get split. 57 * 58 * Revision 1.13 2003/09/16 03:46:21 warmerda 59 * dont use default ellps if any earth model info is set: bug 386 60 * 61 * Revision 1.12 2003/08/21 02:15:59 warmerda 62 * improve MAX_ARG checking 63 * 64 * Revision 1.11 2003/06/09 21:23:16 warmerda 65 * ensure start is initialized at very beginning of pj_init() 66 * 67 * Revision 1.10 2003/03/16 16:38:24 warmerda 68 * Modified get_opt() to terminate reading the definition when a new 69 * definition (a word starting with '<') is encountered, in addition to when 70 * the definition terminator '<>' is encountered, so that unterminated 71 * definitions like those in the distributed esri file will work properly. 72 * http://bugzilla.remotesensing.org/show_bug.cgi?id=302 73 * 74 * Revision 1.9 2002/12/14 20:15:02 warmerda 75 * added geocentric support, updated headers 76 * 77 */ 78 79#define PJ_LIB__ 80#include "projects.h" 81#include <stdio.h> 82#include <string.h> 83#include <errno.h> 84 85PJ_CVSID("$Id: pj_init.c,v 1.19 2007/11/26 00:21:59 fwarmerdam Exp $"); 86 87extern FILE *pj_open_lib(char *, char *); 88 89/************************************************************************/ 90/* get_opt() */ 91/************************************************************************/ 92static paralist * 93get_opt(paralist **start, FILE *fid, char *name, paralist *next) { 94 char sword[302], *word = sword+1; 95 int first = 1, len, c; 96 97 len = strlen(name); 98 *sword = 't'; 99 while (fscanf(fid, "%300s", word) == 1) { 100 if (*word == '#') /* skip comments */ 101 while((c = fgetc(fid)) != EOF && c != '\n') ; 102 else if (*word == '<') { /* control name */ 103 if (first && !strncmp(name, word + 1, len) 104 && word[len + 1] == '>') 105 first = 0; 106 else if (!first && *word == '<') { 107 while((c = fgetc(fid)) != EOF && c != '\n') ; 108 break; 109 } 110 } else if (!first && !pj_param(*start, sword).i) { 111 /* don't default ellipse if datum, ellps or any earth model 112 information is set. */ 113 if( strncmp(word,"ellps=",6) != 0 114 || (!pj_param(*start, "tdatum").i 115 && !pj_param(*start, "tellps").i 116 && !pj_param(*start, "ta").i 117 && !pj_param(*start, "tb").i 118 && !pj_param(*start, "trf").i 119 && !pj_param(*start, "tf").i) ) 120 { 121 next = next->next = pj_mkparam(word); 122 } 123 } 124 } 125 126 if (errno == 25) 127 errno = 0; 128 return next; 129} 130 131/************************************************************************/ 132/* get_defaults() */ 133/************************************************************************/ 134static paralist * 135get_defaults(paralist **start, paralist *next, char *name) { 136 FILE *fid; 137 138 fid = pj_open_lib("proj_def.dat", "rt"); 139 if (fid) { 140 next = get_opt(start, fid, "general", next); 141 rewind(fid); 142 next = get_opt(start, fid, name, next); 143 (void)fclose(fid); 144 } 145 if (errno) 146 errno = 0; /* don't care if can't open file */ 147 return next; 148} 149 150/************************************************************************/ 151/* get_init() */ 152/************************************************************************/ 153static paralist * 154get_init(paralist **start, paralist *next, char *name) { 155 char fname[MAX_PATH_FILENAME+ID_TAG_MAX+3], *opt; 156 FILE *fid; 157 158 (void)strncpy(fname, name, MAX_PATH_FILENAME + ID_TAG_MAX + 1); 159 opt = strrchr(fname, ':'); 160 if (opt) 161 *opt++ = '\0'; 162 else { pj_errno = -3; return(0); } 163 fid = pj_open_lib(fname, "rt"); 164 if (fid) 165 next = get_opt(start, fid, opt, next); 166 else 167 return(0); 168 (void)fclose(fid); 169 if (errno == 25) 170 errno = 0; /* unknown problem with some sys errno<-25 */ 171 return next; 172} 173 174/************************************************************************/ 175/* pj_init_plus() */ 176/* */ 177/* Same as pj_init() except it takes one argument string with */ 178/* individual arguments preceeded by '+', such as "+proj=utm */ 179/* +zone=11 +ellps=WGS84". */ 180/************************************************************************/ 181 182PJ * 183pj_init_plus( const char *definition ) 184 185{ 186#define MAX_ARG 200 187 char *argv[MAX_ARG]; 188 char *defn_copy; 189 int argc = 0, i; 190 PJ *result; 191 192 /* make a copy that we can manipulate */ 193 defn_copy = (char *) pj_malloc( strlen(definition)+1 ); 194 strcpy( defn_copy, definition ); 195 196 /* split into arguments based on '+' and trim white space */ 197 198 for( i = 0; defn_copy[i] != '\0'; i++ ) 199 { 200 switch( defn_copy[i] ) 201 { 202 case '+': 203 if( i == 0 || defn_copy[i-1] == '\0' ) 204 { 205 if( argc+1 == MAX_ARG ) 206 { 207 pj_errno = -44; 208 return NULL; 209 } 210 211 argv[argc++] = defn_copy + i + 1; 212 } 213 break; 214 215 case ' ': 216 case '\t': 217 case '\n': 218 defn_copy[i] = '\0'; 219 break; 220 221 default: 222 /* do nothing */; 223 } 224 } 225 226 /* perform actual initialization */ 227 result = pj_init( argc, argv ); 228 229 pj_dalloc( defn_copy ); 230 231 return result; 232} 233 234/************************************************************************/ 235/* pj_init() */ 236/* */ 237/* Main entry point for initialing a PJ projections */ 238/* definition. Note that the projection specific function is */ 239/* called to do the initial allocation so it can be created */ 240/* large enough to hold projection specific parameters. */ 241/************************************************************************/ 242 243PJ * 244pj_init(int argc, char **argv) { 245 char *s, *name; 246 paralist *start = NULL; 247 PJ *(*proj)(PJ *); 248 paralist *curr = 0; 249 int i; 250 PJ *PIN = 0; 251 252 errno = pj_errno = 0; 253 start = NULL; 254 255 /* put arguments into internal linked list */ 256 if (argc <= 0) { pj_errno = -1; goto bum_call; } 257 for (i = 0; i < argc; ++i) 258 if (i) 259 curr = curr->next = pj_mkparam(argv[i]); 260 else 261 start = curr = pj_mkparam(argv[i]); 262 if (pj_errno) goto bum_call; 263 264 /* check if +init present */ 265 if (pj_param(start, "tinit").i) { 266 paralist *last = curr; 267 268 if (!(curr = get_init(&start, curr, pj_param(start, "sinit").s))) 269 goto bum_call; 270 if (curr == last) { pj_errno = -2; goto bum_call; } 271 } 272 273 /* find projection selection */ 274 if (!(name = pj_param(start, "sproj").s)) 275 { pj_errno = -4; goto bum_call; } 276 for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ; 277 if (!s) { pj_errno = -5; goto bum_call; } 278 279 /* set defaults, unless inhibited */ 280 if (!pj_param(start, "bno_defs").i) 281 get_defaults(&start, curr, name); 282 proj = (PJ *(*)(PJ *)) pj_list[i].proj; 283 284 /* allocate projection structure */ 285 if (!(PIN = (*proj)(0))) goto bum_call; 286 PIN->params = start; 287 PIN->is_latlong = 0; 288 PIN->is_geocent = 0; 289 PIN->long_wrap_center = 0.0; 290 291 /* set datum parameters */ 292 if (pj_datum_set(start, PIN)) goto bum_call; 293 294 /* set ellipsoid/sphere parameters */ 295 if (pj_ell_set(start, &PIN->a, &PIN->es)) goto bum_call; 296 297 PIN->a_orig = PIN->a; 298 PIN->es_orig = PIN->es; 299 300 PIN->e = sqrt(PIN->es); 301 PIN->ra = 1. / PIN->a; 302 PIN->one_es = 1. - PIN->es; 303 if (PIN->one_es == 0.) { pj_errno = -6; goto bum_call; } 304 PIN->rone_es = 1./PIN->one_es; 305 306 /* Now that we have ellipse information check for WGS84 datum */ 307 if( PIN->datum_type == PJD_3PARAM 308 && PIN->datum_params[0] == 0.0 309 && PIN->datum_params[1] == 0.0 310 && PIN->datum_params[2] == 0.0 311 && PIN->a == 6378137.0 312 && ABS(PIN->es - 0.006694379990) < 0.000000000050 )/*WGS84/GRS80*/ 313 { 314 PIN->datum_type = PJD_WGS84; 315 } 316 317 /* set PIN->geoc coordinate system */ 318 PIN->geoc = (PIN->es && pj_param(start, "bgeoc").i); 319 320 /* over-ranging flag */ 321 PIN->over = pj_param(start, "bover").i; 322 323 /* longitude center for wrapping */ 324 PIN->long_wrap_center = pj_param(start, "rlon_wrap").f; 325 326 /* central meridian */ 327 PIN->lam0=pj_param(start, "rlon_0").f; 328 329 /* central latitude */ 330 PIN->phi0 = pj_param(start, "rlat_0").f; 331 332 /* false easting and northing */ 333 PIN->x0 = pj_param(start, "dx_0").f; 334 PIN->y0 = pj_param(start, "dy_0").f; 335 336 /* general scaling factor */ 337 if (pj_param(start, "tk_0").i) 338 PIN->k0 = pj_param(start, "dk_0").f; 339 else if (pj_param(start, "tk").i) 340 PIN->k0 = pj_param(start, "dk").f; 341 else 342 PIN->k0 = 1.; 343 if (PIN->k0 <= 0.) { 344 pj_errno = -31; 345 goto bum_call; 346 } 347 348 /* set units */ 349 s = 0; 350 name = pj_param(start, "sunits").s; 351 if (name) { 352 for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i) ; 353 if (!s) { pj_errno = -7; goto bum_call; } 354 s = pj_units[i].to_meter; 355 } 356 if (s || (s = pj_param(start, "sto_meter").s)) { 357 PIN->to_meter = strtod(s, &s); 358 if (*s == '/') /* ratio number */ 359 PIN->to_meter /= strtod(++s, 0); 360 PIN->fr_meter = 1. / PIN->to_meter; 361 } else 362 PIN->to_meter = PIN->fr_meter = 1.; 363 364 /* prime meridian */ 365 s = 0; 366 name = pj_param(start, "spm").s; 367 if (name) { 368 const char *value = NULL; 369 char *next_str = NULL; 370 371 for (i = 0; pj_prime_meridians[i].id != NULL; ++i ) 372 { 373 if( strcmp(name,pj_prime_meridians[i].id) == 0 ) 374 { 375 value = pj_prime_meridians[i].defn; 376 break; 377 } 378 } 379 380 if( value == NULL 381 && (dmstor(name,&next_str) != 0.0 || *name == '0') 382 && *next_str == '\0' ) 383 value = name; 384 385 if (!value) { pj_errno = -46; goto bum_call; } 386 PIN->from_greenwich = dmstor(value,NULL); 387 } 388 else 389 PIN->from_greenwich = 0.0; 390 391 /* projection specific initialization */ 392 if (!(PIN = (*proj)(PIN)) || errno || pj_errno) { 393bum_call: /* cleanup error return */ 394 if (!pj_errno) 395 pj_errno = errno; 396 if (PIN) 397 pj_free(PIN); 398 else 399 for ( ; start; start = curr) { 400 curr = start->next; 401 pj_dalloc(start); 402 } 403 PIN = 0; 404 } 405 return PIN; 406} 407 408/************************************************************************/ 409/* pj_free() */ 410/* */ 411/* This is the application callable entry point for destroying */ 412/* a projection definition. It does work generic to all */ 413/* projection types, and then calls the projection specific */ 414/* free function (P->pfree()) to do local work. This maps to */ 415/* the FREEUP code in the individual projection source files. */ 416/************************************************************************/ 417 418void 419pj_free(PJ *P) { 420 if (P) { 421 paralist *t, *n; 422 423 /* free parameter list elements */ 424 for (t = P->params; t; t = n) { 425 n = t->next; 426 pj_dalloc(t); 427 } 428 429 /* free projection parameters */ 430 P->pfree(P); 431 } 432} 433 434