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