1 : /******************************************************************************
2 : * $Id: isis3dataset.cpp 10646 2007-01-18 02:38:10Z warmerdam $
3 : *
4 : * Project: ISIS Version 3 Driver
5 : * Purpose: Implementation of ISIS3Dataset
6 : * Author: Trent Hare (thare@usgs.gov)
7 : * Frank Warmerdam (warmerdam@pobox.com)
8 : *
9 : * NOTE: Original code authored by Trent and placed in the public domain as
10 : * per US government policy. I have (within my rights) appropriated it and
11 : * placed it under the following license. This is not intended to diminish
12 : * Trents contribution.
13 : ******************************************************************************
14 : * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
15 : *
16 : * Permission is hereby granted, free of charge, to any person obtaining a
17 : * copy of this software and associated documentation files (the "Software"),
18 : * to deal in the Software without restriction, including without limitation
19 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
20 : * and/or sell copies of the Software, and to permit persons to whom the
21 : * Software is furnished to do so, subject to the following conditions:
22 : *
23 : * The above copyright notice and this permission notice shall be included
24 : * in all copies or substantial portions of the Software.
25 : *
26 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
27 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
28 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
29 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
30 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
31 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
32 : * DEALINGS IN THE SOFTWARE.
33 : ****************************************************************************/
34 :
35 : #define NULL1 0
36 : #define NULL2 -32768
37 : //#define NULL3 0xFF7FFFFB //in hex
38 : //Same as ESRI_GRID_FLOAT_NO_DATA
39 : //#define NULL3 -340282346638528859811704183484516925440.0
40 : #define NULL3 -3.4028226550889044521e+38
41 :
42 : #ifndef PI
43 : # define PI 3.1415926535897932384626433832795
44 : #endif
45 :
46 : #include "rawdataset.h"
47 : #include "ogr_spatialref.h"
48 : #include "cpl_string.h"
49 : #include "nasakeywordhandler.h"
50 :
51 : CPL_CVSID("$Id: isis3dataset.cpp 10646 2007-09-18 02:38:10Z xxxx $");
52 :
53 : CPL_C_START
54 : void GDALRegister_ISIS3(void);
55 : CPL_C_END
56 :
57 : class ISIS3Dataset;
58 :
59 : /************************************************************************/
60 : /* ==================================================================== */
61 : /* ISISTiledBand */
62 : /* ==================================================================== */
63 : /************************************************************************/
64 :
65 : class ISISTiledBand : public GDALPamRasterBand
66 : {
67 : FILE *fpVSIL;
68 : GIntBig nFirstTileOffset;
69 : GIntBig nXTileOffset;
70 : GIntBig nYTileOffset;
71 : int bNativeOrder;
72 :
73 : public:
74 :
75 : ISISTiledBand( GDALDataset *poDS, FILE *fpVSIL,
76 : int nBand, GDALDataType eDT,
77 : int nTileXSize, int nTileYSize,
78 : GIntBig nFirstTileOffset,
79 : GIntBig nXTileOffset,
80 : GIntBig nYTileOffset,
81 : int bNativeOrder );
82 2 : virtual ~ISISTiledBand() {}
83 :
84 : virtual CPLErr IReadBlock( int, int, void * );
85 : };
86 :
87 : /************************************************************************/
88 : /* ISISTiledBand() */
89 : /************************************************************************/
90 :
91 1 : ISISTiledBand::ISISTiledBand( GDALDataset *poDS, FILE *fpVSIL,
92 : int nBand, GDALDataType eDT,
93 : int nTileXSize, int nTileYSize,
94 : GIntBig nFirstTileOffset,
95 : GIntBig nXTileOffset,
96 : GIntBig nYTileOffset,
97 1 : int bNativeOrder )
98 :
99 : {
100 1 : this->poDS = poDS;
101 1 : this->nBand = nBand;
102 1 : this->fpVSIL = fpVSIL;
103 1 : this->bNativeOrder = bNativeOrder;
104 1 : eDataType = eDT;
105 1 : nBlockXSize = nTileXSize;
106 1 : nBlockYSize = nTileYSize;
107 :
108 : int nBlocksPerRow =
109 1 : (poDS->GetRasterXSize() + nTileXSize - 1) / nTileXSize;
110 : int nBlocksPerColumn =
111 1 : (poDS->GetRasterYSize() + nTileYSize - 1) / nTileYSize;
112 :
113 1 : if( nXTileOffset == 0 && nYTileOffset == 0 )
114 : {
115 1 : nXTileOffset = (GDALGetDataTypeSize(eDT)/8) * nTileXSize * nTileYSize;
116 1 : nYTileOffset = nXTileOffset * nBlocksPerRow;
117 : }
118 :
119 : this->nFirstTileOffset = nFirstTileOffset
120 1 : + (nBand-1) * nYTileOffset * nBlocksPerColumn;
121 :
122 1 : this->nXTileOffset = nXTileOffset;
123 1 : this->nYTileOffset = nYTileOffset;
124 1 : }
125 :
126 : /************************************************************************/
127 : /* IReadBlock() */
128 : /************************************************************************/
129 :
130 2 : CPLErr ISISTiledBand::IReadBlock( int nXBlock, int nYBlock, void *pImage )
131 :
132 : {
133 : GIntBig nOffset = nFirstTileOffset +
134 2 : nXBlock * nXTileOffset + nYBlock * nYTileOffset;
135 : size_t nBlockSize =
136 2 : (GDALGetDataTypeSize(eDataType)/8) * nBlockXSize * nBlockYSize;
137 :
138 2 : if( VSIFSeekL( fpVSIL, nOffset, SEEK_SET ) != 0 )
139 : {
140 : CPLError( CE_Failure, CPLE_FileIO,
141 : "Failed to seek to offset %d to read tile %d,%d.",
142 0 : (int) nOffset, nXBlock, nYBlock );
143 0 : return CE_Failure;
144 : }
145 :
146 2 : if( VSIFReadL( pImage, 1, nBlockSize, fpVSIL ) != nBlockSize )
147 : {
148 : CPLError( CE_Failure, CPLE_FileIO,
149 : "Failed to read %d bytes for tile %d,%d.",
150 0 : (int) nBlockSize, nXBlock, nYBlock );
151 0 : return CE_Failure;
152 : }
153 :
154 2 : if( !bNativeOrder )
155 : GDALSwapWords( pImage, GDALGetDataTypeSize(eDataType)/8,
156 : nBlockXSize*nBlockYSize,
157 0 : GDALGetDataTypeSize(eDataType)/8 );
158 :
159 2 : return CE_None;
160 : }
161 :
162 : /************************************************************************/
163 : /* ==================================================================== */
164 : /* ISISDataset */
165 : /* ==================================================================== */
166 : /************************************************************************/
167 :
168 : class ISIS3Dataset : public RawDataset
169 : {
170 : FILE *fpImage; // image data file.
171 :
172 : CPLString osExternalCube;
173 :
174 : NASAKeywordHandler oKeywords;
175 :
176 : int bGotTransform;
177 : double adfGeoTransform[6];
178 :
179 : CPLString osProjection;
180 :
181 : int parse_label(const char *file, char *keyword, char *value);
182 : int strstrip(char instr[], char outstr[], int position);
183 :
184 : CPLString oTempResult;
185 :
186 : const char *GetKeyword( const char *pszPath,
187 : const char *pszDefault = "");
188 : const char *GetKeywordSub( const char *pszPath,
189 : int iSubscript,
190 : const char *pszDefault = "");
191 :
192 : public:
193 : ISIS3Dataset();
194 : ~ISIS3Dataset();
195 :
196 : virtual CPLErr GetGeoTransform( double * padfTransform );
197 : virtual const char *GetProjectionRef(void);
198 :
199 : virtual char **GetFileList();
200 :
201 : static int Identify( GDALOpenInfo * );
202 : static GDALDataset *Open( GDALOpenInfo * );
203 : static GDALDataset *Create( const char * pszFilename,
204 : int nXSize, int nYSize, int nBands,
205 : GDALDataType eType, char ** papszParmList );
206 : };
207 :
208 :
209 :
210 : /************************************************************************/
211 : /* ISIS3Dataset() */
212 : /************************************************************************/
213 :
214 2 : ISIS3Dataset::ISIS3Dataset()
215 : {
216 2 : fpImage = NULL;
217 2 : bGotTransform = FALSE;
218 2 : adfGeoTransform[0] = 0.0;
219 2 : adfGeoTransform[1] = 1.0;
220 2 : adfGeoTransform[2] = 0.0;
221 2 : adfGeoTransform[3] = 0.0;
222 2 : adfGeoTransform[4] = 0.0;
223 2 : adfGeoTransform[5] = 1.0;
224 2 : }
225 :
226 : /************************************************************************/
227 : /* ~ISIS3Dataset() */
228 : /************************************************************************/
229 :
230 4 : ISIS3Dataset::~ISIS3Dataset()
231 :
232 : {
233 2 : FlushCache();
234 2 : if( fpImage != NULL )
235 2 : VSIFCloseL( fpImage );
236 4 : }
237 :
238 : /************************************************************************/
239 : /* GetFileList() */
240 : /************************************************************************/
241 :
242 0 : char **ISIS3Dataset::GetFileList()
243 :
244 : {
245 0 : char **papszFileList = NULL;
246 :
247 0 : papszFileList = GDALPamDataset::GetFileList();
248 :
249 0 : if( strlen(osExternalCube) > 0 )
250 0 : papszFileList = CSLAddString( papszFileList, osExternalCube );
251 :
252 0 : return papszFileList;
253 : }
254 :
255 : /************************************************************************/
256 : /* GetProjectionRef() */
257 : /************************************************************************/
258 :
259 2 : const char *ISIS3Dataset::GetProjectionRef()
260 :
261 : {
262 2 : if( strlen(osProjection) > 0 )
263 2 : return osProjection;
264 : else
265 0 : return GDALPamDataset::GetProjectionRef();
266 : }
267 :
268 : /************************************************************************/
269 : /* GetGeoTransform() */
270 : /************************************************************************/
271 :
272 2 : CPLErr ISIS3Dataset::GetGeoTransform( double * padfTransform )
273 :
274 : {
275 2 : if( bGotTransform )
276 : {
277 2 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
278 2 : return CE_None;
279 : }
280 : else
281 : {
282 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
283 : }
284 : }
285 :
286 : /************************************************************************/
287 : /* Identify() */
288 : /************************************************************************/
289 9031 : int ISIS3Dataset::Identify( GDALOpenInfo * poOpenInfo )
290 :
291 : {
292 9031 : if( poOpenInfo->pabyHeader != NULL
293 : && strstr((const char *)poOpenInfo->pabyHeader,"IsisCube") != NULL )
294 2 : return TRUE;
295 : else
296 9029 : return FALSE;
297 : }
298 :
299 : /************************************************************************/
300 : /* Open() */
301 : /************************************************************************/
302 :
303 1302 : GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo )
304 : {
305 : /* -------------------------------------------------------------------- */
306 : /* Does this look like a CUBE dataset? */
307 : /* -------------------------------------------------------------------- */
308 1302 : if( !Identify( poOpenInfo ) )
309 1300 : return NULL;
310 :
311 : /* -------------------------------------------------------------------- */
312 : /* Open the file using the large file API. */
313 : /* -------------------------------------------------------------------- */
314 2 : FILE *fpQube = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
315 :
316 2 : if( fpQube == NULL )
317 0 : return NULL;
318 :
319 : ISIS3Dataset *poDS;
320 :
321 2 : poDS = new ISIS3Dataset();
322 :
323 2 : if( ! poDS->oKeywords.Ingest( fpQube, 0 ) )
324 : {
325 0 : VSIFCloseL( fpQube );
326 0 : delete poDS;
327 0 : return NULL;
328 : }
329 :
330 2 : VSIFCloseL( fpQube );
331 :
332 : /* -------------------------------------------------------------------- */
333 : /* Assume user is pointing to label (ie .lbl) file for detached option */
334 : /* -------------------------------------------------------------------- */
335 : // Image can be inline or detached and point to an image name
336 : // the Format can be Tiled or Raw
337 : // Object = Core
338 : // StartByte = 65537
339 : // Format = Tile
340 : // TileSamples = 128
341 : // TileLines = 128
342 : //OR-----
343 : // Object = Core
344 : // StartByte = 1
345 : // ^Core = r0200357_detatched.cub
346 : // Format = BandSequential
347 : //OR-----
348 : // Object = Core
349 : // StartByte = 1
350 : // ^Core = r0200357_detached_tiled.cub
351 : // Format = Tile
352 : // TileSamples = 128
353 : // TileLines = 128
354 :
355 : /* -------------------------------------------------------------------- */
356 : /* What file contains the actual data? */
357 : /* -------------------------------------------------------------------- */
358 2 : const char *pszCore = poDS->GetKeyword( "IsisCube.Core.^Core" );
359 2 : CPLString osQubeFile;
360 :
361 2 : if( EQUAL(pszCore,"") )
362 1 : osQubeFile = poOpenInfo->pszFilename;
363 : else
364 : {
365 1 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
366 1 : osQubeFile = CPLFormFilename( osPath, pszCore, NULL );
367 1 : poDS->osExternalCube = osQubeFile;
368 : }
369 :
370 : /* -------------------------------------------------------------------- */
371 : /* Check if file an ISIS3 header file? Read a few lines of text */
372 : /* searching for something starting with nrows or ncols. */
373 : /* -------------------------------------------------------------------- */
374 2 : GDALDataType eDataType = GDT_Byte;
375 2 : OGRSpatialReference oSRS;
376 :
377 2 : int nRows = -1;
378 2 : int nCols = -1;
379 2 : int nBands = 1;
380 2 : int nSkipBytes = 0;
381 2 : int tileSizeX = 0;
382 2 : int tileSizeY = 0;
383 2 : double dfULXMap=0.5;
384 2 : double dfULYMap = 0.5;
385 2 : double dfXDim = 1.0;
386 2 : double dfYDim = 1.0;
387 2 : double scaleFactor = 1.0;
388 2 : double dfNoData = 0.0;
389 2 : int bNoDataSet = FALSE;
390 2 : char chByteOrder = 'M'; //default to MSB
391 2 : char szLayout[32] = "BandSequential"; //default to band seq.
392 : const char *target_name; //planet name
393 : //projection parameters
394 : const char *map_proj_name;
395 2 : int bProjectionSet = TRUE;
396 : char proj_target_name[200];
397 : char geog_name[60];
398 : char datum_name[60];
399 : char sphere_name[60];
400 2 : char bIsGeographic = TRUE;
401 2 : double semi_major = 0.0;
402 2 : double semi_minor = 0.0;
403 2 : double iflattening = 0.0;
404 2 : float center_lat = 0.0;
405 2 : float center_lon = 0.0;
406 2 : float first_std_parallel = 0.0;
407 2 : float second_std_parallel = 0.0;
408 : double radLat, localRadius;
409 : FILE *fp;
410 :
411 : /************* Skipbytes *****************************/
412 2 : nSkipBytes = atoi(poDS->GetKeyword("IsisCube.Core.StartByte","")) - 1;
413 :
414 : /******* Grab format type (BandSequential, Tiled) *******/
415 : const char *value;
416 :
417 2 : value = poDS->GetKeyword( "IsisCube.Core.Format", "" );
418 2 : if (EQUAL(value,"Tile") ) { //Todo
419 1 : strcpy(szLayout,"Tiled");
420 : /******* Get Tile Sizes *********/
421 1 : tileSizeX = atoi(poDS->GetKeyword("IsisCube.Core.TileSamples",""));
422 1 : tileSizeY = atoi(poDS->GetKeyword("IsisCube.Core.TileLines",""));
423 1 : if (tileSizeX <= 0 || tileSizeY <= 0)
424 : {
425 : CPLError( CE_Failure, CPLE_OpenFailed, "Wrong tile dimensions : %d x %d",
426 0 : tileSizeX, tileSizeY);
427 0 : delete poDS;
428 0 : return NULL;
429 : }
430 : }
431 1 : else if (EQUAL(value,"BandSequential") )
432 1 : strcpy(szLayout,"BSQ");
433 : else {
434 : CPLError( CE_Failure, CPLE_OpenFailed,
435 0 : "%s layout not supported. Abort\n\n", value);
436 0 : delete poDS;
437 0 : return NULL;
438 : }
439 :
440 : /*********** Grab samples lines band ************/
441 2 : nCols = atoi(poDS->GetKeyword("IsisCube.Core.Dimensions.Samples",""));
442 2 : nRows = atoi(poDS->GetKeyword("IsisCube.Core.Dimensions.Lines",""));
443 2 : nBands = atoi(poDS->GetKeyword("IsisCube.Core.Dimensions.Bands",""));
444 :
445 : /****** Grab format type - ISIS3 only supports 8,U16,S16,32 *****/
446 : const char *itype;
447 :
448 2 : itype = poDS->GetKeyword( "IsisCube.Core.Pixels.Type" );
449 2 : if (EQUAL(itype,"UnsignedByte") ) {
450 1 : eDataType = GDT_Byte;
451 1 : dfNoData = NULL1;
452 1 : bNoDataSet = TRUE;
453 : }
454 1 : else if (EQUAL(itype,"UnsignedWord") ) {
455 0 : eDataType = GDT_UInt16;
456 0 : dfNoData = NULL1;
457 0 : bNoDataSet = TRUE;
458 : }
459 1 : else if (EQUAL(itype,"SignedWord") ) {
460 1 : eDataType = GDT_Int16;
461 1 : dfNoData = NULL2;
462 1 : bNoDataSet = TRUE;
463 : }
464 0 : else if (EQUAL(itype,"Real") || EQUAL(value,"") ) {
465 0 : eDataType = GDT_Float32;
466 0 : dfNoData = NULL3;
467 0 : bNoDataSet = TRUE;
468 : }
469 : else {
470 : CPLError( CE_Failure, CPLE_OpenFailed,
471 0 : "%s layout type not supported. Abort\n\n", itype);
472 0 : delete poDS;
473 0 : return NULL;
474 : }
475 :
476 : /*********** Grab samples lines band ************/
477 2 : value = poDS->GetKeyword( "IsisCube.Core.Pixels.ByteOrder");
478 2 : if (EQUAL(value,"Lsb"))
479 2 : chByteOrder = 'I';
480 :
481 : /*********** Grab Cellsize ************/
482 2 : value = poDS->GetKeyword("IsisCube.Mapping.PixelResolution");
483 2 : if (strlen(value) > 0 ) {
484 2 : dfXDim = atof(value); /* values are in meters */
485 2 : dfYDim = -atof(value);
486 : }
487 :
488 : /*********** Grab UpperLeftCornerY ************/
489 2 : value = poDS->GetKeyword("IsisCube.Mapping.UpperLeftCornerY");
490 2 : if (strlen(value) > 0) {
491 2 : dfULYMap = atof(value);
492 : }
493 :
494 : /*********** Grab UpperLeftCornerX ************/
495 2 : value = poDS->GetKeyword("IsisCube.Mapping.UpperLeftCornerX");
496 2 : if( strlen(value) > 0 ) {
497 2 : dfULXMap = atof(value);
498 : }
499 :
500 : /*********** Grab TARGET_NAME ************/
501 : /**** This is the planets name i.e. Mars ***/
502 2 : target_name = poDS->GetKeyword("IsisCube.Mapping.TargetName");
503 :
504 : /*********** Grab MAP_PROJECTION_TYPE ************/
505 : map_proj_name =
506 2 : poDS->GetKeyword( "IsisCube.Mapping.ProjectionName");
507 :
508 : /*********** Grab SEMI-MAJOR ************/
509 : semi_major =
510 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.EquatorialRadius"));
511 :
512 : /*********** Grab semi-minor ************/
513 : semi_minor =
514 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.PolarRadius"));
515 :
516 : /*********** Grab CENTER_LAT ************/
517 : center_lat =
518 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.CenterLatitude"));
519 :
520 : /*********** Grab CENTER_LON ************/
521 : center_lon =
522 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.CenterLongitude"));
523 :
524 : /*********** Grab 1st std parallel ************/
525 : first_std_parallel =
526 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.FirstStandardParallel"));
527 :
528 : /*********** Grab 2nd std parallel ************/
529 : second_std_parallel =
530 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.SecondStandardParallel"));
531 :
532 : /*********** Grab scaleFactor ************/
533 : scaleFactor =
534 2 : atof(poDS->GetKeyword( "IsisCube.Mapping.scaleFactor"));
535 :
536 : /*** grab LatitudeType = Planetographic ****/
537 : // Need to further study how ocentric/ographic will effect the gdal library
538 : // So far we will use this fact to define a sphere or ellipse for some
539 : // projections
540 :
541 : // Frank - may need to talk this over
542 2 : value = poDS->GetKeyword("IsisCube.Mapping.LatitudeType");
543 2 : if (EQUAL( value, "\"Planetocentric\"" ))
544 0 : bIsGeographic = FALSE;
545 :
546 : //Set oSRS projection and parameters
547 : //############################################################
548 : //ISIS3 Projection types
549 : // Equirectangular
550 : // LambertConformal
551 : // Mercator
552 : // ObliqueCylindrical //Todo
553 : // Orthographic
554 : // PolarStereographic
555 : // SimpleCylindrical
556 : // Sinusoidal
557 : // TransverseMercator
558 :
559 : #ifdef DEBUG
560 : CPLDebug( "ISIS3", "using projection %s", map_proj_name);
561 : #endif
562 :
563 4 : if ((EQUAL( map_proj_name, "Equirectangular" )) ||
564 : (EQUAL( map_proj_name, "SimpleCylindrical" )) ) {
565 2 : oSRS.OGRSpatialReference::SetEquirectangular2 ( 0.0, center_lon, center_lat, 0, 0 );
566 0 : } else if (EQUAL( map_proj_name, "Orthographic" )) {
567 0 : oSRS.OGRSpatialReference::SetOrthographic ( center_lat, center_lon, 0, 0 );
568 0 : } else if (EQUAL( map_proj_name, "Sinusoidal" )) {
569 0 : oSRS.OGRSpatialReference::SetSinusoidal ( center_lon, 0, 0 );
570 0 : } else if (EQUAL( map_proj_name, "Mercator" )) {
571 0 : oSRS.OGRSpatialReference::SetMercator ( center_lat, center_lon, scaleFactor, 0, 0 );
572 0 : } else if (EQUAL( map_proj_name, "PolarStereographic" )) {
573 0 : oSRS.OGRSpatialReference::SetPS ( center_lat, center_lon, scaleFactor, 0, 0 );
574 0 : } else if (EQUAL( map_proj_name, "TransverseMercator" )) {
575 0 : oSRS.OGRSpatialReference::SetTM ( center_lat, center_lon, scaleFactor, 0, 0 );
576 0 : } else if (EQUAL( map_proj_name, "LambertConformal" )) {
577 0 : oSRS.OGRSpatialReference::SetLCC ( first_std_parallel, second_std_parallel, center_lat, center_lon, 0, 0 );
578 : } else {
579 : CPLDebug( "ISIS3",
580 : "Dataset projection %s is not supported. Continuing...",
581 0 : map_proj_name );
582 0 : bProjectionSet = FALSE;
583 : }
584 :
585 2 : if (bProjectionSet) {
586 : //Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
587 2 : strcpy(proj_target_name, map_proj_name);
588 2 : strcat(proj_target_name, " ");
589 2 : strcat(proj_target_name, target_name);
590 2 : oSRS.SetProjCS(proj_target_name); //set ProjCS keyword
591 :
592 : //The geographic/geocentric name will be the same basic name as the body name
593 : //'GCS' = Geographic/Geocentric Coordinate System
594 2 : strcpy(geog_name, "GCS_");
595 2 : strcat(geog_name, target_name);
596 :
597 : //The datum name will be the same basic name as the planet
598 2 : strcpy(datum_name, "D_");
599 2 : strcat(datum_name, target_name);
600 :
601 2 : strcpy(sphere_name, target_name);
602 : //strcat(sphere_name, "_IAU_IAG"); //Might not be IAU defined so don't add
603 :
604 : //calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
605 2 : if ((semi_major - semi_minor) < 0.0000001)
606 0 : iflattening = 0;
607 : else
608 2 : iflattening = semi_major / (semi_major - semi_minor);
609 :
610 : //Set the body size but take into consideration which proj is being used to help w/ proj4 compatibility
611 : //The use of a Sphere, polar radius or ellipse here is based on how ISIS does it internally
612 2 : if ( ( (EQUAL( map_proj_name, "Stereographic" ) && (fabs(center_lat) == 90)) ) ||
613 : (EQUAL( map_proj_name, "PolarStereographic" )) )
614 : {
615 0 : if (bIsGeographic) {
616 : //Geograpraphic, so set an ellipse
617 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
618 : semi_major, iflattening,
619 0 : "Reference_Meridian", 0.0 );
620 : } else {
621 : //Geocentric, so force a sphere using the semi-minor axis. I hope...
622 0 : strcat(sphere_name, "_polarRadius");
623 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
624 : semi_minor, 0.0,
625 0 : "Reference_Meridian", 0.0 );
626 : }
627 : }
628 2 : else if ( (EQUAL( map_proj_name, "SimpleCylindrical" )) ||
629 : (EQUAL( map_proj_name, "Orthographic" )) ||
630 : (EQUAL( map_proj_name, "Stereographic" )) ||
631 : (EQUAL( map_proj_name, "Sinusoidal" )) ) {
632 : //isis uses the sphereical equation for these projections so force a sphere
633 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
634 : semi_major, 0.0,
635 0 : "Reference_Meridian", 0.0 );
636 : }
637 2 : else if (EQUAL( map_proj_name, "Equirectangular" )) {
638 : //Calculate localRadius using ISIS3 simple elliptical method
639 : // not the more standard Radius of Curvature method
640 : //PI = 4 * atan(1);
641 2 : radLat = center_lat * PI / 180; // in radians
642 : localRadius = semi_major * semi_minor / sqrt(pow(semi_minor*cos(radLat),2)
643 2 : + pow(semi_major*sin(radLat),2) );
644 2 : strcat(sphere_name, "_localRadius");
645 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
646 : localRadius, 0.0,
647 2 : "Reference_Meridian", 0.0 );
648 : }
649 : else {
650 : //All other projections: Mercator, Transverse Mercator, Lambert Conformal, etc.
651 : //Geographic, so set an ellipse
652 0 : if (bIsGeographic) {
653 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
654 : semi_major, iflattening,
655 0 : "Reference_Meridian", 0.0 );
656 : } else {
657 : //Geocentric, so force a sphere. I hope...
658 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
659 : semi_major, 0.0,
660 0 : "Reference_Meridian", 0.0 );
661 : }
662 : }
663 :
664 : // translate back into a projection string.
665 2 : char *pszResult = NULL;
666 2 : oSRS.exportToWkt( &pszResult );
667 2 : poDS->osProjection = pszResult;
668 2 : CPLFree( pszResult );
669 : }
670 :
671 : /* END ISIS3 Label Read */
672 : /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
673 :
674 : /* -------------------------------------------------------------------- */
675 : /* Is the CUB detached - if so, reset name to binary file? */
676 : /* -------------------------------------------------------------------- */
677 : #ifdef notdef
678 : // Frank - is this correct?
679 : //The extension already added on so don't add another. But is this needed?
680 : char *pszPath = CPLStrdup( CPLGetPath( poOpenInfo->pszFilename ) );
681 : char *pszName = CPLStrdup( CPLGetBasename( poOpenInfo->pszFilename ) );
682 : if (bIsDetached)
683 : pszCUBFilename = CPLFormCIFilename( pszPath, detachedCub, "" );
684 : #endif
685 :
686 : /* -------------------------------------------------------------------- */
687 : /* Did we get the required keywords? If not we return with */
688 : /* this never having been considered to be a match. This isn't */
689 : /* an error! */
690 : /* -------------------------------------------------------------------- */
691 2 : if( nRows < 1 || nCols < 1 || nBands < 1 )
692 : {
693 0 : delete poDS;
694 0 : return NULL;
695 : }
696 :
697 : /* -------------------------------------------------------------------- */
698 : /* Capture some information from the file that is of interest. */
699 : /* -------------------------------------------------------------------- */
700 2 : poDS->nRasterXSize = nCols;
701 2 : poDS->nRasterYSize = nRows;
702 :
703 : /* -------------------------------------------------------------------- */
704 : /* Open target binary file. */
705 : /* -------------------------------------------------------------------- */
706 2 : if( poOpenInfo->eAccess == GA_ReadOnly )
707 2 : poDS->fpImage = VSIFOpenL( osQubeFile, "r" );
708 : else
709 0 : poDS->fpImage = VSIFOpenL( osQubeFile, "r+" );
710 :
711 2 : if( poDS->fpImage == NULL )
712 : {
713 : CPLError( CE_Failure, CPLE_OpenFailed,
714 : "Failed to open %s with write permission.\n%s",
715 : osQubeFile.c_str(),
716 0 : VSIStrerror( errno ) );
717 0 : delete poDS;
718 0 : return NULL;
719 : }
720 :
721 2 : poDS->eAccess = poOpenInfo->eAccess;
722 :
723 : /* -------------------------------------------------------------------- */
724 : /* Compute the line offset. */
725 : /* -------------------------------------------------------------------- */
726 2 : int nItemSize = GDALGetDataTypeSize(eDataType)/8;
727 2 : int nLineOffset=0, nPixelOffset=0, nBandOffset=0;
728 :
729 2 : if( EQUAL(szLayout,"BSQ") )
730 : {
731 1 : nPixelOffset = nItemSize;
732 1 : nLineOffset = nPixelOffset * nCols;
733 1 : nBandOffset = nLineOffset * nRows;
734 : }
735 : else /* Tiled */
736 : {
737 : }
738 :
739 : /* -------------------------------------------------------------------- */
740 : /* Create band information objects. */
741 : /* -------------------------------------------------------------------- */
742 : int i;
743 :
744 : #ifdef CPL_LSB
745 2 : int bNativeOrder = !(chByteOrder == 'M');
746 : #else
747 : int bNativeOrder = (chByteOrder == 'M');
748 : #endif
749 :
750 :
751 4 : for( i = 0; i < nBands; i++ )
752 : {
753 : GDALRasterBand *poBand;
754 :
755 2 : if( EQUAL(szLayout,"Tiled") )
756 : {
757 : poBand = new ISISTiledBand( poDS, poDS->fpImage, i+1, eDataType,
758 : tileSizeX, tileSizeY,
759 : nSkipBytes, 0, 0,
760 1 : bNativeOrder );
761 : }
762 : else
763 : {
764 : poBand =
765 : new RawRasterBand( poDS, i+1, poDS->fpImage,
766 : nSkipBytes + nBandOffset * i,
767 : nPixelOffset, nLineOffset, eDataType,
768 : #ifdef CPL_LSB
769 : chByteOrder == 'I' || chByteOrder == 'L',
770 : #else
771 : chByteOrder == 'M',
772 : #endif
773 1 : TRUE );
774 : }
775 :
776 2 : poDS->SetBand( i+1, poBand );
777 :
778 2 : if( bNoDataSet )
779 2 : ((GDALPamRasterBand *) poBand)->SetNoDataValue( dfNoData );
780 :
781 : // Set offset/scale values at the PAM level.
782 : poBand->SetOffset(
783 2 : CPLAtofM(poDS->GetKeyword("IsisCube.Core.Pixels.Base","0.0")));
784 : poBand->SetScale(
785 2 : CPLAtofM(poDS->GetKeyword("IsisCube.Core.Pixels.Multiplier","1.0")));
786 : }
787 :
788 : /* -------------------------------------------------------------------- */
789 : /* Check for a .prj file. For ISIS3 I would like to keep this in */
790 : /* -------------------------------------------------------------------- */
791 2 : CPLString osPath, osName;
792 :
793 2 : osPath = CPLGetPath( poOpenInfo->pszFilename );
794 2 : osName = CPLGetBasename(poOpenInfo->pszFilename);
795 2 : const char *pszPrjFile = CPLFormCIFilename( osPath, osName, "prj" );
796 :
797 2 : fp = VSIFOpen( pszPrjFile, "r" );
798 2 : if( fp != NULL )
799 : {
800 : char **papszLines;
801 0 : OGRSpatialReference oSRS;
802 :
803 0 : VSIFClose( fp );
804 :
805 0 : papszLines = CSLLoad( pszPrjFile );
806 :
807 0 : if( oSRS.importFromESRI( papszLines ) == OGRERR_NONE )
808 : {
809 0 : char *pszResult = NULL;
810 0 : oSRS.exportToWkt( &pszResult );
811 0 : poDS->osProjection = pszResult;
812 0 : CPLFree( pszResult );
813 : }
814 :
815 0 : CSLDestroy( papszLines );
816 : }
817 :
818 :
819 2 : if( dfULYMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0 )
820 : {
821 2 : poDS->bGotTransform = TRUE;
822 2 : poDS->adfGeoTransform[0] = dfULXMap;
823 2 : poDS->adfGeoTransform[1] = dfXDim;
824 2 : poDS->adfGeoTransform[2] = 0.0;
825 2 : poDS->adfGeoTransform[3] = dfULYMap;
826 2 : poDS->adfGeoTransform[4] = 0.0;
827 2 : poDS->adfGeoTransform[5] = dfYDim;
828 : }
829 :
830 2 : if( !poDS->bGotTransform )
831 : poDS->bGotTransform =
832 : GDALReadWorldFile( poOpenInfo->pszFilename, "cbw",
833 0 : poDS->adfGeoTransform );
834 :
835 2 : if( !poDS->bGotTransform )
836 : poDS->bGotTransform =
837 : GDALReadWorldFile( poOpenInfo->pszFilename, "wld",
838 0 : poDS->adfGeoTransform );
839 :
840 : /* -------------------------------------------------------------------- */
841 : /* Initialize any PAM information. */
842 : /* -------------------------------------------------------------------- */
843 2 : poDS->SetDescription( poOpenInfo->pszFilename );
844 2 : poDS->TryLoadXML();
845 :
846 : /* -------------------------------------------------------------------- */
847 : /* Check for overviews. */
848 : /* -------------------------------------------------------------------- */
849 2 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
850 :
851 2 : return( poDS );
852 : }
853 :
854 : /************************************************************************/
855 : /* GetKeyword() */
856 : /************************************************************************/
857 :
858 48 : const char *ISIS3Dataset::GetKeyword( const char *pszPath,
859 : const char *pszDefault )
860 :
861 : {
862 48 : return oKeywords.GetKeyword( pszPath, pszDefault );
863 : }
864 :
865 : /************************************************************************/
866 : /* GetKeywordSub() */
867 : /************************************************************************/
868 :
869 0 : const char *ISIS3Dataset::GetKeywordSub( const char *pszPath,
870 : int iSubscript,
871 : const char *pszDefault )
872 :
873 : {
874 0 : const char *pszResult = oKeywords.GetKeyword( pszPath, NULL );
875 :
876 0 : if( pszResult == NULL )
877 0 : return pszDefault;
878 :
879 0 : if( pszResult[0] != '(' )
880 0 : return pszDefault;
881 :
882 : char **papszTokens = CSLTokenizeString2( pszResult, "(,)",
883 0 : CSLT_HONOURSTRINGS );
884 :
885 0 : if( iSubscript <= CSLCount(papszTokens) )
886 : {
887 0 : oTempResult = papszTokens[iSubscript-1];
888 0 : CSLDestroy( papszTokens );
889 0 : return oTempResult.c_str();
890 : }
891 : else
892 : {
893 0 : CSLDestroy( papszTokens );
894 0 : return pszDefault;
895 : }
896 : }
897 :
898 :
899 : /************************************************************************/
900 : /* GDALRegister_ISIS3() */
901 : /************************************************************************/
902 :
903 338 : void GDALRegister_ISIS3()
904 :
905 : {
906 : GDALDriver *poDriver;
907 :
908 338 : if( GDALGetDriverByName( "ISIS3" ) == NULL )
909 : {
910 336 : poDriver = new GDALDriver();
911 :
912 336 : poDriver->SetDescription( "ISIS3" );
913 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
914 336 : "USGS Astrogeology ISIS cube (Version 3)" );
915 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
916 336 : "frmt_various.html#ISIS3" );
917 :
918 336 : poDriver->pfnOpen = ISIS3Dataset::Open;
919 336 : poDriver->pfnIdentify = ISIS3Dataset::Identify;
920 :
921 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
922 : }
923 338 : }
924 :
|
