1 /******************************************************************************
2 * $Id: geotiff_proj4.c 2594 2014-12-27 16:46:35Z rouault $
3 *
4 * Project: libgeotiff
5 * Purpose: Code to convert a normalized GeoTIFF definition into a PROJ.4
6 * (OGDI) compatible projection string.
7 * Author: Frank Warmerdam, warmerdam@pobox.com
8 *
9 ******************************************************************************
10 * Copyright (c) 1999, Frank Warmerdam
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ******************************************************************************
30 */
31
32 #include "cpl_serv.h"
33 #include "geotiff.h"
34 #include "geo_normalize.h"
35 #include "geovalues.h"
36
37 /************************************************************************/
38 /* OSRProj4Tokenize() */
39 /* */
40 /* Custom tokenizing function for PROJ.4 strings. The main */
41 /* reason we can't just use CSLTokenizeString is to handle */
42 /* strings with a + sign in the exponents of parameter values. */
43 /************************************************************************/
44
OSRProj4Tokenize(const char * pszFull)45 static char **OSRProj4Tokenize( const char *pszFull )
46
47 {
48 char *pszStart = NULL;
49 char *pszFullWrk;
50 char **papszTokens;
51 int i;
52 int nTokens = 0;
53 static const int nMaxTokens = 200;
54
55 if( pszFull == NULL )
56 return NULL;
57
58 papszTokens = (char **) calloc(sizeof(char*),nMaxTokens);
59
60 pszFullWrk = CPLStrdup(pszFull);
61
62 for( i=0; pszFullWrk[i] != '\0' && nTokens != nMaxTokens-1; i++ )
63 {
64 switch( pszFullWrk[i] )
65 {
66 case '+':
67 if( i == 0 || pszFullWrk[i-1] == '\0' )
68 {
69 if( pszStart != NULL )
70 {
71 if( strstr(pszStart,"=") != NULL )
72 {
73 papszTokens[nTokens++] = CPLStrdup(pszStart);
74 }
75 else
76 {
77 char szAsBoolean[100];
78 strncpy( szAsBoolean,pszStart, sizeof(szAsBoolean)-1-4);
79 szAsBoolean[sizeof(szAsBoolean)-1-4] = '\0';
80 strcat( szAsBoolean,"=yes" );
81 papszTokens[nTokens++] = CPLStrdup(szAsBoolean);
82 }
83 }
84 pszStart = pszFullWrk + i + 1;
85 }
86 break;
87
88 case ' ':
89 case '\t':
90 case '\n':
91 pszFullWrk[i] = '\0';
92 break;
93
94 default:
95 break;
96 }
97 }
98
99 if( pszStart != NULL && strlen(pszStart) > 0 )
100 {
101 if (nTokens != 199)
102 papszTokens[nTokens++] = CPLStrdup(pszStart);
103 }
104
105 CPLFree( pszFullWrk );
106
107 return papszTokens;
108 }
109
110
111 /************************************************************************/
112 /* OSR_GSV() */
113 /************************************************************************/
114
OSR_GSV(char ** papszNV,const char * pszField)115 static const char *OSR_GSV( char **papszNV, const char * pszField )
116
117 {
118 size_t field_len = strlen(pszField);
119 int i;
120
121 if( !papszNV )
122 return NULL;
123
124 for( i = 0; papszNV[i] != NULL; i++ )
125 {
126 if( EQUALN(papszNV[i],pszField,field_len) )
127 {
128 if( papszNV[i][field_len] == '=' )
129 return papszNV[i] + field_len + 1;
130
131 if( strlen(papszNV[i]) == field_len )
132 return "";
133 }
134 }
135
136 return NULL;
137 }
138
139 /************************************************************************/
140 /* OSR_GDV() */
141 /* */
142 /* Fetch a particular parameter out of the parameter list, or */
143 /* the indicated default if it isn't available. This is a */
144 /* helper function for importFromProj4(). */
145 /************************************************************************/
146
OSR_GDV(char ** papszNV,const char * pszField,double dfDefaultValue)147 static double OSR_GDV( char **papszNV, const char * pszField,
148 double dfDefaultValue )
149
150 {
151 const char *pszValue = OSR_GSV( papszNV, pszField );
152
153 /* special hack to use k_0 if available. */
154 if( pszValue == NULL && EQUAL(pszField,"k") )
155 return OSR_GDV( papszNV, "k_0", dfDefaultValue );
156
157 if( pszValue == NULL )
158 return dfDefaultValue;
159 else
160 return GTIFAtof(pszValue);
161 }
162
163 /************************************************************************/
164 /* OSRFreeStringList() */
165 /************************************************************************/
166
OSRFreeStringList(char ** list)167 static void OSRFreeStringList( char ** list )
168
169 {
170 int i;
171
172 for( i = 0; list != NULL && list[i] != NULL; i++ )
173 free( list[i] );
174 free(list);
175 }
176
177
178 /************************************************************************/
179 /* GTIFSetFromProj4() */
180 /************************************************************************/
181
GTIFSetFromProj4(GTIF * gtif,const char * proj4)182 int GTIFSetFromProj4( GTIF *gtif, const char *proj4 )
183
184 {
185 char **papszNV = OSRProj4Tokenize( proj4 );
186 short nSpheroid = KvUserDefined;
187 double dfSemiMajor=0.0, dfSemiMinor=0.0, dfInvFlattening=0.0;
188 int nDatum = KvUserDefined;
189 int nGCS = KvUserDefined;
190 const char *value;
191
192 /* -------------------------------------------------------------------- */
193 /* Get the ellipsoid definition. */
194 /* -------------------------------------------------------------------- */
195 value = OSR_GSV( papszNV, "ellps" );
196
197 if( value == NULL )
198 {
199 /* nothing */;
200 }
201 else if( EQUAL(value,"WGS84") )
202 nSpheroid = Ellipse_WGS_84;
203 else if( EQUAL(value,"clrk66") )
204 nSpheroid = Ellipse_Clarke_1866;
205 else if( EQUAL(value,"clrk80") )
206 nSpheroid = Ellipse_Clarke_1880;
207 else if( EQUAL(value,"GRS80") )
208 nSpheroid = Ellipse_GRS_1980;
209
210 if( nSpheroid == KvUserDefined )
211 {
212 dfSemiMajor = OSR_GDV(papszNV,"a",0.0);
213 dfSemiMinor = OSR_GDV(papszNV,"b",0.0);
214 dfInvFlattening = OSR_GDV(papszNV,"rf",0.0);
215 if( dfSemiMinor != 0.0 && dfInvFlattening == 0.0 )
216 dfInvFlattening = -1.0 / (dfSemiMinor/dfSemiMajor - 1.0);
217 }
218
219 /* -------------------------------------------------------------------- */
220 /* Get the GCS/Datum code. */
221 /* -------------------------------------------------------------------- */
222 value = OSR_GSV( papszNV, "datum" );
223
224 if( value == NULL )
225 {
226 }
227 else if( EQUAL(value,"WGS84") )
228 {
229 nGCS = GCS_WGS_84;
230 nDatum = Datum_WGS84;
231 }
232 else if( EQUAL(value,"NAD83") )
233 {
234 nGCS = GCS_NAD83;
235 nDatum = Datum_North_American_Datum_1983;
236 }
237 else if( EQUAL(value,"NAD27") )
238 {
239 nGCS = GCS_NAD27;
240 nDatum = Datum_North_American_Datum_1927;
241 }
242
243 /* -------------------------------------------------------------------- */
244 /* Operate on the basis of the projection name. */
245 /* -------------------------------------------------------------------- */
246 value = OSR_GSV(papszNV,"proj");
247
248 if( value == NULL )
249 {
250 OSRFreeStringList( papszNV );
251 return FALSE;
252 }
253
254 else if( EQUAL(value,"longlat") || EQUAL(value,"latlong") )
255 {
256 }
257
258 else if( EQUAL(value,"tmerc") )
259 {
260 GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1,
261 ModelTypeProjected);
262 GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
263 KvUserDefined );
264 GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1,
265 KvUserDefined );
266
267 GTIFKeySet(gtif, ProjCoordTransGeoKey, TYPE_SHORT, 1,
268 CT_TransverseMercator );
269
270 GTIFKeySet(gtif, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
271 OSR_GDV( papszNV, "lat_0", 0.0 ) );
272
273 GTIFKeySet(gtif, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
274 OSR_GDV( papszNV, "lon_0", 0.0 ) );
275
276 GTIFKeySet(gtif, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
277 OSR_GDV( papszNV, "k", 1.0 ) );
278
279 GTIFKeySet(gtif, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
280 OSR_GDV( papszNV, "x_0", 0.0 ) );
281
282 GTIFKeySet(gtif, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
283 OSR_GDV( papszNV, "y_0", 0.0 ) );
284 }
285
286 else if( EQUAL(value,"utm") )
287 {
288 int nZone = (int) OSR_GDV(papszNV,"zone",0);
289 const char *south = OSR_GSV(papszNV,"south");
290
291 GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1,
292 ModelTypeProjected);
293 GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
294 KvUserDefined );
295 GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1,
296 KvUserDefined );
297
298 GTIFKeySet(gtif, ProjCoordTransGeoKey, TYPE_SHORT, 1,
299 CT_TransverseMercator );
300
301 GTIFKeySet(gtif, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
302 0.0 );
303
304 GTIFKeySet(gtif, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
305 nZone * 6 - 183.0 );
306
307 GTIFKeySet(gtif, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
308 0.9996 );
309
310 GTIFKeySet(gtif, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
311 500000.0 );
312
313 if( south != NULL )
314 GTIFKeySet(gtif, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
315 10000000.0 );
316 else
317 GTIFKeySet(gtif, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
318 0.0 );
319 }
320
321 else if( EQUAL(value,"lcc")
322 && OSR_GDV(papszNV, "lat_0", 0.0 )
323 == OSR_GDV(papszNV, "lat_1", 0.0 ) )
324 {
325 GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1,
326 ModelTypeProjected);
327 GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
328 KvUserDefined );
329 GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1,
330 KvUserDefined );
331
332 GTIFKeySet(gtif, ProjCoordTransGeoKey, TYPE_SHORT, 1,
333 CT_LambertConfConic_1SP );
334
335 GTIFKeySet(gtif, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
336 OSR_GDV( papszNV, "lat_0", 0.0 ) );
337
338 GTIFKeySet(gtif, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
339 OSR_GDV( papszNV, "lon_0", 0.0 ) );
340
341 GTIFKeySet(gtif, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
342 OSR_GDV( papszNV, "k", 1.0 ) );
343
344 GTIFKeySet(gtif, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
345 OSR_GDV( papszNV, "x_0", 0.0 ) );
346
347 GTIFKeySet(gtif, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
348 OSR_GDV( papszNV, "y_0", 0.0 ) );
349 }
350
351 else if( EQUAL(value,"lcc") )
352 {
353 GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1,
354 ModelTypeProjected);
355 GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
356 KvUserDefined );
357 GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1,
358 KvUserDefined );
359
360 GTIFKeySet(gtif, ProjCoordTransGeoKey, TYPE_SHORT, 1,
361 CT_LambertConfConic_2SP );
362
363 GTIFKeySet(gtif, ProjFalseOriginLatGeoKey, TYPE_DOUBLE, 1,
364 OSR_GDV( papszNV, "lat_0", 0.0 ) );
365
366 GTIFKeySet(gtif, ProjFalseOriginLongGeoKey, TYPE_DOUBLE, 1,
367 OSR_GDV( papszNV, "lon_0", 0.0 ) );
368
369 GTIFKeySet(gtif, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
370 OSR_GDV( papszNV, "lat_1", 0.0 ) );
371
372 GTIFKeySet(gtif, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
373 OSR_GDV( papszNV, "lat_2", 0.0 ) );
374
375 GTIFKeySet(gtif, ProjFalseOriginEastingGeoKey, TYPE_DOUBLE, 1,
376 OSR_GDV( papszNV, "x_0", 0.0 ) );
377
378 GTIFKeySet(gtif, ProjFalseOriginNorthingGeoKey, TYPE_DOUBLE, 1,
379 OSR_GDV( papszNV, "y_0", 0.0 ) );
380 }
381
382 #ifdef notdef
383 else if( EQUAL(value,"bonne") )
384 {
385 SetBonne( OSR_GDV( papszNV, "lat_1", 0.0 ),
386 OSR_GDV( papszNV, "lon_0", 0.0 ),
387 OSR_GDV( papszNV, "x_0", 0.0 ),
388 OSR_GDV( papszNV, "y_0", 0.0 ) );
389 }
390
391 else if( EQUAL(value,"cass") )
392 {
393 SetCS( OSR_GDV( papszNV, "lat_0", 0.0 ),
394 OSR_GDV( papszNV, "lon_0", 0.0 ),
395 OSR_GDV( papszNV, "x_0", 0.0 ),
396 OSR_GDV( papszNV, "y_0", 0.0 ) );
397 }
398
399 else if( EQUAL(value,"nzmg") )
400 {
401 SetNZMG( OSR_GDV( papszNV, "lat_0", -41.0 ),
402 OSR_GDV( papszNV, "lon_0", 173.0 ),
403 OSR_GDV( papszNV, "x_0", 2510000.0 ),
404 OSR_GDV( papszNV, "y_0", 6023150.0 ) );
405 }
406
407 else if( EQUAL(value,"cea") )
408 {
409 SetCEA( OSR_GDV( papszNV, "lat_ts", 0.0 ),
410 OSR_GDV( papszNV, "lon_0", 0.0 ),
411 OSR_GDV( papszNV, "x_0", 0.0 ),
412 OSR_GDV( papszNV, "y_0", 0.0 ) );
413 }
414
415 else if( EQUAL(value,"merc") /* 2SP form */
416 && OSR_GDV(papszNV, "lat_ts", 1000.0) < 999.0 )
417 {
418 SetMercator2SP( OSR_GDV( papszNV, "lat_ts", 0.0 ),
419 0.0,
420 OSR_GDV( papszNV, "lon_0", 0.0 ),
421 OSR_GDV( papszNV, "x_0", 0.0 ),
422 OSR_GDV( papszNV, "y_0", 0.0 ) );
423 }
424
425 else if( EQUAL(value,"merc") ) /* 1SP form */
426 {
427 SetMercator( 0.0,
428 OSR_GDV( papszNV, "lon_0", 0.0 ),
429 OSR_GDV( papszNV, "k", 1.0 ),
430 OSR_GDV( papszNV, "x_0", 0.0 ),
431 OSR_GDV( papszNV, "y_0", 0.0 ) );
432 }
433
434 else if( EQUAL(value,"stere")
435 && ABS(OSR_GDV( papszNV, "lat_0", 0.0 ) - 90) < 0.001 )
436 {
437 SetPS( OSR_GDV( papszNV, "lat_ts", 90.0 ),
438 OSR_GDV( papszNV, "lon_0", 0.0 ),
439 OSR_GDV( papszNV, "k", 1.0 ),
440 OSR_GDV( papszNV, "x_0", 0.0 ),
441 OSR_GDV( papszNV, "y_0", 0.0 ) );
442 }
443
444 else if( EQUAL(value,"stere")
445 && ABS(OSR_GDV( papszNV, "lat_0", 0.0 ) + 90) < 0.001 )
446 {
447 SetPS( OSR_GDV( papszNV, "lat_ts", -90.0 ),
448 OSR_GDV( papszNV, "lon_0", 0.0 ),
449 OSR_GDV( papszNV, "k", 1.0 ),
450 OSR_GDV( papszNV, "x_0", 0.0 ),
451 OSR_GDV( papszNV, "y_0", 0.0 ) );
452 }
453
454 else if( EQUALN(value,"stere",5) /* mostly sterea */
455 && CSLFetchNameValue(papszNV,"k") != NULL )
456 {
457 SetOS( OSR_GDV( papszNV, "lat_0", 0.0 ),
458 OSR_GDV( papszNV, "lon_0", 0.0 ),
459 OSR_GDV( papszNV, "k", 1.0 ),
460 OSR_GDV( papszNV, "x_0", 0.0 ),
461 OSR_GDV( papszNV, "y_0", 0.0 ) );
462 }
463
464 else if( EQUAL(value,"stere") )
465 {
466 SetStereographic( OSR_GDV( papszNV, "lat_0", 0.0 ),
467 OSR_GDV( papszNV, "lon_0", 0.0 ),
468 1.0,
469 OSR_GDV( papszNV, "x_0", 0.0 ),
470 OSR_GDV( papszNV, "y_0", 0.0 ) );
471 }
472
473 else if( EQUAL(value,"eqc") )
474 {
475 if( OSR_GDV( papszNV, "lat_0", 0.0 ) != OSR_GDV( papszNV, "lat_ts", 0.0 ) )
476 SetEquirectangular2( OSR_GDV( papszNV, "lat_0", 0.0 ),
477 OSR_GDV( papszNV, "lon_0", 0.0 )+dfFromGreenwich,
478 OSR_GDV( papszNV, "lat_ts", 0.0 ),
479 OSR_GDV( papszNV, "x_0", 0.0 ),
480 OSR_GDV( papszNV, "y_0", 0.0 ) );
481 else
482 SetEquirectangular( OSR_GDV( papszNV, "lat_ts", 0.0 ),
483 OSR_GDV( papszNV, "lon_0", 0.0 )+dfFromGreenwich,
484 OSR_GDV( papszNV, "x_0", 0.0 ),
485 OSR_GDV( papszNV, "y_0", 0.0 ) );
486 }
487
488 else if( EQUAL(value,"glabsgm") )
489 {
490 SetGaussLabordeReunion( OSR_GDV( papszNV, "lat_0", -21.116666667 ),
491 OSR_GDV( papszNV, "lon_0", 55.53333333309)+dfFromGreenwich,
492 OSR_GDV( papszNV, "k_0", 1.0 ),
493 OSR_GDV( papszNV, "x_0", 160000.000 ),
494 OSR_GDV( papszNV, "y_0", 50000.000 ) );
495 }
496
497 else if( EQUAL(value,"gnom") )
498 {
499 SetGnomonic( OSR_GDV( papszNV, "lat_0", 0.0 ),
500 OSR_GDV( papszNV, "lon_0", 0.0 ),
501 OSR_GDV( papszNV, "x_0", 0.0 ),
502 OSR_GDV( papszNV, "y_0", 0.0 ) );
503 }
504
505 else if( EQUAL(value,"ortho") )
506 {
507 SetOrthographic( OSR_GDV( papszNV, "lat_0", 0.0 ),
508 OSR_GDV( papszNV, "lon_0", 0.0 ),
509 OSR_GDV( papszNV, "x_0", 0.0 ),
510 OSR_GDV( papszNV, "y_0", 0.0 ) );
511 }
512
513 else if( EQUAL(value,"laea") )
514 {
515 SetLAEA( OSR_GDV( papszNV, "lat_0", 0.0 ),
516 OSR_GDV( papszNV, "lon_0", 0.0 ),
517 OSR_GDV( papszNV, "x_0", 0.0 ),
518 OSR_GDV( papszNV, "y_0", 0.0 ) );
519 }
520
521 else if( EQUAL(value,"aeqd") )
522 {
523 SetAE( OSR_GDV( papszNV, "lat_0", 0.0 ),
524 OSR_GDV( papszNV, "lon_0", 0.0 ),
525 OSR_GDV( papszNV, "x_0", 0.0 ),
526 OSR_GDV( papszNV, "y_0", 0.0 ) );
527 }
528
529 else if( EQUAL(value,"eqdc") )
530 {
531 SetEC( OSR_GDV( papszNV, "lat_1", 0.0 ),
532 OSR_GDV( papszNV, "lat_2", 0.0 ),
533 OSR_GDV( papszNV, "lat_0", 0.0 ),
534 OSR_GDV( papszNV, "lon_0", 0.0 ),
535 OSR_GDV( papszNV, "x_0", 0.0 ),
536 OSR_GDV( papszNV, "y_0", 0.0 ) );
537 }
538
539 else if( EQUAL(value,"mill") )
540 {
541 SetMC( OSR_GDV( papszNV, "lat_0", 0.0 ),
542 OSR_GDV( papszNV, "lon_0", 0.0 ),
543 OSR_GDV( papszNV, "x_0", 0.0 ),
544 OSR_GDV( papszNV, "y_0", 0.0 ) );
545 }
546
547 else if( EQUAL(value,"moll") )
548 {
549 SetMollweide( OSR_GDV( papszNV, "lon_0", 0.0 ),
550 OSR_GDV( papszNV, "x_0", 0.0 ),
551 OSR_GDV( papszNV, "y_0", 0.0 ) );
552 }
553
554 else if( EQUAL(value,"eck4") )
555 {
556 SetEckertIV( OSR_GDV( papszNV, "lon_0", 0.0 ),
557 OSR_GDV( papszNV, "x_0", 0.0 ),
558 OSR_GDV( papszNV, "y_0", 0.0 ) );
559 }
560
561 else if( EQUAL(value,"eck6") )
562 {
563 SetEckertVI( OSR_GDV( papszNV, "lon_0", 0.0 ),
564 OSR_GDV( papszNV, "x_0", 0.0 ),
565 OSR_GDV( papszNV, "y_0", 0.0 ) );
566 }
567
568 else if( EQUAL(value,"poly") )
569 {
570 SetPolyconic( OSR_GDV( papszNV, "lat_0", 0.0 ),
571 OSR_GDV( papszNV, "lon_0", 0.0 ),
572 OSR_GDV( papszNV, "x_0", 0.0 ),
573 OSR_GDV( papszNV, "y_0", 0.0 ) );
574 }
575
576 else if( EQUAL(value,"aea") )
577 {
578 SetACEA( OSR_GDV( papszNV, "lat_1", 0.0 ),
579 OSR_GDV( papszNV, "lat_2", 0.0 ),
580 OSR_GDV( papszNV, "lat_0", 0.0 ),
581 OSR_GDV( papszNV, "lon_0", 0.0 ),
582 OSR_GDV( papszNV, "x_0", 0.0 ),
583 OSR_GDV( papszNV, "y_0", 0.0 ) );
584 }
585
586 else if( EQUAL(value,"robin") )
587 {
588 SetRobinson( OSR_GDV( papszNV, "lon_0", 0.0 ),
589 OSR_GDV( papszNV, "x_0", 0.0 ),
590 OSR_GDV( papszNV, "y_0", 0.0 ) );
591 }
592
593 else if( EQUAL(value,"vandg") )
594 {
595 SetVDG( OSR_GDV( papszNV, "lon_0", 0.0 ),
596 OSR_GDV( papszNV, "x_0", 0.0 ),
597 OSR_GDV( papszNV, "y_0", 0.0 ) );
598 }
599
600 else if( EQUAL(value,"sinu") )
601 {
602 SetSinusoidal( OSR_GDV( papszNV, "lon_0", 0.0 ),
603 OSR_GDV( papszNV, "x_0", 0.0 ),
604 OSR_GDV( papszNV, "y_0", 0.0 ) );
605 }
606
607 else if( EQUAL(value,"gall") )
608 {
609 SetGS( OSR_GDV( papszNV, "lon_0", 0.0 ),
610 OSR_GDV( papszNV, "x_0", 0.0 ),
611 OSR_GDV( papszNV, "y_0", 0.0 ) );
612 }
613
614 else if( EQUAL(value,"goode") )
615 {
616 SetGH( OSR_GDV( papszNV, "lon_0", 0.0 ),
617 OSR_GDV( papszNV, "x_0", 0.0 ),
618 OSR_GDV( papszNV, "y_0", 0.0 ) );
619 }
620
621 else if( EQUAL(value,"geos") )
622 {
623 SetGEOS( OSR_GDV( papszNV, "lon_0", 0.0 ),
624 OSR_GDV( papszNV, "h", 35785831.0 ),
625 OSR_GDV( papszNV, "x_0", 0.0 ),
626 OSR_GDV( papszNV, "y_0", 0.0 ) );
627 }
628
629 else if( EQUAL(value,"lcc") )
630 {
631 if( OSR_GDV(papszNV, "lat_0", 0.0 )
632 == OSR_GDV(papszNV, "lat_1", 0.0 ) )
633 {
634 /* 1SP form */
635 SetLCC1SP( OSR_GDV( papszNV, "lat_0", 0.0 ),
636 OSR_GDV( papszNV, "lon_0", 0.0 ),
637 OSR_GDV( papszNV, "k_0", 1.0 ),
638 OSR_GDV( papszNV, "x_0", 0.0 ),
639 OSR_GDV( papszNV, "y_0", 0.0 ) );
640 }
641 else
642 {
643 /* 2SP form */
644 SetLCC( OSR_GDV( papszNV, "lat_1", 0.0 ),
645 OSR_GDV( papszNV, "lat_2", 0.0 ),
646 OSR_GDV( papszNV, "lat_0", 0.0 ),
647 OSR_GDV( papszNV, "lon_0", 0.0 ),
648 OSR_GDV( papszNV, "x_0", 0.0 ),
649 OSR_GDV( papszNV, "y_0", 0.0 ) );
650 }
651 }
652
653 else if( EQUAL(value,"omerc") )
654 {
655 SetHOM( OSR_GDV( papszNV, "lat_0", 0.0 ),
656 OSR_GDV( papszNV, "lonc", 0.0 ),
657 OSR_GDV( papszNV, "alpha", 0.0 ),
658 0.0, /* ??? */
659 OSR_GDV( papszNV, "k", 1.0 ),
660 OSR_GDV( papszNV, "x_0", 0.0 ),
661 OSR_GDV( papszNV, "y_0", 0.0 ) );
662 }
663
664 else if( EQUAL(value,"somerc") )
665 {
666 SetHOM( OSR_GDV( papszNV, "lat_0", 0.0 ),
667 OSR_GDV( papszNV, "lon_0", 0.0 ),
668 90.0, 90.0,
669 OSR_GDV( papszNV, "k", 1.0 ),
670 OSR_GDV( papszNV, "x_0", 0.0 ),
671 OSR_GDV( papszNV, "y_0", 0.0 ) );
672 }
673
674 else if( EQUAL(value,"krovak") )
675 {
676 SetKrovak( OSR_GDV( papszNV, "lat_0", 0.0 ),
677 OSR_GDV( papszNV, "lon_0", 0.0 ),
678 OSR_GDV( papszNV, "alpha", 0.0 ),
679 0.0, /* pseudo_standard_parallel_1 */
680 OSR_GDV( papszNV, "k", 1.0 ),
681 OSR_GDV( papszNV, "x_0", 0.0 ),
682 OSR_GDV( papszNV, "y_0", 0.0 ) );
683 }
684
685 else if( EQUAL(value, "iwm_p") )
686 {
687 SetIWMPolyconic( OSR_GDV( papszNV, "lat_1", 0.0 ),
688 OSR_GDV( papszNV, "lat_2", 0.0 ),
689 OSR_GDV( papszNV, "lon_0", 0.0 ),
690 OSR_GDV( papszNV, "x_0", 0.0 ),
691 OSR_GDV( papszNV, "y_0", 0.0 ) );
692 }
693
694 else if( EQUAL(value, "wag1") )
695 {
696 SetWagner( 1, 0.0,
697 OSR_GDV( papszNV, "x_0", 0.0 ),
698 OSR_GDV( papszNV, "y_0", 0.0 ) );
699 }
700
701 else if( EQUAL(value, "wag2") )
702 {
703 SetWagner( 2, 0.0,
704 OSR_GDV( papszNV, "x_0", 0.0 ),
705 OSR_GDV( papszNV, "y_0", 0.0 ) );
706 }
707
708 else if( EQUAL(value, "wag3") )
709 {
710 SetWagner( 3,
711 OSR_GDV( papszNV, "lat_ts", 0.0 ),
712 OSR_GDV( papszNV, "x_0", 0.0 ),
713 OSR_GDV( papszNV, "y_0", 0.0 ) );
714 }
715
716 else if( EQUAL(value, "wag1") )
717 {
718 SetWagner( 4, 0.0,
719 OSR_GDV( papszNV, "x_0", 0.0 ),
720 OSR_GDV( papszNV, "y_0", 0.0 ) );
721 }
722
723 else if( EQUAL(value, "wag1") )
724 {
725 SetWagner( 5, 0.0,
726 OSR_GDV( papszNV, "x_0", 0.0 ),
727 OSR_GDV( papszNV, "y_0", 0.0 ) );
728 }
729
730 else if( EQUAL(value, "wag1") )
731 {
732 SetWagner( 6, 0.0,
733 OSR_GDV( papszNV, "x_0", 0.0 ),
734 OSR_GDV( papszNV, "y_0", 0.0 ) );
735 }
736
737 else if( EQUAL(value, "wag1") )
738 {
739 SetWagner( 7, 0.0,
740 OSR_GDV( papszNV, "x_0", 0.0 ),
741 OSR_GDV( papszNV, "y_0", 0.0 ) );
742 }
743
744 else if( EQUAL(value,"tpeqd") )
745 {
746 SetTPED( OSR_GDV( papszNV, "lat_1", 0.0 ),
747 OSR_GDV( papszNV, "lon_1", 0.0 ),
748 OSR_GDV( papszNV, "lat_2", 0.0 ),
749 OSR_GDV( papszNV, "lon_2", 0.0 ),
750 OSR_GDV( papszNV, "x_0", 0.0 ),
751 OSR_GDV( papszNV, "y_0", 0.0 ) );
752 }
753 #endif
754 else
755 {
756 /* unsupported coordinate system */
757 OSRFreeStringList( papszNV );
758 return FALSE;
759 }
760
761 /* -------------------------------------------------------------------- */
762 /* Write the GCS if we have it, otherwise write the datum. */
763 /* -------------------------------------------------------------------- */
764 if( nGCS != KvUserDefined )
765 {
766 GTIFKeySet( gtif, GeographicTypeGeoKey, TYPE_SHORT,
767 1, nGCS );
768 }
769 else
770 {
771 GTIFKeySet( gtif, GeographicTypeGeoKey, TYPE_SHORT, 1,
772 KvUserDefined );
773 GTIFKeySet( gtif, GeogGeodeticDatumGeoKey, TYPE_SHORT,
774 1, nDatum );
775 }
776
777 /* -------------------------------------------------------------------- */
778 /* Write the ellipsoid if we don't know the GCS. */
779 /* -------------------------------------------------------------------- */
780 if( nGCS == KvUserDefined )
781 {
782 if( nSpheroid != KvUserDefined )
783 GTIFKeySet( gtif, GeogEllipsoidGeoKey, TYPE_SHORT, 1,
784 nSpheroid );
785 else
786 {
787 GTIFKeySet( gtif, GeogEllipsoidGeoKey, TYPE_SHORT, 1,
788 KvUserDefined );
789 GTIFKeySet( gtif, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
790 dfSemiMajor );
791 if( dfInvFlattening == 0.0 )
792 GTIFKeySet( gtif, GeogSemiMinorAxisGeoKey, TYPE_DOUBLE, 1,
793 dfSemiMajor );
794 else
795 GTIFKeySet( gtif, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
796 dfInvFlattening );
797 }
798
799 }
800
801 /* -------------------------------------------------------------------- */
802 /* Linear units translation */
803 /* -------------------------------------------------------------------- */
804 value = OSR_GSV( papszNV, "units" );
805
806 if( value == NULL )
807 {
808 value = OSR_GSV( papszNV, "to_meter" );
809 if( value )
810 {
811 GTIFKeySet( gtif, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
812 KvUserDefined );
813 GTIFKeySet( gtif, ProjLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
814 GTIFAtof(value) );
815 }
816 }
817 else if( EQUAL(value,"meter") || EQUAL(value,"m") )
818 {
819 GTIFKeySet( gtif, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
820 Linear_Meter );
821 }
822 else if( EQUAL(value,"us-ft") )
823 {
824 GTIFKeySet( gtif, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
825 Linear_Foot_US_Survey );
826 }
827 else if( EQUAL(value,"ft") )
828 {
829 GTIFKeySet( gtif, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
830 Linear_Foot );
831 }
832
833
834 OSRFreeStringList( papszNV );
835
836 return TRUE;
837 }
838
839 /************************************************************************/
840 /* GTIFGetProj4Defn() */
841 /************************************************************************/
842
GTIFGetProj4Defn(GTIFDefn * psDefn)843 char * GTIFGetProj4Defn( GTIFDefn * psDefn )
844
845 {
846 char szProjection[512];
847 char szUnits[64];
848 double dfFalseEasting, dfFalseNorthing;
849
850 if( psDefn == NULL || !psDefn->DefnSet )
851 return CPLStrdup("");
852
853 szProjection[0] = '\0';
854
855 /* ==================================================================== */
856 /* Translate the units of measure. */
857 /* */
858 /* Note that even with a +units, or +to_meter in effect, it is */
859 /* still assumed that all the projection parameters are in */
860 /* meters. */
861 /* ==================================================================== */
862 if( psDefn->UOMLength == Linear_Meter )
863 {
864 strcpy( szUnits, "+units=m " );
865 }
866 else if( psDefn->UOMLength == Linear_Foot )
867 {
868 strcpy( szUnits, "+units=ft " );
869 }
870 else if( psDefn->UOMLength == Linear_Foot_US_Survey )
871 {
872 strcpy( szUnits, "+units=us-ft " );
873 }
874 else if( psDefn->UOMLength == Linear_Foot_Indian )
875 {
876 strcpy( szUnits, "+units=ind-ft " );
877 }
878 else if( psDefn->UOMLength == Linear_Link )
879 {
880 strcpy( szUnits, "+units=link " );
881 }
882 else if( psDefn->UOMLength == Linear_Yard_Indian)
883 {
884 strcpy( szUnits, "+units=ind-yd " );
885 }
886 else if( psDefn->UOMLength == Linear_Fathom )
887 {
888 strcpy( szUnits, "+units=fath " );
889 }
890 else if( psDefn->UOMLength == Linear_Mile_International_Nautical )
891 {
892 strcpy( szUnits, "+units=kmi " );
893 }
894 else
895 {
896 sprintf( szUnits, "+to_meter=%.10f", psDefn->UOMLengthInMeters );
897 }
898
899 /* -------------------------------------------------------------------- */
900 /* false easting and northing are in meters and that is what */
901 /* PROJ.4 wants regardless of the linear units. */
902 /* -------------------------------------------------------------------- */
903 dfFalseEasting = psDefn->ProjParm[5];
904 dfFalseNorthing = psDefn->ProjParm[6];
905
906 /* ==================================================================== */
907 /* Handle general projection methods. */
908 /* ==================================================================== */
909
910 /* -------------------------------------------------------------------- */
911 /* Geographic. */
912 /* -------------------------------------------------------------------- */
913 if(psDefn->Model==ModelTypeGeographic)
914 {
915 sprintf(szProjection+strlen(szProjection),"+proj=latlong ");
916
917 }
918
919 /* -------------------------------------------------------------------- */
920 /* UTM - special case override on transverse mercator so things */
921 /* will be more meaningful to the user. */
922 /* -------------------------------------------------------------------- */
923 else if( psDefn->MapSys == MapSys_UTM_North )
924 {
925 sprintf( szProjection+strlen(szProjection),
926 "+proj=utm +zone=%d ",
927 psDefn->Zone );
928 }
929
930 /* -------------------------------------------------------------------- */
931 /* Transverse Mercator */
932 /* -------------------------------------------------------------------- */
933 else if( psDefn->CTProjection == CT_TransverseMercator )
934 {
935 sprintf( szProjection+strlen(szProjection),
936 "+proj=tmerc +lat_0=%.9f +lon_0=%.9f +k=%f +x_0=%.3f +y_0=%.3f ",
937 psDefn->ProjParm[0],
938 psDefn->ProjParm[1],
939 psDefn->ProjParm[4],
940 dfFalseEasting,
941 dfFalseNorthing );
942 }
943
944 /* -------------------------------------------------------------------- */
945 /* Mercator */
946 /* -------------------------------------------------------------------- */
947 else if( psDefn->CTProjection == CT_Mercator )
948 {
949 if( psDefn->ProjParm[2] != 0.0 ) /* Mercator 2SP: FIXME we need a better way of detecting it */
950 sprintf( szProjection+strlen(szProjection),
951 "+proj=merc +lat_ts=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
952 psDefn->ProjParm[2],
953 psDefn->ProjParm[1],
954 dfFalseEasting,
955 dfFalseNorthing );
956 else
957 sprintf( szProjection+strlen(szProjection),
958 "+proj=merc +lat_ts=%.9f +lon_0=%.9f +k=%f +x_0=%.3f +y_0=%.3f ",
959 psDefn->ProjParm[0],
960 psDefn->ProjParm[1],
961 psDefn->ProjParm[4],
962 dfFalseEasting,
963 dfFalseNorthing );
964 }
965
966 /* -------------------------------------------------------------------- */
967 /* Cassini/Soldner */
968 /* -------------------------------------------------------------------- */
969 else if( psDefn->CTProjection == CT_CassiniSoldner )
970 {
971 sprintf( szProjection+strlen(szProjection),
972 "+proj=cass +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
973 psDefn->ProjParm[0],
974 psDefn->ProjParm[1],
975 dfFalseEasting,
976 dfFalseNorthing );
977 }
978
979 /* -------------------------------------------------------------------- */
980 /* Oblique Stereographic - Should this really map onto */
981 /* Stereographic? */
982 /* -------------------------------------------------------------------- */
983 else if( psDefn->CTProjection == CT_ObliqueStereographic )
984 {
985 sprintf( szProjection+strlen(szProjection),
986 "+proj=stere +lat_0=%.9f +lon_0=%.9f +k=%f +x_0=%.3f +y_0=%.3f ",
987 psDefn->ProjParm[0],
988 psDefn->ProjParm[1],
989 psDefn->ProjParm[4],
990 dfFalseEasting,
991 dfFalseNorthing );
992 }
993
994 /* -------------------------------------------------------------------- */
995 /* Stereographic */
996 /* -------------------------------------------------------------------- */
997 else if( psDefn->CTProjection == CT_Stereographic )
998 {
999 sprintf( szProjection+strlen(szProjection),
1000 "+proj=stere +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1001 psDefn->ProjParm[0],
1002 psDefn->ProjParm[1],
1003 dfFalseEasting,
1004 dfFalseNorthing );
1005 }
1006
1007 /* -------------------------------------------------------------------- */
1008 /* Polar Stereographic */
1009 /* -------------------------------------------------------------------- */
1010 else if( psDefn->CTProjection == CT_PolarStereographic )
1011 {
1012 if( psDefn->ProjParm[0] > 0.0 )
1013 sprintf( szProjection+strlen(szProjection),
1014 "+proj=stere +lat_0=90 +lat_ts=%.9f +lon_0=%.9f "
1015 "+k=%.9f +x_0=%.3f +y_0=%.3f ",
1016 psDefn->ProjParm[0],
1017 psDefn->ProjParm[1],
1018 psDefn->ProjParm[4],
1019 dfFalseEasting,
1020 dfFalseNorthing );
1021 else
1022 sprintf( szProjection+strlen(szProjection),
1023 "+proj=stere +lat_0=-90 +lat_ts=%.9f +lon_0=%.9f "
1024 "+k=%.9f +x_0=%.3f +y_0=%.3f ",
1025 psDefn->ProjParm[0],
1026 psDefn->ProjParm[1],
1027 psDefn->ProjParm[4],
1028 dfFalseEasting,
1029 dfFalseNorthing );
1030 }
1031
1032 /* -------------------------------------------------------------------- */
1033 /* Equirectangular */
1034 /* -------------------------------------------------------------------- */
1035 else if( psDefn->CTProjection == CT_Equirectangular )
1036 {
1037 sprintf( szProjection+strlen(szProjection),
1038 "+proj=eqc +lat_0=%.9f +lon_0=%.9f +lat_ts=%.9f +x_0=%.3f +y_0=%.3f ",
1039 psDefn->ProjParm[0],
1040 psDefn->ProjParm[1],
1041 psDefn->ProjParm[2],
1042 dfFalseEasting,
1043 dfFalseNorthing );
1044 }
1045
1046 /* -------------------------------------------------------------------- */
1047 /* Gnomonic */
1048 /* -------------------------------------------------------------------- */
1049 else if( psDefn->CTProjection == CT_Gnomonic )
1050 {
1051 sprintf( szProjection+strlen(szProjection),
1052 "+proj=gnom +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1053 psDefn->ProjParm[0],
1054 psDefn->ProjParm[1],
1055 dfFalseEasting,
1056 dfFalseNorthing );
1057 }
1058
1059 /* -------------------------------------------------------------------- */
1060 /* Orthographic */
1061 /* -------------------------------------------------------------------- */
1062 else if( psDefn->CTProjection == CT_Orthographic )
1063 {
1064 sprintf( szProjection+strlen(szProjection),
1065 "+proj=ortho +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1066 psDefn->ProjParm[0],
1067 psDefn->ProjParm[1],
1068 dfFalseEasting,
1069 dfFalseNorthing );
1070 }
1071
1072 /* -------------------------------------------------------------------- */
1073 /* Lambert Azimuthal Equal Area */
1074 /* -------------------------------------------------------------------- */
1075 else if( psDefn->CTProjection == CT_LambertAzimEqualArea )
1076 {
1077 sprintf( szProjection+strlen(szProjection),
1078 "+proj=laea +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1079 psDefn->ProjParm[0],
1080 psDefn->ProjParm[1],
1081 dfFalseEasting,
1082 dfFalseNorthing );
1083 }
1084
1085 /* -------------------------------------------------------------------- */
1086 /* Azimuthal Equidistant */
1087 /* -------------------------------------------------------------------- */
1088 else if( psDefn->CTProjection == CT_AzimuthalEquidistant )
1089 {
1090 sprintf( szProjection+strlen(szProjection),
1091 "+proj=aeqd +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1092 psDefn->ProjParm[0],
1093 psDefn->ProjParm[1],
1094 dfFalseEasting,
1095 dfFalseNorthing );
1096 }
1097
1098 /* -------------------------------------------------------------------- */
1099 /* Miller Cylindrical */
1100 /* -------------------------------------------------------------------- */
1101 else if( psDefn->CTProjection == CT_MillerCylindrical )
1102 {
1103 sprintf( szProjection+strlen(szProjection),
1104 "+proj=mill +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f +R_A ",
1105 psDefn->ProjParm[0],
1106 psDefn->ProjParm[1],
1107 dfFalseEasting,
1108 dfFalseNorthing );
1109 }
1110
1111 /* -------------------------------------------------------------------- */
1112 /* Polyconic */
1113 /* -------------------------------------------------------------------- */
1114 else if( psDefn->CTProjection == CT_Polyconic )
1115 {
1116 sprintf( szProjection+strlen(szProjection),
1117 "+proj=poly +lat_0=%.9f +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1118 psDefn->ProjParm[0],
1119 psDefn->ProjParm[1],
1120 dfFalseEasting,
1121 dfFalseNorthing );
1122 }
1123
1124 /* -------------------------------------------------------------------- */
1125 /* AlbersEqualArea */
1126 /* -------------------------------------------------------------------- */
1127 else if( psDefn->CTProjection == CT_AlbersEqualArea )
1128 {
1129 sprintf( szProjection+strlen(szProjection),
1130 "+proj=aea +lat_1=%.9f +lat_2=%.9f +lat_0=%.9f +lon_0=%.9f"
1131 " +x_0=%.3f +y_0=%.3f ",
1132 psDefn->ProjParm[0],
1133 psDefn->ProjParm[1],
1134 psDefn->ProjParm[2],
1135 psDefn->ProjParm[3],
1136 dfFalseEasting,
1137 dfFalseNorthing );
1138 }
1139
1140 /* -------------------------------------------------------------------- */
1141 /* EquidistantConic */
1142 /* -------------------------------------------------------------------- */
1143 else if( psDefn->CTProjection == CT_EquidistantConic )
1144 {
1145 sprintf( szProjection+strlen(szProjection),
1146 "+proj=eqdc +lat_1=%.9f +lat_2=%.9f +lat_0=%.9f +lon_0=%.9f"
1147 " +x_0=%.3f +y_0=%.3f ",
1148 psDefn->ProjParm[0],
1149 psDefn->ProjParm[1],
1150 psDefn->ProjParm[2],
1151 psDefn->ProjParm[3],
1152 dfFalseEasting,
1153 dfFalseNorthing );
1154 }
1155
1156 /* -------------------------------------------------------------------- */
1157 /* Robinson */
1158 /* -------------------------------------------------------------------- */
1159 else if( psDefn->CTProjection == CT_Robinson )
1160 {
1161 sprintf( szProjection+strlen(szProjection),
1162 "+proj=robin +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1163 psDefn->ProjParm[1],
1164 dfFalseEasting,
1165 dfFalseNorthing );
1166 }
1167
1168 /* -------------------------------------------------------------------- */
1169 /* VanDerGrinten */
1170 /* -------------------------------------------------------------------- */
1171 else if( psDefn->CTProjection == CT_VanDerGrinten )
1172 {
1173 sprintf( szProjection+strlen(szProjection),
1174 "+proj=vandg +lon_0=%.9f +x_0=%.3f +y_0=%.3f +R_A ",
1175 psDefn->ProjParm[1],
1176 dfFalseEasting,
1177 dfFalseNorthing );
1178 }
1179
1180 /* -------------------------------------------------------------------- */
1181 /* Sinusoidal */
1182 /* -------------------------------------------------------------------- */
1183 else if( psDefn->CTProjection == CT_Sinusoidal )
1184 {
1185 sprintf( szProjection+strlen(szProjection),
1186 "+proj=sinu +lon_0=%.9f +x_0=%.3f +y_0=%.3f ",
1187 psDefn->ProjParm[1],
1188 dfFalseEasting,
1189 dfFalseNorthing );
1190 }
1191
1192 /* -------------------------------------------------------------------- */
1193 /* LambertConfConic_2SP */
1194 /* -------------------------------------------------------------------- */
1195 else if( psDefn->CTProjection == CT_LambertConfConic_2SP )
1196 {
1197 sprintf( szProjection+strlen(szProjection),
1198 "+proj=lcc +lat_0=%.9f +lon_0=%.9f +lat_1=%.9f +lat_2=%.9f "
1199 " +x_0=%.3f +y_0=%.3f ",
1200 psDefn->ProjParm[0],
1201 psDefn->ProjParm[1],
1202 psDefn->ProjParm[2],
1203 psDefn->ProjParm[3],
1204 dfFalseEasting,
1205 dfFalseNorthing );
1206 }
1207
1208 /* -------------------------------------------------------------------- */
1209 /* LambertConfConic_1SP */
1210 /* -------------------------------------------------------------------- */
1211 else if( psDefn->CTProjection == CT_LambertConfConic_1SP )
1212 {
1213 sprintf( szProjection+strlen(szProjection),
1214 "+proj=lcc +lat_0=%.9f +lat_1=%.9f +lon_0=%.9f"
1215 " +k_0=%.9f +x_0=%.3f +y_0=%.3f ",
1216 psDefn->ProjParm[0],
1217 psDefn->ProjParm[0],
1218 psDefn->ProjParm[1],
1219 psDefn->ProjParm[4],
1220 psDefn->ProjParm[5],
1221 psDefn->ProjParm[6] );
1222 }
1223
1224 /* -------------------------------------------------------------------- */
1225 /* CT_CylindricalEqualArea */
1226 /* -------------------------------------------------------------------- */
1227 else if( psDefn->CTProjection == CT_CylindricalEqualArea )
1228 {
1229 sprintf( szProjection+strlen(szProjection),
1230 "+proj=cea +lat_ts=%.9f +lon_0=%.9f "
1231 " +x_0=%.3f +y_0=%.3f ",
1232 psDefn->ProjParm[0],
1233 psDefn->ProjParm[1],
1234 psDefn->ProjParm[5],
1235 psDefn->ProjParm[6] );
1236 }
1237
1238 /* -------------------------------------------------------------------- */
1239 /* NewZealandMapGrid */
1240 /* -------------------------------------------------------------------- */
1241 else if( psDefn->CTProjection == CT_NewZealandMapGrid )
1242 {
1243 sprintf( szProjection+strlen(szProjection),
1244 "+proj=nzmg +lat_0=%.9f +lon_0=%.9f"
1245 " +x_0=%.3f +y_0=%.3f ",
1246 psDefn->ProjParm[0],
1247 psDefn->ProjParm[1],
1248 psDefn->ProjParm[5],
1249 psDefn->ProjParm[6] );
1250 }
1251
1252 /* -------------------------------------------------------------------- */
1253 /* Transverse Mercator - south oriented. */
1254 /* -------------------------------------------------------------------- */
1255 else if( psDefn->CTProjection == CT_TransvMercator_SouthOriented )
1256 {
1257 /* this appears to be an unsupported formulation with PROJ.4 */
1258 }
1259
1260 /* -------------------------------------------------------------------- */
1261 /* ObliqueMercator (Hotine) */
1262 /* -------------------------------------------------------------------- */
1263 else if( psDefn->CTProjection == CT_ObliqueMercator )
1264 {
1265 /* not clear how ProjParm[3] - angle from rectified to skewed grid -
1266 should be applied ... see the +not_rot flag for PROJ.4.
1267 Just ignoring for now. */
1268
1269 sprintf( szProjection+strlen(szProjection),
1270 "+proj=omerc +lat_0=%.9f +lonc=%.9f +alpha=%.9f"
1271 " +k=%.9f +x_0=%.3f +y_0=%.3f ",
1272 psDefn->ProjParm[0],
1273 psDefn->ProjParm[1],
1274 psDefn->ProjParm[2],
1275 psDefn->ProjParm[4],
1276 psDefn->ProjParm[5],
1277 psDefn->ProjParm[6] );
1278 }
1279
1280 /* ==================================================================== */
1281 /* Handle ellipsoid information. */
1282 /* ==================================================================== */
1283 if( psDefn->Ellipsoid == Ellipse_WGS_84 )
1284 strcat( szProjection, "+ellps=WGS84 " );
1285 else if( psDefn->Ellipsoid == Ellipse_Clarke_1866 )
1286 strcat( szProjection, "+ellps=clrk66 " );
1287 else if( psDefn->Ellipsoid == Ellipse_Clarke_1880 )
1288 strcat( szProjection, "+ellps=clrk80 " );
1289 else if( psDefn->Ellipsoid == Ellipse_GRS_1980 )
1290 strcat( szProjection, "+ellps=GRS80 " );
1291 else
1292 {
1293 if( psDefn->SemiMajor != 0.0 && psDefn->SemiMinor != 0.0 )
1294 {
1295 sprintf( szProjection+strlen(szProjection),
1296 "+a=%.3f +b=%.3f ",
1297 psDefn->SemiMajor,
1298 psDefn->SemiMinor );
1299 }
1300 }
1301
1302 strcat( szProjection, szUnits );
1303
1304 /* If we don't have anything, reset */
1305 if (strstr(szProjection, "+proj=") == NULL) { return CPLStrdup(""); }
1306
1307 return( CPLStrdup( szProjection ) );
1308 }
1309
1310 #if !defined(HAVE_LIBPROJ)
1311
GTIFProj4ToLatLong(GTIFDefn * psDefn,int nPoints,double * padfX,double * padfY)1312 int GTIFProj4ToLatLong( GTIFDefn * psDefn, int nPoints,
1313 double *padfX, double *padfY )
1314 {
1315 (void) psDefn;
1316 (void) nPoints;
1317 (void) padfX;
1318 (void) padfY;
1319 #ifdef DEBUG
1320 fprintf( stderr,
1321 "GTIFProj4ToLatLong() - PROJ.4 support not compiled in.\n" );
1322 #endif
1323 return FALSE;
1324 }
1325
GTIFProj4FromLatLong(GTIFDefn * psDefn,int nPoints,double * padfX,double * padfY)1326 int GTIFProj4FromLatLong( GTIFDefn * psDefn, int nPoints,
1327 double *padfX, double *padfY )
1328 {
1329 (void) psDefn;
1330 (void) nPoints;
1331 (void) padfX;
1332 (void) padfY;
1333 #ifdef DEBUG
1334 fprintf( stderr,
1335 "GTIFProj4FromLatLong() - PROJ.4 support not compiled in.\n" );
1336 #endif
1337 return FALSE;
1338 }
1339 #else
1340
1341 #include "proj_api.h"
1342
1343 /************************************************************************/
1344 /* GTIFProj4FromLatLong() */
1345 /* */
1346 /* Convert lat/long values to projected coordinate for a */
1347 /* particular definition. */
1348 /************************************************************************/
1349
GTIFProj4FromLatLong(GTIFDefn * psDefn,int nPoints,double * padfX,double * padfY)1350 int GTIFProj4FromLatLong( GTIFDefn * psDefn, int nPoints,
1351 double *padfX, double *padfY )
1352
1353 {
1354 char *pszProjection, **papszArgs;
1355 projPJ *psPJ;
1356 int i;
1357
1358 /* -------------------------------------------------------------------- */
1359 /* Get a projection definition. */
1360 /* -------------------------------------------------------------------- */
1361 pszProjection = GTIFGetProj4Defn( psDefn );
1362
1363 if( pszProjection == NULL )
1364 return FALSE;
1365
1366 /* -------------------------------------------------------------------- */
1367 /* Parse into tokens for pj_init(), and initialize the projection. */
1368 /* -------------------------------------------------------------------- */
1369
1370 papszArgs = CSLTokenizeStringComplex( pszProjection, " +", TRUE, FALSE );
1371 free( pszProjection );
1372
1373 psPJ = pj_init( CSLCount(papszArgs), papszArgs );
1374
1375 CSLDestroy( papszArgs );
1376
1377 if( psPJ == NULL )
1378 {
1379 return FALSE;
1380 }
1381
1382 /* -------------------------------------------------------------------- */
1383 /* Process each of the points. */
1384 /* -------------------------------------------------------------------- */
1385 for( i = 0; i < nPoints; i++ )
1386 {
1387 projUV sUV;
1388
1389 sUV.u = padfX[i] * DEG_TO_RAD;
1390 sUV.v = padfY[i] * DEG_TO_RAD;
1391
1392 sUV = pj_fwd( sUV, psPJ );
1393
1394 padfX[i] = sUV.u;
1395 padfY[i] = sUV.v;
1396 }
1397
1398 pj_free( psPJ );
1399
1400 return TRUE;
1401 }
1402
1403 /************************************************************************/
1404 /* GTIFProj4ToLatLong() */
1405 /* */
1406 /* Convert projection coordinates to lat/long for a particular */
1407 /* definition. */
1408 /************************************************************************/
1409
GTIFProj4ToLatLong(GTIFDefn * psDefn,int nPoints,double * padfX,double * padfY)1410 int GTIFProj4ToLatLong( GTIFDefn * psDefn, int nPoints,
1411 double *padfX, double *padfY )
1412
1413 {
1414 char *pszProjection, **papszArgs;
1415 projPJ *psPJ;
1416 int i;
1417
1418 /* -------------------------------------------------------------------- */
1419 /* Get a projection definition. */
1420 /* -------------------------------------------------------------------- */
1421 pszProjection = GTIFGetProj4Defn( psDefn );
1422
1423 if( pszProjection == NULL )
1424 return FALSE;
1425
1426 /* -------------------------------------------------------------------- */
1427 /* Parse into tokens for pj_init(), and initialize the projection. */
1428 /* -------------------------------------------------------------------- */
1429
1430 papszArgs = CSLTokenizeStringComplex( pszProjection, " +", TRUE, FALSE );
1431 free( pszProjection );
1432
1433 psPJ = pj_init( CSLCount(papszArgs), papszArgs );
1434
1435 CSLDestroy( papszArgs );
1436
1437 if( psPJ == NULL )
1438 {
1439 return FALSE;
1440 }
1441
1442 /* -------------------------------------------------------------------- */
1443 /* Process each of the points. */
1444 /* -------------------------------------------------------------------- */
1445 for( i = 0; i < nPoints; i++ )
1446 {
1447 projUV sUV;
1448
1449 sUV.u = padfX[i];
1450 sUV.v = padfY[i];
1451
1452 sUV = pj_inv( sUV, psPJ );
1453
1454 padfX[i] = sUV.u * RAD_TO_DEG;
1455 padfY[i] = sUV.v * RAD_TO_DEG;
1456 }
1457
1458 pj_free( psPJ );
1459
1460 return TRUE;
1461 }
1462
1463 #endif /* has proj_api.h and -lproj */
1464
1465