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