1 : /******************************************************************************
2 : * $Id: mffdataset.cpp 18182 2009-12-05 01:12:13Z warmerdam $
3 : *
4 : * Project: GView
5 : * Purpose: Implementation of Atlantis MFF Support
6 : * Author: Frank Warmerdam, warmerda@home.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2000, Frank Warmerdam
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "cpl_string.h"
32 : #include <ctype.h>
33 : #include "ogr_spatialref.h"
34 : #include "atlsci_spheroid.h"
35 :
36 : CPL_CVSID("$Id: mffdataset.cpp 18182 2009-12-05 01:12:13Z warmerdam $");
37 :
38 : CPL_C_START
39 : void GDALRegister_MFF(void);
40 : CPL_C_END
41 :
42 : enum {
43 : MFFPRJ_NONE,
44 : MFFPRJ_LL,
45 : MFFPRJ_UTM,
46 : MFFPRJ_UNRECOGNIZED
47 : } ;
48 :
49 : static int GetMFFProjectionType(const char * pszNewProjection);
50 :
51 : /************************************************************************/
52 : /* ==================================================================== */
53 : /* MFFDataset */
54 : /* ==================================================================== */
55 : /************************************************************************/
56 :
57 : class MFFDataset : public RawDataset
58 : {
59 : FILE *fpImage; // image data file.
60 :
61 : int nGCPCount;
62 : GDAL_GCP *pasGCPList;
63 :
64 : char *pszProjection;
65 : char *pszGCPProjection;
66 : double adfGeoTransform[6];
67 :
68 : void ScanForGCPs();
69 : void ScanForProjectionInfo();
70 :
71 :
72 : public:
73 : MFFDataset();
74 : ~MFFDataset();
75 :
76 : char **papszHdrLines;
77 :
78 : FILE **pafpBandFiles;
79 :
80 : virtual int GetGCPCount();
81 : virtual const char *GetGCPProjection();
82 : virtual const GDAL_GCP *GetGCPs();
83 :
84 : virtual const char *GetProjectionRef();
85 : virtual CPLErr GetGeoTransform( double * );
86 :
87 : static GDALDataset *Open( GDALOpenInfo * );
88 : static GDALDataset *Create( const char * pszFilename,
89 : int nXSize, int nYSize, int nBands,
90 : GDALDataType eType, char ** papszParmList );
91 : static GDALDataset *CreateCopy( const char * pszFilename,
92 : GDALDataset *poSrcDS,
93 : int bStrict, char ** papszOptions,
94 : GDALProgressFunc pfnProgress,
95 : void * pProgressData );
96 :
97 : };
98 :
99 : /************************************************************************/
100 : /* ==================================================================== */
101 : /* MFFTiledBand */
102 : /* ==================================================================== */
103 : /************************************************************************/
104 :
105 : class MFFTiledBand : public GDALRasterBand
106 : {
107 : friend class MFFDataset;
108 :
109 : FILE *fpRaw;
110 : int bNative;
111 :
112 : public:
113 :
114 : MFFTiledBand( MFFDataset *, int, FILE *, int, int,
115 : GDALDataType, int );
116 : ~MFFTiledBand();
117 :
118 : virtual CPLErr IReadBlock( int, int, void * );
119 : };
120 :
121 :
122 : /************************************************************************/
123 : /* MFFTiledBand() */
124 : /************************************************************************/
125 :
126 : MFFTiledBand::MFFTiledBand( MFFDataset *poDS, int nBand, FILE *fp,
127 : int nTileXSize, int nTileYSize,
128 1 : GDALDataType eDataType, int bNative )
129 :
130 : {
131 1 : this->poDS = poDS;
132 1 : this->nBand = nBand;
133 :
134 1 : this->eDataType = eDataType;
135 :
136 1 : this->bNative = bNative;
137 :
138 1 : this->nBlockXSize = nTileXSize;
139 1 : this->nBlockYSize = nTileYSize;
140 :
141 1 : this->fpRaw = fp;
142 1 : }
143 :
144 : /************************************************************************/
145 : /* ~MFFTiledBand() */
146 : /************************************************************************/
147 :
148 1 : MFFTiledBand::~MFFTiledBand()
149 :
150 : {
151 1 : VSIFCloseL( fpRaw );
152 1 : }
153 :
154 :
155 : /************************************************************************/
156 : /* IReadBlock() */
157 : /************************************************************************/
158 :
159 : CPLErr MFFTiledBand::IReadBlock( int nBlockXOff, int nBlockYOff,
160 1 : void * pImage )
161 :
162 : {
163 : long nOffset;
164 : int nTilesPerRow;
165 : int nWordSize, nBlockSize;
166 :
167 1 : nTilesPerRow = (nRasterXSize + nBlockXSize - 1) / nBlockXSize;
168 1 : nWordSize = GDALGetDataTypeSize( eDataType ) / 8;
169 1 : nBlockSize = nWordSize * nBlockXSize * nBlockYSize;
170 :
171 1 : nOffset = nBlockSize * (nBlockXOff + nBlockYOff*nTilesPerRow);
172 :
173 1 : if( VSIFSeekL( fpRaw, nOffset, SEEK_SET ) == -1
174 : || VSIFReadL( pImage, 1, nBlockSize, fpRaw ) < 1 )
175 : {
176 : CPLError( CE_Failure, CPLE_FileIO,
177 : "Read of tile %d/%d failed with fseek or fread error.",
178 0 : nBlockXOff, nBlockYOff );
179 0 : return CE_Failure;
180 : }
181 :
182 1 : if( !bNative && nWordSize > 1 )
183 : {
184 0 : if( GDALDataTypeIsComplex( eDataType ) )
185 : {
186 : GDALSwapWords( pImage, nWordSize/2, nBlockXSize*nBlockYSize,
187 0 : nWordSize );
188 : GDALSwapWords( ((GByte *) pImage)+nWordSize/2,
189 0 : nWordSize/2, nBlockXSize*nBlockYSize, nWordSize );
190 : }
191 : else
192 : GDALSwapWords( pImage, nWordSize,
193 0 : nBlockXSize * nBlockYSize, nWordSize );
194 : }
195 :
196 1 : return CE_None;
197 : }
198 :
199 : /************************************************************************/
200 : /* MFF Spheroids */
201 : /************************************************************************/
202 :
203 : class MFFSpheroidList : public SpheroidList
204 : {
205 :
206 : public:
207 :
208 : MFFSpheroidList();
209 : ~MFFSpheroidList();
210 :
211 : };
212 :
213 12 : MFFSpheroidList :: MFFSpheroidList()
214 : {
215 12 : num_spheroids = 18;
216 :
217 12 : epsilonR = 0.1;
218 12 : epsilonI = 0.000001;
219 :
220 12 : spheroids[0].SetValuesByRadii("SPHERE",6371007.0,6371007.0);
221 12 : spheroids[1].SetValuesByRadii("EVEREST",6377304.0,6356103.0);
222 12 : spheroids[2].SetValuesByRadii("BESSEL",6377397.0,6356082.0);
223 12 : spheroids[3].SetValuesByRadii("AIRY",6377563.0,6356300.0);
224 12 : spheroids[4].SetValuesByRadii("CLARKE_1858",6378294.0,6356621.0);
225 12 : spheroids[5].SetValuesByRadii("CLARKE_1866",6378206.4,6356583.8);
226 12 : spheroids[6].SetValuesByRadii("CLARKE_1880",6378249.0,6356517.0);
227 12 : spheroids[7].SetValuesByRadii("HAYFORD",6378388.0,6356915.0);
228 12 : spheroids[8].SetValuesByRadii("KRASOVSKI",6378245.0,6356863.0);
229 12 : spheroids[9].SetValuesByRadii("HOUGH",6378270.0,6356794.0);
230 12 : spheroids[10].SetValuesByRadii("FISHER_60",6378166.0,6356784.0);
231 12 : spheroids[11].SetValuesByRadii("KAULA",6378165.0,6356345.0);
232 12 : spheroids[12].SetValuesByRadii("IUGG_67",6378160.0,6356775.0);
233 12 : spheroids[13].SetValuesByRadii("FISHER_68",6378150.0,6356330.0);
234 12 : spheroids[14].SetValuesByRadii("WGS_72",6378135.0,6356751.0);
235 12 : spheroids[15].SetValuesByRadii("IUGG_75",6378140.0,6356755.0);
236 12 : spheroids[16].SetValuesByRadii("WGS_84",6378137.0,6356752.0);
237 12 : spheroids[17].SetValuesByRadii("HUGHES",6378273.0,6356889.4);
238 12 : }
239 :
240 12 : MFFSpheroidList::~MFFSpheroidList()
241 :
242 : {
243 12 : }
244 :
245 : /************************************************************************/
246 : /* ==================================================================== */
247 : /* MFFDataset */
248 : /* ==================================================================== */
249 : /************************************************************************/
250 :
251 : /************************************************************************/
252 : /* MFFDataset() */
253 : /************************************************************************/
254 :
255 26 : MFFDataset::MFFDataset()
256 : {
257 26 : papszHdrLines = NULL;
258 26 : pafpBandFiles = NULL;
259 26 : nGCPCount = 0;
260 26 : pasGCPList = NULL;
261 :
262 26 : pszProjection = CPLStrdup("");
263 26 : pszGCPProjection = CPLStrdup("");
264 26 : adfGeoTransform[0] = 0.0;
265 26 : adfGeoTransform[1] = 1.0;
266 26 : adfGeoTransform[2] = 0.0;
267 26 : adfGeoTransform[3] = 0.0;
268 26 : adfGeoTransform[4] = 0.0;
269 26 : adfGeoTransform[5] = 1.0;
270 26 : }
271 :
272 : /************************************************************************/
273 : /* ~MFFDataset() */
274 : /************************************************************************/
275 :
276 26 : MFFDataset::~MFFDataset()
277 :
278 : {
279 26 : FlushCache();
280 26 : CSLDestroy( papszHdrLines );
281 26 : if( pafpBandFiles != NULL )
282 : {
283 0 : for( int i = 0; i < GetRasterCount(); i++ )
284 : {
285 0 : if( pafpBandFiles[i] != NULL )
286 0 : VSIFCloseL( pafpBandFiles[i] );
287 : }
288 0 : CPLFree( pafpBandFiles );
289 : }
290 :
291 26 : if( nGCPCount > 0 )
292 : {
293 3 : GDALDeinitGCPs( nGCPCount, pasGCPList );
294 : }
295 26 : CPLFree( pasGCPList );
296 26 : CPLFree( pszProjection );
297 26 : CPLFree( pszGCPProjection );
298 :
299 26 : }
300 :
301 : /************************************************************************/
302 : /* GetGCPCount() */
303 : /************************************************************************/
304 :
305 0 : int MFFDataset::GetGCPCount()
306 :
307 : {
308 0 : return nGCPCount;
309 : }
310 :
311 : /************************************************************************/
312 : /* GetGCPProjection() */
313 : /************************************************************************/
314 :
315 0 : const char *MFFDataset::GetGCPProjection()
316 :
317 : {
318 0 : if( nGCPCount > 0 )
319 0 : return pszGCPProjection;
320 : else
321 0 : return "";
322 : }
323 :
324 : /************************************************************************/
325 : /* GetProjectionRef() */
326 : /************************************************************************/
327 :
328 18 : const char *MFFDataset::GetProjectionRef()
329 :
330 : {
331 18 : return ( pszProjection );
332 : }
333 :
334 : /************************************************************************/
335 : /* GetGeoTransform() */
336 : /************************************************************************/
337 :
338 9 : CPLErr MFFDataset::GetGeoTransform( double * padfTransform )
339 :
340 : {
341 9 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
342 9 : return( CE_None );
343 : }
344 :
345 : /************************************************************************/
346 : /* GetGCP() */
347 : /************************************************************************/
348 :
349 0 : const GDAL_GCP *MFFDataset::GetGCPs()
350 :
351 : {
352 0 : return pasGCPList;
353 : }
354 :
355 : /************************************************************************/
356 : /* ScanForGCPs() */
357 : /************************************************************************/
358 :
359 26 : void MFFDataset::ScanForGCPs()
360 :
361 : {
362 : int nCorner;
363 26 : int NUM_GCPS = 0;
364 :
365 26 : if( CSLFetchNameValue(papszHdrLines, "NUM_GCPS") != NULL )
366 2 : NUM_GCPS = atoi(CSLFetchNameValue(papszHdrLines, "NUM_GCPS"));
367 26 : if (NUM_GCPS < 0)
368 0 : return;
369 :
370 26 : nGCPCount = 0;
371 26 : pasGCPList = (GDAL_GCP *) VSICalloc(sizeof(GDAL_GCP),5+NUM_GCPS);
372 26 : if (pasGCPList == NULL)
373 0 : return;
374 :
375 156 : for( nCorner = 0; nCorner < 5; nCorner++ )
376 : {
377 130 : const char * pszBase=NULL;
378 130 : double dfRasterX=0.0, dfRasterY=0.0;
379 : char szLatName[40], szLongName[40];
380 :
381 130 : if( nCorner == 0 )
382 : {
383 26 : dfRasterX = 0.5;
384 26 : dfRasterY = 0.5;
385 26 : pszBase = "TOP_LEFT_CORNER";
386 : }
387 104 : else if( nCorner == 1 )
388 : {
389 26 : dfRasterX = GetRasterXSize()-0.5;
390 26 : dfRasterY = 0.5;
391 26 : pszBase = "TOP_RIGHT_CORNER";
392 : }
393 78 : else if( nCorner == 2 )
394 : {
395 26 : dfRasterX = GetRasterXSize()-0.5;
396 26 : dfRasterY = GetRasterYSize()-0.5;
397 26 : pszBase = "BOTTOM_RIGHT_CORNER";
398 : }
399 52 : else if( nCorner == 3 )
400 : {
401 26 : dfRasterX = 0.5;
402 26 : dfRasterY = GetRasterYSize()-0.5;
403 26 : pszBase = "BOTTOM_LEFT_CORNER";
404 : }
405 26 : else if( nCorner == 4 )
406 : {
407 26 : dfRasterX = GetRasterXSize()/2.0;
408 26 : dfRasterY = GetRasterYSize()/2.0;
409 26 : pszBase = "CENTRE";
410 : }
411 :
412 130 : sprintf( szLatName, "%s_LATITUDE", pszBase );
413 130 : sprintf( szLongName, "%s_LONGITUDE", pszBase );
414 :
415 130 : if( CSLFetchNameValue(papszHdrLines, szLatName) != NULL
416 : && CSLFetchNameValue(papszHdrLines, szLongName) != NULL )
417 : {
418 5 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
419 :
420 5 : CPLFree( pasGCPList[nGCPCount].pszId );
421 :
422 5 : pasGCPList[nGCPCount].pszId = CPLStrdup( pszBase );
423 :
424 : pasGCPList[nGCPCount].dfGCPX =
425 5 : atof(CSLFetchNameValue(papszHdrLines, szLongName));
426 : pasGCPList[nGCPCount].dfGCPY =
427 5 : atof(CSLFetchNameValue(papszHdrLines, szLatName));
428 5 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
429 :
430 5 : pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
431 5 : pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
432 :
433 5 : nGCPCount++;
434 : }
435 :
436 : }
437 :
438 : /* -------------------------------------------------------------------- */
439 : /* Collect standalone GCPs. They look like: */
440 : /* */
441 : /* GCPn = row, col, lat, long */
442 : /* GCP1 = 1, 1, 45.0, -75.0 */
443 : /* -------------------------------------------------------------------- */
444 : int i;
445 :
446 28 : for( i = 0; i < NUM_GCPS; i++ )
447 : {
448 : char szName[25];
449 : char **papszTokens;
450 :
451 2 : sprintf( szName, "GCP%d", i+1 );
452 2 : if( CSLFetchNameValue( papszHdrLines, szName ) == NULL )
453 0 : continue;
454 :
455 : papszTokens = CSLTokenizeStringComplex(
456 : CSLFetchNameValue( papszHdrLines, szName ),
457 2 : ",", FALSE, FALSE );
458 2 : if( CSLCount(papszTokens) == 4 )
459 : {
460 2 : GDALInitGCPs( 1, pasGCPList + nGCPCount );
461 :
462 2 : CPLFree( pasGCPList[nGCPCount].pszId );
463 2 : pasGCPList[nGCPCount].pszId = CPLStrdup( szName );
464 :
465 2 : pasGCPList[nGCPCount].dfGCPX = atof(papszTokens[3]);
466 2 : pasGCPList[nGCPCount].dfGCPY = atof(papszTokens[2]);
467 2 : pasGCPList[nGCPCount].dfGCPZ = 0.0;
468 2 : pasGCPList[nGCPCount].dfGCPPixel = atof(papszTokens[1])+0.5;
469 2 : pasGCPList[nGCPCount].dfGCPLine = atof(papszTokens[0])+0.5;
470 :
471 2 : nGCPCount++;
472 : }
473 :
474 2 : CSLDestroy(papszTokens);
475 : }
476 : }
477 :
478 : /************************************************************************/
479 : /* ScanForProjectionInfo */
480 : /************************************************************************/
481 :
482 26 : void MFFDataset::ScanForProjectionInfo()
483 : {
484 : const char *pszProjName, *pszOriginLong, *pszSpheroidName;
485 : const char *pszSpheroidEqRadius, *pszSpheroidPolarRadius;
486 : double eq_radius, polar_radius;
487 26 : OGRSpatialReference oProj;
488 26 : OGRSpatialReference oLL;
489 : MFFSpheroidList *mffEllipsoids;
490 :
491 : pszProjName = CSLFetchNameValue(papszHdrLines,
492 26 : "PROJECTION_NAME");
493 : pszOriginLong = CSLFetchNameValue(papszHdrLines,
494 26 : "PROJECTION_ORIGIN_LONGITUDE");
495 : pszSpheroidName = CSLFetchNameValue(papszHdrLines,
496 26 : "SPHEROID_NAME");
497 :
498 26 : if (pszProjName == NULL)
499 : {
500 23 : CPLFree( pszProjection );
501 23 : CPLFree( pszGCPProjection );
502 23 : pszProjection=CPLStrdup("");
503 23 : pszGCPProjection=CPLStrdup("");
504 23 : return;
505 : }
506 3 : else if ((!EQUAL(pszProjName,"utm")) && (!EQUAL(pszProjName,"ll")))
507 : {
508 : CPLError(CE_Warning,CPLE_AppDefined,
509 0 : "Warning- only utm and lat/long projections are currently supported.");
510 0 : CPLFree( pszProjection );
511 0 : CPLFree( pszGCPProjection );
512 0 : pszProjection=CPLStrdup("");
513 0 : pszGCPProjection=CPLStrdup("");
514 : return;
515 : }
516 3 : mffEllipsoids = new MFFSpheroidList;
517 :
518 3 : if( EQUAL(pszProjName,"utm") )
519 : {
520 : int nZone;
521 :
522 3 : if (pszOriginLong == NULL)
523 : {
524 : /* If origin not specified, assume 0.0 */
525 : CPLError(CE_Warning,CPLE_AppDefined,
526 2 : "Warning- no projection origin longitude specified. Assuming 0.0.");
527 2 : nZone = 31;
528 : }
529 : else
530 1 : nZone = 31 + (int) floor(atof(pszOriginLong)/6.0);
531 :
532 :
533 3 : if( nGCPCount >= 5 && pasGCPList[4].dfGCPY < 0 )
534 0 : oProj.SetUTM( nZone, 0 );
535 : else
536 3 : oProj.SetUTM( nZone, 1 );
537 :
538 3 : if (pszOriginLong != NULL)
539 1 : oProj.SetProjParm(SRS_PP_CENTRAL_MERIDIAN,atof(pszOriginLong));
540 :
541 : }
542 :
543 3 : if (pszOriginLong != NULL)
544 1 : oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN,atof(pszOriginLong));
545 :
546 3 : if ((pszSpheroidName == NULL))
547 : {
548 : CPLError(CE_Warning,CPLE_AppDefined,
549 2 : "Warning- unspecified ellipsoid. Using wgs-84 parameters.\n");
550 :
551 2 : oProj.SetWellKnownGeogCS( "WGS84" );
552 2 : oLL.SetWellKnownGeogCS( "WGS84" );
553 : }
554 : else
555 : {
556 1 : if (mffEllipsoids->SpheroidInList(pszSpheroidName))
557 : {
558 : oProj.SetGeogCS( "unknown","unknown",pszSpheroidName,
559 : mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
560 : mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
561 1 : );
562 : oLL.SetGeogCS( "unknown","unknown",pszSpheroidName,
563 : mffEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
564 : mffEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName)
565 1 : );
566 : }
567 0 : else if (EQUAL(pszSpheroidName,"USER_DEFINED"))
568 : {
569 : pszSpheroidEqRadius = CSLFetchNameValue(papszHdrLines,
570 0 : "SPHEROID_EQUATORIAL_RADIUS");
571 : pszSpheroidPolarRadius = CSLFetchNameValue(papszHdrLines,
572 0 : "SPHEROID_POLAR_RADIUS");
573 0 : if ((pszSpheroidEqRadius != NULL) && (pszSpheroidPolarRadius != NULL))
574 : {
575 0 : eq_radius = atof( pszSpheroidEqRadius );
576 0 : polar_radius = atof( pszSpheroidPolarRadius );
577 : oProj.SetGeogCS( "unknown","unknown","unknown",
578 0 : eq_radius, eq_radius/(eq_radius - polar_radius));
579 : oLL.SetGeogCS( "unknown","unknown","unknown",
580 0 : eq_radius, eq_radius/(eq_radius - polar_radius));
581 : }
582 : else
583 : {
584 : CPLError(CE_Warning,CPLE_AppDefined,
585 0 : "Warning- radii not specified for user-defined ellipsoid. Using wgs-84 parameters. \n");
586 0 : oProj.SetWellKnownGeogCS( "WGS84" );
587 0 : oLL.SetWellKnownGeogCS( "WGS84" );
588 : }
589 : }
590 : else
591 : {
592 : CPLError(CE_Warning,CPLE_AppDefined,
593 0 : "Warning- unrecognized ellipsoid. Using wgs-84 parameters.\n");
594 0 : oProj.SetWellKnownGeogCS( "WGS84" );
595 0 : oLL.SetWellKnownGeogCS( "WGS84" );
596 : }
597 : }
598 :
599 : /* If a geotransform is sufficient to represent the GCP's (ie. each */
600 : /* estimated gcp is within 0.25*pixel size of the actual value- this */
601 : /* is the test applied by GDALGCPsToGeoTransform), store the */
602 : /* geotransform. */
603 3 : int transform_ok = FALSE;
604 :
605 3 : if (EQUAL(pszProjName,"LL"))
606 : {
607 0 : transform_ok = GDALGCPsToGeoTransform(nGCPCount,pasGCPList,adfGeoTransform,0);
608 : }
609 : else
610 : {
611 3 : OGRCoordinateTransformation *poTransform = NULL;
612 : double *dfPrjX, *dfPrjY;
613 : int gcp_index;
614 3 : int bSuccess = TRUE;
615 :
616 3 : dfPrjX = (double *) CPLMalloc(nGCPCount*sizeof(double));
617 3 : dfPrjY = (double *) CPLMalloc(nGCPCount*sizeof(double));
618 :
619 :
620 3 : poTransform = OGRCreateCoordinateTransformation( &oLL, &oProj );
621 3 : if( poTransform == NULL )
622 : {
623 0 : CPLErrorReset();
624 0 : bSuccess = FALSE;
625 : }
626 :
627 10 : for(gcp_index=0;gcp_index<nGCPCount;gcp_index++)
628 : {
629 7 : dfPrjX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
630 7 : dfPrjY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
631 :
632 7 : if( bSuccess && !poTransform->Transform( 1, &(dfPrjX[gcp_index]), &(dfPrjY[gcp_index]) ) )
633 0 : bSuccess = FALSE;
634 :
635 : }
636 :
637 3 : if( bSuccess )
638 : {
639 :
640 10 : for(gcp_index=0;gcp_index<nGCPCount;gcp_index++)
641 : {
642 7 : pasGCPList[gcp_index].dfGCPX = dfPrjX[gcp_index];
643 7 : pasGCPList[gcp_index].dfGCPY = dfPrjY[gcp_index];
644 :
645 : }
646 3 : transform_ok = GDALGCPsToGeoTransform(nGCPCount,pasGCPList,adfGeoTransform,0);
647 :
648 : }
649 :
650 3 : if (poTransform)
651 3 : delete poTransform;
652 :
653 3 : CPLFree(dfPrjX);
654 3 : CPLFree(dfPrjY);
655 :
656 : }
657 :
658 3 : CPLFree( pszProjection );
659 3 : CPLFree( pszGCPProjection );
660 3 : pszProjection = NULL;
661 3 : pszGCPProjection = NULL;
662 3 : oProj.exportToWkt( &pszProjection );
663 3 : oProj.exportToWkt( &pszGCPProjection );
664 :
665 3 : if (transform_ok == FALSE)
666 : {
667 : /* transform is sufficient in some cases (slant range, standalone gcps) */
668 2 : adfGeoTransform[0] = 0.0;
669 2 : adfGeoTransform[1] = 1.0;
670 2 : adfGeoTransform[2] = 0.0;
671 2 : adfGeoTransform[3] = 0.0;
672 2 : adfGeoTransform[4] = 0.0;
673 2 : adfGeoTransform[5] = 1.0;
674 2 : CPLFree( pszProjection );
675 2 : pszProjection = CPLStrdup("");
676 : }
677 :
678 3 : delete mffEllipsoids;
679 :
680 : }
681 :
682 :
683 : /************************************************************************/
684 : /* Open() */
685 : /************************************************************************/
686 :
687 9756 : GDALDataset *MFFDataset::Open( GDALOpenInfo * poOpenInfo )
688 :
689 : {
690 9756 : int i, bNative = TRUE;
691 : char **papszHdrLines;
692 :
693 : /* -------------------------------------------------------------------- */
694 : /* We assume the user is pointing to the header file. */
695 : /* -------------------------------------------------------------------- */
696 9756 : if( poOpenInfo->nHeaderBytes < 17 || poOpenInfo->fp == NULL )
697 9368 : return NULL;
698 :
699 388 : if( !EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"hdr") )
700 360 : return NULL;
701 :
702 : /* -------------------------------------------------------------------- */
703 : /* Load the .hdr file, and compress white space out around the */
704 : /* equal sign. */
705 : /* -------------------------------------------------------------------- */
706 28 : papszHdrLines = CSLLoad( poOpenInfo->pszFilename );
707 28 : if( papszHdrLines == NULL )
708 0 : return NULL;
709 :
710 207 : for( i = 0; papszHdrLines[i] != NULL; i++ )
711 : {
712 179 : int bAfterEqual = FALSE;
713 : int iSrc, iDst;
714 179 : char *pszLine = papszHdrLines[i];
715 :
716 3393 : for( iSrc=0, iDst=0; pszLine[iSrc] != '\0'; iSrc++ )
717 : {
718 3214 : if( bAfterEqual || pszLine[iSrc] != ' ' )
719 : {
720 2864 : pszLine[iDst++] = pszLine[iSrc];
721 : }
722 :
723 3214 : if( iDst > 0 && pszLine[iDst-1] == '=' )
724 281 : bAfterEqual = FALSE;
725 : }
726 179 : pszLine[iDst] = '\0';
727 : }
728 :
729 : /* -------------------------------------------------------------------- */
730 : /* Verify it is an MFF file. */
731 : /* -------------------------------------------------------------------- */
732 28 : if( CSLFetchNameValue( papszHdrLines, "IMAGE_FILE_FORMAT" ) != NULL
733 : && !EQUAL(CSLFetchNameValue(papszHdrLines,"IMAGE_FILE_FORMAT"),"MFF") )
734 : {
735 0 : CSLDestroy( papszHdrLines );
736 0 : return NULL;
737 : }
738 :
739 28 : if( (CSLFetchNameValue( papszHdrLines, "IMAGE_LINES" ) == NULL
740 : || CSLFetchNameValue(papszHdrLines,"LINE_SAMPLES") == NULL)
741 : && (CSLFetchNameValue( papszHdrLines, "no_rows" ) == NULL
742 : || CSLFetchNameValue(papszHdrLines,"no_columns") == NULL) )
743 : {
744 2 : CSLDestroy( papszHdrLines );
745 2 : return NULL;
746 : }
747 :
748 : /* -------------------------------------------------------------------- */
749 : /* Create a corresponding GDALDataset. */
750 : /* -------------------------------------------------------------------- */
751 : MFFDataset *poDS;
752 :
753 26 : poDS = new MFFDataset();
754 :
755 26 : poDS->papszHdrLines = papszHdrLines;
756 :
757 26 : poDS->eAccess = poOpenInfo->eAccess;
758 :
759 : /* -------------------------------------------------------------------- */
760 : /* Set some dataset wide information. */
761 : /* -------------------------------------------------------------------- */
762 27 : if( CSLFetchNameValue(papszHdrLines,"no_rows") != NULL
763 : && CSLFetchNameValue(papszHdrLines,"no_columns") != NULL )
764 : {
765 0 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszHdrLines,"no_columns"));
766 0 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszHdrLines,"no_rows"));
767 : }
768 : else
769 : {
770 26 : poDS->nRasterXSize = atoi(CSLFetchNameValue(papszHdrLines,"LINE_SAMPLES"));
771 26 : poDS->nRasterYSize = atoi(CSLFetchNameValue(papszHdrLines,"IMAGE_LINES"));
772 : }
773 :
774 26 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
775 : {
776 0 : delete poDS;
777 0 : return NULL;
778 : }
779 :
780 26 : if( CSLFetchNameValue( papszHdrLines, "BYTE_ORDER" ) != NULL )
781 : {
782 : #ifdef CPL_MSB
783 : bNative = EQUAL(CSLFetchNameValue(papszHdrLines,"BYTE_ORDER"),"MSB");
784 : #else
785 24 : bNative = EQUAL(CSLFetchNameValue(papszHdrLines,"BYTE_ORDER"),"LSB");
786 : #endif
787 : }
788 :
789 : /* -------------------------------------------------------------------- */
790 : /* Get some information specific to APP tiled files. */
791 : /* -------------------------------------------------------------------- */
792 26 : int bTiled, nTileXSize=0, nTileYSize=0;
793 26 : const char *pszRefinedType = NULL;
794 :
795 26 : pszRefinedType = CSLFetchNameValue(papszHdrLines, "type" );
796 :
797 26 : bTiled = CSLFetchNameValue(papszHdrLines,"no_rows") != NULL;
798 26 : if( bTiled )
799 : {
800 1 : if( CSLFetchNameValue(papszHdrLines,"tile_size_rows") )
801 : nTileYSize =
802 1 : atoi(CSLFetchNameValue(papszHdrLines,"tile_size_rows"));
803 1 : if( CSLFetchNameValue(papszHdrLines,"tile_size_columns") )
804 : nTileXSize =
805 1 : atoi(CSLFetchNameValue(papszHdrLines,"tile_size_columns"));
806 :
807 1 : if (nTileXSize <= 0 || nTileYSize <= 0 ||
808 : poDS->nRasterXSize > INT_MAX - (nTileXSize - 1) ||
809 : poDS->nRasterYSize > INT_MAX - (nTileYSize - 1))
810 : {
811 0 : delete poDS;
812 0 : return NULL;
813 : }
814 : }
815 :
816 : /* -------------------------------------------------------------------- */
817 : /* Read the directory to find matching band files. */
818 : /* -------------------------------------------------------------------- */
819 : char **papszDirFiles;
820 : char *pszTargetBase, *pszTargetPath;
821 26 : int nRawBand, nSkipped=0;
822 :
823 26 : pszTargetPath = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
824 26 : pszTargetBase = CPLStrdup(CPLGetBasename( poOpenInfo->pszFilename ));
825 26 : papszDirFiles = CPLReadDir( CPLGetPath( poOpenInfo->pszFilename ) );
826 26 : if( papszDirFiles == NULL )
827 : {
828 0 : CPLFree(pszTargetPath);
829 0 : CPLFree(pszTargetBase);
830 0 : delete poDS;
831 0 : return NULL;
832 : }
833 :
834 80 : for( nRawBand = 0; TRUE; nRawBand++ )
835 : {
836 : const char *pszExtension;
837 : int nBand;
838 : GDALDataType eDataType;
839 :
840 : /* Find the next raw band file. */
841 1800 : for( i = 0; papszDirFiles[i] != NULL; i++ )
842 : {
843 1774 : if( !EQUAL(CPLGetBasename(papszDirFiles[i]),pszTargetBase) )
844 1535 : continue;
845 :
846 239 : pszExtension = CPLGetExtension(papszDirFiles[i]);
847 239 : if( strlen(pszExtension) >= 2
848 : && isdigit(pszExtension[1])
849 : && atoi(pszExtension+1) == nRawBand
850 : && strchr("bBcCiIjJrRxXzZ",pszExtension[0]) != NULL )
851 54 : break;
852 : }
853 :
854 80 : if( papszDirFiles[i] == NULL )
855 26 : break;
856 :
857 : /* open the file for required level of access */
858 : FILE *fpRaw;
859 : const char *pszRawFilename = CPLFormFilename(pszTargetPath,
860 54 : papszDirFiles[i], NULL );
861 :
862 54 : if( poOpenInfo->eAccess == GA_Update )
863 50 : fpRaw = VSIFOpenL( pszRawFilename, "rb+" );
864 : else
865 4 : fpRaw = VSIFOpenL( pszRawFilename, "rb" );
866 :
867 54 : if( fpRaw == NULL )
868 : {
869 : CPLError( CE_Warning, CPLE_OpenFailed,
870 : "Unable to open %s ... skipping.\n",
871 0 : pszRawFilename );
872 0 : nSkipped++;
873 0 : continue;
874 : }
875 :
876 54 : pszExtension = CPLGetExtension(papszDirFiles[i]);
877 54 : if( pszRefinedType != NULL )
878 : {
879 0 : if( EQUAL(pszRefinedType,"C*4") )
880 0 : eDataType = GDT_CFloat32;
881 0 : else if( EQUAL(pszRefinedType,"C*8") )
882 0 : eDataType = GDT_CFloat64;
883 0 : else if( EQUAL(pszRefinedType,"R*4") )
884 0 : eDataType = GDT_Float32;
885 0 : else if( EQUAL(pszRefinedType,"R*8") )
886 0 : eDataType = GDT_Float64;
887 0 : else if( EQUAL(pszRefinedType,"I*1") )
888 0 : eDataType = GDT_Byte;
889 0 : else if( EQUAL(pszRefinedType,"I*2") )
890 0 : eDataType = GDT_Int16;
891 0 : else if( EQUAL(pszRefinedType,"I*4") )
892 0 : eDataType = GDT_Int32;
893 0 : else if( EQUAL(pszRefinedType,"U*2") )
894 0 : eDataType = GDT_UInt16;
895 0 : else if( EQUAL(pszRefinedType,"U*4") )
896 0 : eDataType = GDT_UInt32;
897 0 : else if( EQUAL(pszRefinedType,"J*1") )
898 : {
899 : CPLError( CE_Warning, CPLE_OpenFailed,
900 : "Unable to open band %d because type J*1 is not handled ... skipping.\n",
901 0 : nRawBand + 1 );
902 0 : nSkipped++;
903 0 : VSIFCloseL(fpRaw);
904 0 : continue; /* we don't support 1 byte complex */
905 : }
906 0 : else if( EQUAL(pszRefinedType,"J*2") )
907 0 : eDataType = GDT_CInt16;
908 0 : else if( EQUAL(pszRefinedType,"K*4") )
909 0 : eDataType = GDT_CInt32;
910 : else
911 : {
912 : CPLError( CE_Warning, CPLE_OpenFailed,
913 : "Unable to open band %d because type %s is not handled ... skipping.\n",
914 0 : nRawBand + 1, pszRefinedType );
915 0 : nSkipped++;
916 0 : VSIFCloseL(fpRaw);
917 0 : continue;
918 : }
919 : }
920 54 : else if( EQUALN(pszExtension,"b",1) )
921 : {
922 33 : eDataType = GDT_Byte;
923 : }
924 21 : else if( EQUALN(pszExtension,"i",1) )
925 : {
926 5 : eDataType = GDT_UInt16;
927 : }
928 16 : else if( EQUALN(pszExtension,"j",1) )
929 : {
930 5 : eDataType = GDT_CInt16;
931 : }
932 11 : else if( EQUALN(pszExtension,"r",1) )
933 : {
934 6 : eDataType = GDT_Float32;
935 : }
936 5 : else if( EQUALN(pszExtension,"x",1) )
937 : {
938 5 : eDataType = GDT_CFloat32;
939 : }
940 : else
941 : {
942 : CPLError( CE_Warning, CPLE_OpenFailed,
943 : "Unable to open band %d because extension %s is not handled ... skipping.\n",
944 0 : nRawBand + 1, pszExtension );
945 0 : nSkipped++;
946 0 : VSIFCloseL(fpRaw);
947 0 : continue;
948 : }
949 :
950 54 : nBand = poDS->GetRasterCount() + 1;
951 :
952 54 : int nPixelOffset = GDALGetDataTypeSize(eDataType)/8;
953 54 : GDALRasterBand *poBand = NULL;
954 :
955 54 : if( bTiled )
956 : {
957 : poBand =
958 : new MFFTiledBand( poDS, nBand, fpRaw, nTileXSize, nTileYSize,
959 1 : eDataType, bNative );
960 : }
961 : else
962 : {
963 53 : if (poDS->GetRasterXSize() > INT_MAX / nPixelOffset)
964 : {
965 0 : CPLError( CE_Warning, CPLE_AppDefined, "Int overflow occured... skipping");
966 0 : nSkipped++;
967 0 : VSIFCloseL(fpRaw);
968 0 : continue;
969 : }
970 :
971 : poBand =
972 : new RawRasterBand( poDS, nBand, fpRaw, 0, nPixelOffset,
973 : nPixelOffset * poDS->GetRasterXSize(),
974 53 : eDataType, bNative, TRUE, TRUE );
975 : }
976 :
977 54 : poDS->SetBand( nBand, poBand );
978 : }
979 :
980 26 : CPLFree(pszTargetPath);
981 26 : CPLFree(pszTargetBase);
982 26 : CSLDestroy(papszDirFiles);
983 :
984 : /* -------------------------------------------------------------------- */
985 : /* Check if we have bands. */
986 : /* -------------------------------------------------------------------- */
987 26 : if( poDS->GetRasterCount() == 0 )
988 : {
989 0 : if( nSkipped > 0 && poOpenInfo->eAccess )
990 : {
991 : CPLError( CE_Failure, CPLE_OpenFailed,
992 : "Failed to open %d files that were apparently bands.\n"
993 : "Perhaps this dataset is readonly?\n",
994 0 : nSkipped );
995 0 : delete poDS;
996 0 : return NULL;
997 : }
998 : else
999 : {
1000 : CPLError( CE_Failure, CPLE_OpenFailed,
1001 : "MFF header file read successfully, but no bands\n"
1002 0 : "were successfully found and opened." );
1003 0 : delete poDS;
1004 0 : return NULL;
1005 : }
1006 : }
1007 :
1008 : /* -------------------------------------------------------------------- */
1009 : /* Set all information from the .hdr that isn't well know to be */
1010 : /* metadata. */
1011 : /* -------------------------------------------------------------------- */
1012 189 : for( i = 0; papszHdrLines[i] != NULL; i++ )
1013 : {
1014 : const char *pszValue;
1015 : char *pszName;
1016 :
1017 163 : pszValue = CPLParseNameValue(papszHdrLines[i], &pszName);
1018 163 : if( pszName == NULL || pszValue == NULL )
1019 15 : continue;
1020 :
1021 148 : if( !EQUAL(pszName,"END")
1022 : && !EQUAL(pszName,"FILE_TYPE")
1023 : && !EQUAL(pszName,"BYTE_ORDER")
1024 : && !EQUAL(pszName,"no_columns")
1025 : && !EQUAL(pszName,"no_rows")
1026 : && !EQUAL(pszName,"type")
1027 : && !EQUAL(pszName,"tile_size_rows")
1028 : && !EQUAL(pszName,"tile_size_columns")
1029 : && !EQUAL(pszName,"IMAGE_FILE_FORMAT")
1030 : && !EQUAL(pszName,"IMAGE_LINES")
1031 : && !EQUAL(pszName,"LINE_SAMPLES") )
1032 : {
1033 19 : poDS->SetMetadataItem( pszName, pszValue );
1034 : }
1035 :
1036 148 : CPLFree( pszName );
1037 : }
1038 :
1039 : /* -------------------------------------------------------------------- */
1040 : /* Any GCPs in header file? */
1041 : /* -------------------------------------------------------------------- */
1042 26 : poDS->ScanForGCPs();
1043 26 : poDS->ScanForProjectionInfo();
1044 :
1045 : /* -------------------------------------------------------------------- */
1046 : /* Initialize any PAM information. */
1047 : /* -------------------------------------------------------------------- */
1048 26 : poDS->SetDescription( poOpenInfo->pszFilename );
1049 26 : poDS->TryLoadXML();
1050 :
1051 : /* -------------------------------------------------------------------- */
1052 : /* Check for overviews. */
1053 : /* -------------------------------------------------------------------- */
1054 26 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1055 :
1056 26 : return( poDS );
1057 : }
1058 :
1059 9 : int GetMFFProjectionType(const char *pszNewProjection)
1060 : {
1061 9 : OGRSpatialReference oSRS(pszNewProjection);
1062 :
1063 9 : if( !EQUALN(pszNewProjection,"GEOGCS",6)
1064 : && !EQUALN(pszNewProjection,"PROJCS",6)
1065 : && !EQUAL(pszNewProjection,"") )
1066 : {
1067 0 : return MFFPRJ_UNRECOGNIZED;
1068 : }
1069 9 : else if (EQUAL(pszNewProjection,""))
1070 : {
1071 0 : return MFFPRJ_NONE;
1072 : }
1073 : else
1074 : {
1075 9 : if ((oSRS.GetAttrValue("PROJECTION") != NULL) &&
1076 : (EQUAL(oSRS.GetAttrValue("PROJECTION"),SRS_PT_TRANSVERSE_MERCATOR)))
1077 : {
1078 0 : return MFFPRJ_UTM;
1079 : }
1080 9 : else if ((oSRS.GetAttrValue("PROJECTION") == NULL) && (oSRS.IsGeographic()))
1081 : {
1082 9 : return MFFPRJ_LL;
1083 : }
1084 : else
1085 : {
1086 0 : return MFFPRJ_UNRECOGNIZED;
1087 : }
1088 0 : }
1089 : }
1090 :
1091 : /************************************************************************/
1092 : /* Create() */
1093 : /************************************************************************/
1094 :
1095 : GDALDataset *MFFDataset::Create( const char * pszFilenameIn,
1096 : int nXSize, int nYSize, int nBands,
1097 : GDALDataType eType,
1098 41 : char ** papszParmList )
1099 :
1100 : {
1101 : /* -------------------------------------------------------------------- */
1102 : /* Verify input options. */
1103 : /* -------------------------------------------------------------------- */
1104 41 : if (nBands <= 0)
1105 : {
1106 : CPLError( CE_Failure, CPLE_NotSupported,
1107 1 : "MFF driver does not support %d bands.\n", nBands);
1108 1 : return NULL;
1109 : }
1110 :
1111 40 : if( eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16
1112 : && eType != GDT_CInt16 && eType != GDT_CFloat32 )
1113 : {
1114 : CPLError( CE_Failure, CPLE_AppDefined,
1115 : "Attempt to create MFF file with currently unsupported\n"
1116 : "data type (%s).\n",
1117 18 : GDALGetDataTypeName(eType) );
1118 :
1119 18 : return NULL;
1120 : }
1121 :
1122 : /* -------------------------------------------------------------------- */
1123 : /* Establish the base filename (path+filename, less extension). */
1124 : /* -------------------------------------------------------------------- */
1125 : char *pszBaseFilename;
1126 : int i;
1127 :
1128 22 : pszBaseFilename = (char *) CPLMalloc(strlen(pszFilenameIn)+5);
1129 22 : strcpy( pszBaseFilename, pszFilenameIn );
1130 :
1131 88 : for( i = strlen(pszBaseFilename)-1; i > 0; i-- )
1132 : {
1133 88 : if( pszBaseFilename[i] == '.' )
1134 : {
1135 0 : pszBaseFilename[i] = '\0';
1136 0 : break;
1137 : }
1138 :
1139 88 : if( pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\' )
1140 22 : break;
1141 : }
1142 :
1143 : /* -------------------------------------------------------------------- */
1144 : /* Create the header file. */
1145 : /* -------------------------------------------------------------------- */
1146 : FILE *fp;
1147 : const char *pszFilename;
1148 :
1149 22 : pszFilename = CPLFormFilename( NULL, pszBaseFilename, "hdr" );
1150 :
1151 22 : fp = VSIFOpen( pszFilename, "wt" );
1152 22 : if( fp == NULL )
1153 : {
1154 : CPLError( CE_Failure, CPLE_OpenFailed,
1155 0 : "Couldn't create %s.\n", pszFilename );
1156 0 : CPLFree(pszBaseFilename);
1157 0 : return NULL;
1158 : }
1159 :
1160 22 : fprintf( fp, "IMAGE_FILE_FORMAT = MFF\n" );
1161 22 : fprintf( fp, "FILE_TYPE = IMAGE\n" );
1162 22 : fprintf( fp, "IMAGE_LINES = %d\n", nYSize );
1163 22 : fprintf( fp, "LINE_SAMPLES = %d\n", nXSize );
1164 : #ifdef CPL_MSB
1165 : fprintf( fp, "BYTE_ORDER = MSB\n" );
1166 : #else
1167 22 : fprintf( fp, "BYTE_ORDER = LSB\n" );
1168 : #endif
1169 :
1170 22 : if (CSLFetchNameValue(papszParmList,"NO_END") == NULL)
1171 13 : fprintf( fp, "END\n" );
1172 :
1173 22 : VSIFClose( fp );
1174 :
1175 : /* -------------------------------------------------------------------- */
1176 : /* Create the data files, but don't bother writing any data to them.*/
1177 : /* -------------------------------------------------------------------- */
1178 72 : for( int iBand = 0; iBand < nBands; iBand++ )
1179 : {
1180 : char szExtension[4];
1181 :
1182 50 : if( eType == GDT_Byte )
1183 30 : sprintf( szExtension, "b%02d", iBand );
1184 20 : else if( eType == GDT_UInt16 )
1185 5 : sprintf( szExtension, "i%02d", iBand );
1186 15 : else if( eType == GDT_Float32 )
1187 5 : sprintf( szExtension, "r%02d", iBand );
1188 10 : else if( eType == GDT_CInt16 )
1189 5 : sprintf( szExtension, "j%02d", iBand );
1190 5 : else if( eType == GDT_CFloat32 )
1191 5 : sprintf( szExtension, "x%02d", iBand );
1192 :
1193 50 : pszFilename = CPLFormFilename( NULL, pszBaseFilename, szExtension );
1194 50 : fp = VSIFOpen( pszFilename, "wb" );
1195 50 : if( fp == NULL )
1196 : {
1197 : CPLError( CE_Failure, CPLE_OpenFailed,
1198 0 : "Couldn't create %s.\n", pszFilename );
1199 0 : CPLFree(pszBaseFilename);
1200 0 : return NULL;
1201 : }
1202 :
1203 50 : VSIFWrite( (void *) "", 1, 1, fp );
1204 50 : VSIFClose( fp );
1205 : }
1206 :
1207 : /* -------------------------------------------------------------------- */
1208 : /* Open the dataset normally. */
1209 : /* -------------------------------------------------------------------- */
1210 : GDALDataset *poDS;
1211 :
1212 22 : strcat( pszBaseFilename, ".hdr" );
1213 22 : poDS = (GDALDataset *) GDALOpen( pszBaseFilename, GA_Update );
1214 22 : CPLFree( pszBaseFilename );
1215 :
1216 22 : return poDS;
1217 : }
1218 :
1219 : /************************************************************************/
1220 : /* CreateCopy() */
1221 : /************************************************************************/
1222 :
1223 : GDALDataset *
1224 : MFFDataset::CreateCopy( const char * pszFilename, GDALDataset *poSrcDS,
1225 : int bStrict, char ** papszOptions,
1226 16 : GDALProgressFunc pfnProgress, void * pProgressData )
1227 :
1228 : {
1229 : MFFDataset *poDS;
1230 : GDALDataType eType;
1231 : int iBand;
1232 16 : char **newpapszOptions=NULL;
1233 :
1234 16 : int nBands = poSrcDS->GetRasterCount();
1235 16 : if (nBands == 0)
1236 : {
1237 : CPLError( CE_Failure, CPLE_NotSupported,
1238 1 : "MFF driver does not support source dataset with zero band.\n");
1239 1 : return NULL;
1240 : }
1241 :
1242 15 : eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
1243 15 : if( !pfnProgress( 0.0, NULL, pProgressData ) )
1244 0 : return NULL;
1245 :
1246 : /* check that other bands match type- sets type */
1247 : /* to unknown if they differ. */
1248 25 : for( iBand = 1; iBand < poSrcDS->GetRasterCount(); iBand++ )
1249 : {
1250 10 : GDALRasterBand *poBand = poSrcDS->GetRasterBand( iBand+1 );
1251 10 : eType = GDALDataTypeUnion( eType, poBand->GetRasterDataType() );
1252 : }
1253 :
1254 15 : newpapszOptions=CSLDuplicate(papszOptions);
1255 15 : newpapszOptions=CSLSetNameValue(newpapszOptions,"NO_END","TRUE");
1256 :
1257 : poDS = (MFFDataset *) Create( pszFilename,
1258 : poSrcDS->GetRasterXSize(),
1259 : poSrcDS->GetRasterYSize(),
1260 : poSrcDS->GetRasterCount(),
1261 15 : eType, newpapszOptions );
1262 :
1263 15 : CSLDestroy(newpapszOptions);
1264 :
1265 :
1266 : /* Check that Create worked- return Null if it didn't */
1267 15 : if (poDS == NULL)
1268 6 : return NULL;
1269 :
1270 :
1271 : /* -------------------------------------------------------------------- */
1272 : /* Copy the image data. */
1273 : /* -------------------------------------------------------------------- */
1274 9 : int nXSize = poDS->GetRasterXSize();
1275 9 : int nYSize = poDS->GetRasterYSize();
1276 : int nBlockXSize, nBlockYSize, nBlockTotal, nBlocksDone;
1277 :
1278 9 : poDS->GetRasterBand(1)->GetBlockSize( &nBlockXSize, &nBlockYSize );
1279 :
1280 : nBlockTotal = ((nXSize + nBlockXSize - 1) / nBlockXSize)
1281 : * ((nYSize + nBlockYSize - 1) / nBlockYSize)
1282 9 : * poSrcDS->GetRasterCount();
1283 :
1284 9 : nBlocksDone = 0;
1285 28 : for( iBand = 0; iBand < poSrcDS->GetRasterCount(); iBand++ )
1286 : {
1287 19 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand+1 );
1288 19 : GDALRasterBand *poDstBand = poDS->GetRasterBand( iBand+1 );
1289 : int iYOffset, iXOffset;
1290 : void *pData;
1291 : CPLErr eErr;
1292 :
1293 :
1294 : pData = CPLMalloc(nBlockXSize * nBlockYSize
1295 19 : * GDALGetDataTypeSize(eType) / 8);
1296 :
1297 209 : for( iYOffset = 0; iYOffset < nYSize; iYOffset += nBlockYSize )
1298 : {
1299 380 : for( iXOffset = 0; iXOffset < nXSize; iXOffset += nBlockXSize )
1300 : {
1301 : int nTBXSize, nTBYSize;
1302 :
1303 190 : if( !pfnProgress( (nBlocksDone++) / (float) nBlockTotal,
1304 : NULL, pProgressData ) )
1305 : {
1306 : CPLError( CE_Failure, CPLE_UserInterrupt,
1307 0 : "User terminated" );
1308 0 : delete poDS;
1309 :
1310 : GDALDriver *poMFFDriver =
1311 0 : (GDALDriver *) GDALGetDriverByName( "MFF" );
1312 0 : poMFFDriver->Delete( pszFilename );
1313 0 : return NULL;
1314 : }
1315 :
1316 190 : nTBXSize = MIN(nBlockXSize,nXSize-iXOffset);
1317 190 : nTBYSize = MIN(nBlockYSize,nYSize-iYOffset);
1318 :
1319 : eErr = poSrcBand->RasterIO( GF_Read,
1320 : iXOffset, iYOffset,
1321 : nTBXSize, nTBYSize,
1322 : pData, nTBXSize, nTBYSize,
1323 190 : eType, 0, 0 );
1324 190 : if( eErr != CE_None )
1325 : {
1326 0 : return NULL;
1327 : }
1328 :
1329 : eErr = poDstBand->RasterIO( GF_Write,
1330 : iXOffset, iYOffset,
1331 : nTBXSize, nTBYSize,
1332 : pData, nTBXSize, nTBYSize,
1333 190 : eType, 0, 0 );
1334 :
1335 190 : if( eErr != CE_None )
1336 : {
1337 0 : return NULL;
1338 : }
1339 : }
1340 : }
1341 :
1342 19 : CPLFree( pData );
1343 : }
1344 :
1345 : /* -------------------------------------------------------------------- */
1346 : /* Copy georeferencing information, if enough is available. */
1347 : /* -------------------------------------------------------------------- */
1348 :
1349 :
1350 : /* -------------------------------------------------------------------- */
1351 : /* Establish the base filename (path+filename, less extension). */
1352 : /* -------------------------------------------------------------------- */
1353 : char *pszBaseFilename;
1354 : int i;
1355 : FILE *fp;
1356 : const char *pszFilenameGEO;
1357 :
1358 9 : pszBaseFilename = (char *) CPLMalloc(strlen(pszFilename)+5);
1359 9 : strcpy( pszBaseFilename, pszFilename );
1360 :
1361 36 : for( i = strlen(pszBaseFilename)-1; i > 0; i-- )
1362 : {
1363 36 : if( pszBaseFilename[i] == '.' )
1364 : {
1365 0 : pszBaseFilename[i] = '\0';
1366 0 : break;
1367 : }
1368 :
1369 36 : if( pszBaseFilename[i] == '/' || pszBaseFilename[i] == '\\' )
1370 9 : break;
1371 : }
1372 :
1373 9 : pszFilenameGEO = CPLFormFilename( NULL, pszBaseFilename, "hdr" );
1374 :
1375 9 : fp = VSIFOpen( pszFilenameGEO, "at" );
1376 9 : if( fp == NULL )
1377 : {
1378 : CPLError( CE_Failure, CPLE_OpenFailed,
1379 0 : "Couldn't open %s for appending.\n", pszFilenameGEO );
1380 0 : CPLFree(pszBaseFilename);
1381 0 : return NULL;
1382 : }
1383 :
1384 :
1385 : /* MFF requires corner and center gcps */
1386 : double *padfTiepoints;
1387 : int src_prj;
1388 9 : int georef_created = FALSE;
1389 :
1390 9 : padfTiepoints = (double *) CPLMalloc(2*sizeof(double)*5);
1391 :
1392 9 : src_prj = GetMFFProjectionType(poSrcDS->GetProjectionRef());
1393 :
1394 9 : if ((src_prj != MFFPRJ_NONE) && (src_prj != MFFPRJ_UNRECOGNIZED))
1395 : {
1396 9 : double *tempGeoTransform = NULL;
1397 :
1398 9 : tempGeoTransform = (double *) CPLMalloc(6*sizeof(double));
1399 :
1400 9 : if (( poSrcDS->GetGeoTransform( tempGeoTransform ) == CE_None)
1401 : && (tempGeoTransform[0] != 0.0 || tempGeoTransform[1] != 1.0
1402 : || tempGeoTransform[2] != 0.0 || tempGeoTransform[3] != 0.0
1403 : || tempGeoTransform[4] != 0.0 || ABS(tempGeoTransform[5]) != 1.0 ))
1404 : {
1405 9 : OGRCoordinateTransformation *poTransform = NULL;
1406 9 : char *newGCPProjection=NULL;
1407 :
1408 9 : padfTiepoints[0]=tempGeoTransform[0] + tempGeoTransform[1]*0.5 +\
1409 : tempGeoTransform[2]*0.5;
1410 :
1411 9 : padfTiepoints[1]=tempGeoTransform[3] + tempGeoTransform[4]*0.5 +\
1412 : tempGeoTransform[5]*0.5;
1413 :
1414 9 : padfTiepoints[2]=tempGeoTransform[0] + tempGeoTransform[2]*0.5 +\
1415 : tempGeoTransform[1]*(poSrcDS->GetRasterXSize()-0.5);
1416 :
1417 9 : padfTiepoints[3]=tempGeoTransform[3] + tempGeoTransform[5]*0.5 +\
1418 : tempGeoTransform[4]*(poSrcDS->GetRasterXSize()-0.5);
1419 :
1420 9 : padfTiepoints[4]=tempGeoTransform[0] + tempGeoTransform[1]*0.5 +\
1421 : tempGeoTransform[2]*(poSrcDS->GetRasterYSize()-0.5);
1422 :
1423 9 : padfTiepoints[5]=tempGeoTransform[3] + tempGeoTransform[4]*0.5 +\
1424 : tempGeoTransform[5]*(poSrcDS->GetRasterYSize()-0.5);
1425 :
1426 9 : padfTiepoints[6]=tempGeoTransform[0] +\
1427 : tempGeoTransform[1]*(poSrcDS->GetRasterXSize()-0.5) +\
1428 : tempGeoTransform[2]*(poSrcDS->GetRasterYSize()-0.5);
1429 :
1430 9 : padfTiepoints[7]=tempGeoTransform[3]+\
1431 : tempGeoTransform[4]*(poSrcDS->GetRasterXSize()-0.5)+\
1432 : tempGeoTransform[5]*(poSrcDS->GetRasterYSize()-0.5);
1433 :
1434 9 : padfTiepoints[8]=tempGeoTransform[0]+\
1435 : tempGeoTransform[1]*(poSrcDS->GetRasterXSize())/2.0+\
1436 : tempGeoTransform[2]*(poSrcDS->GetRasterYSize())/2.0;
1437 :
1438 9 : padfTiepoints[9]=tempGeoTransform[3]+\
1439 : tempGeoTransform[4]*(poSrcDS->GetRasterXSize())/2.0+\
1440 : tempGeoTransform[5]*(poSrcDS->GetRasterYSize())/2.0;
1441 :
1442 9 : OGRSpatialReference oUTMorLL(poSrcDS->GetProjectionRef());
1443 9 : (oUTMorLL.GetAttrNode("GEOGCS"))->exportToWkt(&newGCPProjection);
1444 9 : OGRSpatialReference oLL(newGCPProjection);
1445 9 : CPLFree(newGCPProjection);
1446 9 : newGCPProjection = NULL;
1447 :
1448 9 : if EQUALN(poSrcDS->GetProjectionRef(),"PROJCS",6)
1449 : {
1450 : // projected coordinate system- need to translate gcps */
1451 0 : int bSuccess=TRUE;
1452 : int index;
1453 :
1454 0 : poTransform = OGRCreateCoordinateTransformation( &oUTMorLL, &oLL );
1455 0 : if( poTransform == NULL )
1456 0 : bSuccess = FALSE;
1457 :
1458 0 : for (index=0;index<5;index++)
1459 : {
1460 0 : if( !bSuccess || !poTransform->Transform( 1, &(padfTiepoints[index*2]), &(padfTiepoints[index*2+1]) ) )
1461 0 : bSuccess = FALSE;
1462 : }
1463 0 : if (bSuccess == TRUE)
1464 0 : georef_created = TRUE;
1465 : }
1466 : else
1467 : {
1468 9 : georef_created = TRUE;
1469 9 : }
1470 : }
1471 9 : CPLFree(tempGeoTransform);
1472 : }
1473 :
1474 9 : if (georef_created == TRUE)
1475 : {
1476 : /* -------------------------------------------------------------------- */
1477 : /* top left */
1478 : /* -------------------------------------------------------------------- */
1479 9 : fprintf( fp, "TOP_LEFT_CORNER_LATITUDE = %.10f\n", padfTiepoints[1] );
1480 9 : fprintf( fp, "TOP_LEFT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[0] );
1481 : /* -------------------------------------------------------------------- */
1482 : /* top_right */
1483 : /* -------------------------------------------------------------------- */
1484 9 : fprintf( fp, "TOP_RIGHT_CORNER_LATITUDE = %.10f\n", padfTiepoints[3] );
1485 9 : fprintf( fp, "TOP_RIGHT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[2] );
1486 : /* -------------------------------------------------------------------- */
1487 : /* bottom_left */
1488 : /* -------------------------------------------------------------------- */
1489 9 : fprintf( fp, "BOTTOM_LEFT_CORNER_LATITUDE = %.10f\n", padfTiepoints[5] );
1490 9 : fprintf( fp, "BOTTOM_LEFT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[4] );
1491 : /* -------------------------------------------------------------------- */
1492 : /* bottom_right */
1493 : /* -------------------------------------------------------------------- */
1494 9 : fprintf( fp, "BOTTOM_RIGHT_CORNER_LATITUDE = %.10f\n", padfTiepoints[7] );
1495 9 : fprintf( fp, "BOTTOM_RIGHT_CORNER_LONGITUDE = %.10f\n", padfTiepoints[6] );
1496 : /* -------------------------------------------------------------------- */
1497 : /* Center */
1498 : /* -------------------------------------------------------------------- */
1499 9 : fprintf( fp, "CENTRE_LATITUDE = %.10f\n", padfTiepoints[9] );
1500 9 : fprintf( fp, "CENTRE_LONGITUDE = %.10f\n", padfTiepoints[8] );
1501 : /* ------------------------------------------------------------------- */
1502 : /* Ellipsoid/projection */
1503 : /* --------------------------------------------------------------------*/
1504 :
1505 :
1506 : MFFSpheroidList *mffEllipsoids;
1507 : double eq_radius, inv_flattening;
1508 9 : OGRErr ogrerrorEq=OGRERR_NONE;
1509 9 : OGRErr ogrerrorInvf=OGRERR_NONE;
1510 9 : OGRErr ogrerrorOl=OGRERR_NONE;
1511 9 : const char *pszSrcProjection = poSrcDS->GetProjectionRef();
1512 9 : char *spheroid_name = NULL;
1513 :
1514 9 : if( !EQUALN(pszSrcProjection,"GEOGCS",6)
1515 : && !EQUALN(pszSrcProjection,"PROJCS",6)
1516 : && !EQUAL(pszSrcProjection,"") )
1517 : {
1518 : CPLError( CE_Warning, CPLE_AppDefined,
1519 : "Only OGC WKT Projections supported for writing to MFF.\n"
1520 : "%s not supported.",
1521 0 : pszSrcProjection );
1522 : }
1523 9 : else if (!EQUAL(pszSrcProjection,""))
1524 : {
1525 9 : OGRSpatialReference oSRS(pszSrcProjection);
1526 :
1527 9 : if ((oSRS.GetAttrValue("PROJECTION") != NULL) &&
1528 : (EQUAL(oSRS.GetAttrValue("PROJECTION"),SRS_PT_TRANSVERSE_MERCATOR)))
1529 : {
1530 0 : fprintf(fp,"PROJECTION_NAME = UTM\n");
1531 : fprintf(fp,"PROJECTION_ORIGIN_LONGITUDE = %f\n",
1532 0 : oSRS.GetProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0,&ogrerrorOl));
1533 : }
1534 9 : else if ((oSRS.GetAttrValue("PROJECTION") == NULL) && (oSRS.IsGeographic()))
1535 : {
1536 9 : fprintf(fp,"PROJECTION_NAME = LL\n");
1537 : }
1538 : else
1539 : {
1540 : CPLError( CE_Warning, CPLE_AppDefined,
1541 0 : "Unrecognized projection- no georeferencing information transferred.");
1542 0 : fprintf(fp,"PROJECTION_NAME = LL\n");
1543 : }
1544 9 : eq_radius = oSRS.GetSemiMajor(&ogrerrorEq);
1545 9 : inv_flattening = oSRS.GetInvFlattening(&ogrerrorInvf);
1546 9 : if ((ogrerrorEq == OGRERR_NONE) && (ogrerrorInvf == OGRERR_NONE))
1547 : {
1548 9 : mffEllipsoids = new MFFSpheroidList;
1549 18 : spheroid_name = mffEllipsoids->GetSpheroidNameByEqRadiusAndInvFlattening(eq_radius,inv_flattening);
1550 9 : if (spheroid_name != NULL)
1551 : {
1552 0 : fprintf(fp,"SPHEROID_NAME = %s\n",spheroid_name );
1553 : }
1554 : else
1555 : {
1556 : fprintf(fp,
1557 : "SPHEROID_NAME = USER_DEFINED\nSPHEROID_EQUATORIAL_RADIUS = %.10f\nSPHEROID_POLAR_RADIUS = %.10f\n",
1558 9 : eq_radius,eq_radius*(1-1.0/inv_flattening) );
1559 : }
1560 9 : delete mffEllipsoids;
1561 9 : CPLFree(spheroid_name);
1562 9 : }
1563 : }
1564 : }
1565 :
1566 9 : CPLFree( padfTiepoints );
1567 9 : fprintf( fp, "END\n" );
1568 9 : VSIFClose( fp );
1569 :
1570 : /* End of georeferencing stuff */
1571 :
1572 : /* Make sure image data gets flushed */
1573 28 : for( iBand = 0; iBand < poDS->GetRasterCount(); iBand++ )
1574 : {
1575 19 : RawRasterBand *poDstBand = (RawRasterBand *) poDS->GetRasterBand( iBand+1 );
1576 19 : poDstBand->FlushCache();
1577 : }
1578 :
1579 :
1580 9 : if( !pfnProgress( 1.0, NULL, pProgressData ) )
1581 : {
1582 : CPLError( CE_Failure, CPLE_UserInterrupt,
1583 0 : "User terminated" );
1584 0 : delete poDS;
1585 :
1586 : GDALDriver *poMFFDriver =
1587 0 : (GDALDriver *) GDALGetDriverByName( "MFF" );
1588 0 : poMFFDriver->Delete( pszFilename );
1589 0 : CPLFree(pszBaseFilename);
1590 0 : return NULL;
1591 : }
1592 :
1593 9 : poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
1594 9 : CPLFree(pszBaseFilename);
1595 :
1596 9 : return poDS;
1597 : }
1598 :
1599 :
1600 : /************************************************************************/
1601 : /* GDALRegister_MFF() */
1602 : /************************************************************************/
1603 :
1604 409 : void GDALRegister_MFF()
1605 :
1606 : {
1607 : GDALDriver *poDriver;
1608 :
1609 409 : if( GDALGetDriverByName( "MFF" ) == NULL )
1610 : {
1611 392 : poDriver = new GDALDriver();
1612 :
1613 392 : poDriver->SetDescription( "MFF" );
1614 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1615 392 : "Vexcel MFF Raster" );
1616 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1617 392 : "frmt_various.html#MFF" );
1618 392 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "hdr" );
1619 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1620 392 : "Byte UInt16 Float32 CInt16 CFloat32" );
1621 :
1622 392 : poDriver->pfnOpen = MFFDataset::Open;
1623 392 : poDriver->pfnCreate = MFFDataset::Create;
1624 392 : poDriver->pfnCreateCopy = MFFDataset::CreateCopy;
1625 :
1626 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
1627 : }
1628 409 : }
1629 :
|
