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