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