PageRenderTime 37ms CodeModel.GetById 15ms app.highlight 17ms RepoModel.GetById 1ms app.codeStats 0ms

/Proj4/pj_init.c

http://github.com/route-me/route-me
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