1 : /******************************************************************************
2 : * $Id: isis2dataset.cpp 17786 2009-10-09 19:35:19Z warmerdam $
3 : *
4 : * Project: ISIS Version 2 Driver
5 : * Purpose: Implementation of ISIS2Dataset
6 : * Author: Trent Hare (thare@usgs.gov),
7 : * Robert Soricone (rsoricone@usgs.gov)
8 : *
9 : * NOTE: Original code authored by Trent and Robert and placed in the public
10 : * domain as per US government policy. I have (within my rights) appropriated
11 : * it and placed it under the following license. This is not intended to
12 : * diminish Trent and Roberts contribution.
13 : ******************************************************************************
14 : * Copyright (c) 2006, 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: isis2dataset.cpp 17786 2009-10-09 19:35:19Z warmerdam $");
52 :
53 : CPL_C_START
54 : void GDALRegister_ISIS2(void);
55 : CPL_C_END
56 :
57 : /************************************************************************/
58 : /* ==================================================================== */
59 : /* ISISDataset version2 */
60 : /* ==================================================================== */
61 : /************************************************************************/
62 :
63 : class ISIS2Dataset : public RawDataset
64 : {
65 : FILE *fpImage; // image data file.
66 :
67 : NASAKeywordHandler oKeywords;
68 :
69 : int bGotTransform;
70 : double adfGeoTransform[6];
71 :
72 : CPLString osProjection;
73 :
74 : int parse_label(const char *file, char *keyword, char *value);
75 : int strstrip(char instr[], char outstr[], int position);
76 :
77 : CPLString oTempResult;
78 :
79 : void CleanString( CPLString &osInput );
80 :
81 : const char *GetKeyword( const char *pszPath,
82 : const char *pszDefault = "");
83 : const char *GetKeywordSub( const char *pszPath,
84 : int iSubscript,
85 : const char *pszDefault = "");
86 :
87 : public:
88 : ISIS2Dataset();
89 : ~ISIS2Dataset();
90 :
91 : virtual CPLErr GetGeoTransform( double * padfTransform );
92 : virtual const char *GetProjectionRef(void);
93 :
94 : static GDALDataset *Open( GDALOpenInfo * );
95 : static GDALDataset *Create( const char * pszFilename,
96 : int nXSize, int nYSize, int nBands,
97 : GDALDataType eType, char ** papszParmList );
98 : };
99 :
100 : /************************************************************************/
101 : /* ISIS2Dataset() */
102 : /************************************************************************/
103 :
104 1 : ISIS2Dataset::ISIS2Dataset()
105 : {
106 1 : fpImage = NULL;
107 1 : bGotTransform = FALSE;
108 1 : adfGeoTransform[0] = 0.0;
109 1 : adfGeoTransform[1] = 1.0;
110 1 : adfGeoTransform[2] = 0.0;
111 1 : adfGeoTransform[3] = 0.0;
112 1 : adfGeoTransform[4] = 0.0;
113 1 : adfGeoTransform[5] = 1.0;
114 1 : }
115 :
116 : /************************************************************************/
117 : /* ~ISIS2Dataset() */
118 : /************************************************************************/
119 :
120 2 : ISIS2Dataset::~ISIS2Dataset()
121 :
122 : {
123 1 : FlushCache();
124 1 : if( fpImage != NULL )
125 1 : VSIFCloseL( fpImage );
126 2 : }
127 :
128 : /************************************************************************/
129 : /* GetProjectionRef() */
130 : /************************************************************************/
131 :
132 1 : const char *ISIS2Dataset::GetProjectionRef()
133 :
134 : {
135 1 : if( strlen(osProjection) > 0 )
136 1 : return osProjection;
137 : else
138 0 : return GDALPamDataset::GetProjectionRef();
139 : }
140 :
141 : /************************************************************************/
142 : /* GetGeoTransform() */
143 : /************************************************************************/
144 :
145 1 : CPLErr ISIS2Dataset::GetGeoTransform( double * padfTransform )
146 :
147 : {
148 1 : if( bGotTransform )
149 : {
150 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
151 1 : return CE_None;
152 : }
153 : else
154 : {
155 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
156 : }
157 : }
158 :
159 : /************************************************************************/
160 : /* Open() */
161 : /************************************************************************/
162 :
163 9029 : GDALDataset *ISIS2Dataset::Open( GDALOpenInfo * poOpenInfo )
164 : {
165 : /* -------------------------------------------------------------------- */
166 : /* Does this look like a CUBE dataset? */
167 : /* -------------------------------------------------------------------- */
168 9029 : if( poOpenInfo->pabyHeader == NULL
169 : || strstr((const char *)poOpenInfo->pabyHeader,"^QUBE") == NULL )
170 9028 : return NULL;
171 :
172 : /* -------------------------------------------------------------------- */
173 : /* Open the file using the large file API. */
174 : /* -------------------------------------------------------------------- */
175 1 : FILE *fpQube = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
176 :
177 1 : if( fpQube == NULL )
178 0 : return NULL;
179 :
180 : ISIS2Dataset *poDS;
181 :
182 1 : poDS = new ISIS2Dataset();
183 :
184 1 : if( ! poDS->oKeywords.Ingest( fpQube, 0 ) )
185 : {
186 0 : VSIFCloseL( fpQube );
187 0 : delete poDS;
188 0 : return NULL;
189 : }
190 :
191 1 : VSIFCloseL( fpQube );
192 :
193 : /* -------------------------------------------------------------------- */
194 : /* We assume the user is pointing to the label (ie. .lab) file. */
195 : /* -------------------------------------------------------------------- */
196 : // QUBE can be inline or detached and point to an image name
197 : // ^QUBE = 76
198 : // ^QUBE = ("ui31s015.img",6441<BYTES>) - has another label on the image
199 : // ^QUBE = "ui31s015.img" - which implies no label or skip value
200 :
201 1 : const char *pszQube = poDS->GetKeyword( "^QUBE" );
202 1 : int nQube = atoi(pszQube);
203 :
204 1 : if( pszQube[0] == '"' || pszQube[0] == '(' )
205 : {
206 : CPLError( CE_Failure, CPLE_AppDefined,
207 0 : "ISIS2 driver does not support detached images." );
208 0 : return NULL;
209 : }
210 :
211 : /* -------------------------------------------------------------------- */
212 : /* Check if file an ISIS2 header file? Read a few lines of text */
213 : /* searching for something starting with nrows or ncols. */
214 : /* -------------------------------------------------------------------- */
215 1 : GDALDataType eDataType = GDT_Byte;
216 1 : OGRSpatialReference oSRS;
217 :
218 : //image parameters
219 1 : int nRows, nCols, nBands = 1;
220 1 : int nSkipBytes = 0;
221 : int itype;
222 : int s_ix, s_iy, s_iz; // check SUFFIX_ITEMS params.
223 : int record_bytes;
224 1 : int bNoDataSet = FALSE;
225 1 : char chByteOrder = 'M'; //default to MSB
226 :
227 : //Georef parameters
228 1 : double dfULXMap=0.5;
229 1 : double dfULYMap = 0.5;
230 1 : double dfXDim = 1.0;
231 1 : double dfYDim = 1.0;
232 1 : double dfNoData = 0.0;
233 1 : double xulcenter = 0.0;
234 1 : double yulcenter = 0.0;
235 :
236 : //projection parameters
237 1 : int bProjectionSet = TRUE;
238 1 : double semi_major = 0.0;
239 1 : double semi_minor = 0.0;
240 1 : double iflattening = 0.0;
241 1 : float center_lat = 0.0;
242 1 : float center_lon = 0.0;
243 1 : float first_std_parallel = 0.0;
244 1 : float second_std_parallel = 0.0;
245 : FILE *fp;
246 :
247 : /* -------------------------------------------------------------------- */
248 : /* Checks to see if this is valid ISIS2 cube */
249 : /* SUFFIX_ITEM tag in .cub file should be (0,0,0); no side-planes */
250 : /* -------------------------------------------------------------------- */
251 1 : s_ix = atoi(poDS->GetKeywordSub( "QUBE.SUFFIX_ITEMS", 1 ));
252 1 : s_iy = atoi(poDS->GetKeywordSub( "QUBE.SUFFIX_ITEMS", 2 ));
253 1 : s_iz = atoi(poDS->GetKeywordSub( "QUBE.SUFFIX_ITEMS", 3 ));
254 :
255 1 : if( s_ix != 0 || s_iy != 0 || s_iz != 0 )
256 : {
257 : CPLError( CE_Failure, CPLE_OpenFailed,
258 : "*** ISIS 2 cube file has invalid SUFFIX_ITEMS parameters:\n"
259 : "*** gdal isis2 driver requires (0, 0, 0), thus no sideplanes or backplanes\n"
260 0 : "found: (%i, %i, %i)\n\n", s_ix, s_iy, s_iz );
261 0 : return NULL;
262 : }
263 :
264 : /**************** end SUFFIX_ITEM check ***********************/
265 :
266 :
267 : /*********** Grab layout type (BSQ, BIP, BIL) ************/
268 : // AXIS_NAME = (SAMPLE,LINE,BAND)
269 : /***********************************************************/
270 : const char *value;
271 :
272 1 : char szLayout[10] = "BSQ"; //default to band seq.
273 1 : value = poDS->GetKeyword( "QUBE.AXIS_NAME", "" );
274 1 : if (EQUAL(value,"(SAMPLE,LINE,BAND)") )
275 1 : strcpy(szLayout,"BSQ");
276 0 : else if (EQUAL(value,"(BAND,LINE,SAMPLE)") )
277 0 : strcpy(szLayout,"BIP");
278 0 : else if (EQUAL(value,"(SAMPLE,BAND,LINE)") || EQUAL(value,"") )
279 0 : strcpy(szLayout,"BSQ");
280 : else {
281 : CPLError( CE_Failure, CPLE_OpenFailed,
282 0 : "%s layout not supported. Abort\n\n", value);
283 0 : return NULL;
284 : }
285 :
286 : /*********** Grab samples lines band ************/
287 1 : nCols = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS",1));
288 1 : nRows = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS",2));
289 1 : nBands = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS",3));
290 :
291 : /*********** Grab Qube record bytes **********/
292 1 : record_bytes = atoi(poDS->GetKeyword("RECORD_BYTES"));
293 :
294 1 : if (nQube > 0)
295 1 : nSkipBytes = (nQube - 1) * record_bytes;
296 : else
297 0 : nSkipBytes = 0;
298 :
299 : /******** Grab format type - isis2 only supports 8,16,32 *******/
300 1 : itype = atoi(poDS->GetKeyword("QUBE.CORE_ITEM_BYTES",""));
301 1 : switch(itype) {
302 : case 1 :
303 0 : eDataType = GDT_Byte;
304 0 : dfNoData = NULL1;
305 0 : bNoDataSet = TRUE;
306 0 : break;
307 : case 2 :
308 0 : eDataType = GDT_Int16;
309 0 : dfNoData = NULL2;
310 0 : bNoDataSet = TRUE;
311 0 : break;
312 : case 4 :
313 1 : eDataType = GDT_Float32;
314 1 : dfNoData = NULL3;
315 1 : bNoDataSet = TRUE;
316 1 : break;
317 : default :
318 : CPLError( CE_Failure, CPLE_AppDefined,
319 : "Itype of %d is not supported in ISIS 2.",
320 0 : itype);
321 0 : delete poDS;
322 0 : return NULL;
323 : }
324 :
325 : /*********** Grab samples lines band ************/
326 1 : value = poDS->GetKeyword( "QUBE.CORE_ITEM_TYPE" );
327 1 : if( (EQUAL(value,"PC_INTEGER")) ||
328 : (EQUAL(value,"PC_UNSIGNED_INTEGER")) ||
329 : (EQUAL(value,"PC_REAL")) ) {
330 0 : chByteOrder = 'I';
331 : }
332 :
333 : /*********** Grab Cellsize ************/
334 1 : value = poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.MAP_SCALE");
335 1 : if (strlen(value) > 0 ) {
336 1 : dfXDim = (float) atof(value) * 1000.0; /* convert from km to m */
337 1 : dfYDim = (float) atof(value) * 1000.0 * -1;
338 : }
339 :
340 : /*********** Grab LINE_PROJECTION_OFFSET ************/
341 1 : value = poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.LINE_PROJECTION_OFFSET");
342 1 : if (strlen(value) > 0) {
343 1 : yulcenter = (float) atof(value);
344 1 : yulcenter = ((yulcenter) * dfYDim);
345 1 : dfULYMap = yulcenter - (dfYDim/2);
346 : }
347 :
348 : /*********** Grab SAMPLE_PROJECTION_OFFSET ************/
349 1 : value = poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.SAMPLE_PROJECTION_OFFSET");
350 1 : if( strlen(value) > 0 ) {
351 1 : xulcenter = (float) atof(value);
352 1 : xulcenter = ((xulcenter) * dfXDim);
353 1 : dfULXMap = xulcenter - (dfXDim/2);
354 : }
355 :
356 : /*********** Grab TARGET_NAME ************/
357 : /**** This is the planets name i.e. MARS ***/
358 1 : CPLString target_name = poDS->GetKeyword("QUBE.TARGET_NAME");
359 :
360 : /*********** Grab MAP_PROJECTION_TYPE ************/
361 : CPLString map_proj_name =
362 1 : poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.MAP_PROJECTION_TYPE");
363 1 : poDS->CleanString( map_proj_name );
364 :
365 : /*********** Grab SEMI-MAJOR ************/
366 : semi_major =
367 1 : atof(poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.A_AXIS_RADIUS")) * 1000.0;
368 :
369 : /*********** Grab semi-minor ************/
370 : semi_minor =
371 1 : atof(poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.C_AXIS_RADIUS")) * 1000.0;
372 :
373 : /*********** Grab CENTER_LAT ************/
374 : center_lat =
375 1 : atof(poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.CENTER_LATITUDE"));
376 :
377 : /*********** Grab CENTER_LON ************/
378 : center_lon =
379 1 : atof(poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.CENTER_LONGITUDE"));
380 :
381 : /*********** Grab 1st std parallel ************/
382 : first_std_parallel =
383 1 : atof(poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.FIRST_STANDARD_PARALLEL"));
384 :
385 : /*********** Grab 2nd std parallel ************/
386 : second_std_parallel =
387 1 : atof(poDS->GetKeyword( "QUBE.IMAGE_MAP_PROJECTION.SECOND_STANDARD_PARALLEL"));
388 :
389 : /*** grab PROJECTION_LATITUDE_TYPE = "PLANETOCENTRIC" ****/
390 : // Need to further study how ocentric/ographic will effect the gdal library.
391 : // So far we will use this fact to define a sphere or ellipse for some projections
392 : // Frank - may need to talk this over
393 1 : char bIsGeographic = TRUE;
394 1 : value = poDS->GetKeyword("CUBE.IMAGE_MAP_PROJECTION.PROJECTION_LATITUDE_TYPE");
395 1 : if (EQUAL( value, "\"PLANETOCENTRIC\"" ))
396 0 : bIsGeographic = FALSE;
397 :
398 1 : CPLDebug("ISIS2","using projection %s", map_proj_name.c_str() );
399 :
400 : //Set oSRS projection and parameters
401 1 : if ((EQUAL( map_proj_name, "EQUIRECTANGULAR_CYLINDRICAL" )) ||
402 : (EQUAL( map_proj_name, "EQUIRECTANGULAR" )) ||
403 : (EQUAL( map_proj_name, "SIMPLE_CYLINDRICAL" )) ) {
404 1 : oSRS.OGRSpatialReference::SetEquirectangular2 ( 0.0, center_lon, center_lat, 0, 0 );
405 0 : } else if (EQUAL( map_proj_name, "ORTHOGRAPHIC" )) {
406 0 : oSRS.OGRSpatialReference::SetOrthographic ( center_lat, center_lon, 0, 0 );
407 0 : } else if ((EQUAL( map_proj_name, "SINUSOIDAL" )) ||
408 : (EQUAL( map_proj_name, "SINUSOIDAL_EQUAL-AREA" ))) {
409 0 : oSRS.OGRSpatialReference::SetSinusoidal ( center_lon, 0, 0 );
410 0 : } else if (EQUAL( map_proj_name, "MERCATOR" )) {
411 0 : oSRS.OGRSpatialReference::SetMercator ( center_lat, center_lon, 1, 0, 0 );
412 0 : } else if (EQUAL( map_proj_name, "POLAR_STEREOGRAPHIC" )) {
413 0 : oSRS.OGRSpatialReference::SetPS ( center_lat, center_lon, 1, 0, 0 );
414 0 : } else if (EQUAL( map_proj_name, "TRANSVERSE_MERCATOR" )) {
415 0 : oSRS.OGRSpatialReference::SetTM ( center_lat, center_lon, 1, 0, 0 );
416 0 : } else if (EQUAL( map_proj_name, "LAMBERT_CONFORMAL_CONIC" )) {
417 0 : oSRS.OGRSpatialReference::SetLCC ( first_std_parallel, second_std_parallel, center_lat, center_lon, 0, 0 );
418 : } else {
419 : CPLDebug( "ISIS2",
420 : "Dataset projection %s is not supported. Continuing...",
421 0 : map_proj_name.c_str() );
422 0 : bProjectionSet = FALSE;
423 : }
424 :
425 1 : if (bProjectionSet) {
426 : //Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
427 1 : CPLString proj_target_name = map_proj_name + " " + target_name;
428 1 : oSRS.SetProjCS(proj_target_name); //set ProjCS keyword
429 :
430 : //The geographic/geocentric name will be the same basic name as the body name
431 : //'GCS' = Geographic/Geocentric Coordinate System
432 1 : CPLString geog_name = "GCS_" + target_name;
433 :
434 : //The datum and sphere names will be the same basic name aas the planet
435 1 : CPLString datum_name = "D_" + target_name;
436 1 : CPLString sphere_name = target_name; // + "_IAU_IAG"); //Might not be IAU defined so don't add
437 :
438 : //calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
439 1 : if ((semi_major - semi_minor) < 0.0000001)
440 1 : iflattening = 0;
441 : else
442 0 : iflattening = semi_major / (semi_major - semi_minor);
443 :
444 : //Set the body size but take into consideration which proj is being used to help w/ proj4 compatibility
445 : //The use of a Sphere, polar radius or ellipse here is based on how ISIS does it internally
446 1 : if ( ( (EQUAL( map_proj_name, "STEREOGRAPHIC" ) && (fabs(center_lat) == 90)) ) ||
447 : (EQUAL( map_proj_name, "POLAR_STEREOGRAPHIC" )))
448 : {
449 0 : if (bIsGeographic) {
450 : //Geograpraphic, so set an ellipse
451 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
452 : semi_major, iflattening,
453 0 : "Reference_Meridian", 0.0 );
454 : } else {
455 : //Geocentric, so force a sphere using the semi-minor axis. I hope...
456 0 : sphere_name += "_polarRadius";
457 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
458 : semi_minor, 0.0,
459 0 : "Reference_Meridian", 0.0 );
460 : }
461 : }
462 1 : else if ( (EQUAL( map_proj_name, "SIMPLE_CYLINDRICAL" )) ||
463 : (EQUAL( map_proj_name, "ORTHOGRAPHIC" )) ||
464 : (EQUAL( map_proj_name, "STEREOGRAPHIC" )) ||
465 : (EQUAL( map_proj_name, "SINUSOIDAL_EQUAL-AREA" )) ||
466 : (EQUAL( map_proj_name, "SINUSOIDAL" )) ) {
467 : //isis uses the sphereical equation for these projections so force a sphere
468 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
469 : semi_major, 0.0,
470 1 : "Reference_Meridian", 0.0 );
471 : }
472 0 : else if ((EQUAL( map_proj_name, "EQUIRECTANGULAR_CYLINDRICAL" )) ||
473 : (EQUAL( map_proj_name, "EQUIRECTANGULAR" )) ) {
474 : //Calculate localRadius using ISIS3 simple elliptical method
475 : // not the more standard Radius of Curvature method
476 : //PI = 4 * atan(1);
477 : double radLat, localRadius;
478 0 : if (center_lon == 0) { //No need to calculate local radius
479 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
480 : semi_major, 0.0,
481 0 : "Reference_Meridian", 0.0 );
482 : } else {
483 0 : radLat = center_lat * PI / 180; // in radians
484 : localRadius = semi_major * semi_minor / sqrt(pow(semi_minor*cos(radLat),2)
485 0 : + pow(semi_major*sin(radLat),2) );
486 0 : sphere_name += "_localRadius";
487 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
488 : localRadius, 0.0,
489 0 : "Reference_Meridian", 0.0 );
490 0 : CPLDebug( "ISIS2", "local radius: %f", localRadius);
491 : }
492 : }
493 : else {
494 : //All other projections: Mercator, Transverse Mercator, Lambert Conformal, etc.
495 : //Geographic, so set an ellipse
496 0 : if (bIsGeographic) {
497 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
498 : semi_major, iflattening,
499 0 : "Reference_Meridian", 0.0 );
500 : } else {
501 : //Geocentric, so force a sphere. I hope...
502 : oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
503 : semi_major, 0.0,
504 0 : "Reference_Meridian", 0.0 );
505 : }
506 : }
507 :
508 :
509 : // translate back into a projection string.
510 1 : char *pszResult = NULL;
511 1 : oSRS.exportToWkt( &pszResult );
512 1 : poDS->osProjection = pszResult;
513 1 : CPLFree( pszResult );
514 : }
515 :
516 : /* END ISIS2 Label Read */
517 : /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
518 :
519 : /* -------------------------------------------------------------------- */
520 : /* Did we get the required keywords? If not we return with */
521 : /* this never having been considered to be a match. This isn't */
522 : /* an error! */
523 : /* -------------------------------------------------------------------- */
524 1 : if( nRows < 1 || nCols < 1 || nBands < 1 )
525 : {
526 0 : return NULL;
527 : }
528 :
529 : /* -------------------------------------------------------------------- */
530 : /* Capture some information from the file that is of interest. */
531 : /* -------------------------------------------------------------------- */
532 1 : poDS->nRasterXSize = nCols;
533 1 : poDS->nRasterYSize = nRows;
534 :
535 : /* -------------------------------------------------------------------- */
536 : /* Open target binary file. */
537 : /* -------------------------------------------------------------------- */
538 :
539 1 : if( poOpenInfo->eAccess == GA_ReadOnly )
540 1 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
541 : else
542 0 : poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "r+b" );
543 :
544 1 : if( poDS->fpImage == NULL )
545 : {
546 : CPLError( CE_Failure, CPLE_OpenFailed,
547 : "Failed to open %s with write permission.\n%s",
548 0 : poOpenInfo->pszFilename, VSIStrerror( errno ) );
549 0 : delete poDS;
550 0 : return NULL;
551 : }
552 :
553 1 : poDS->eAccess = poOpenInfo->eAccess;
554 :
555 : /* -------------------------------------------------------------------- */
556 : /* Compute the line offset. */
557 : /* -------------------------------------------------------------------- */
558 1 : int nItemSize = GDALGetDataTypeSize(eDataType)/8;
559 : int nLineOffset, nPixelOffset, nBandOffset;
560 :
561 1 : if( EQUAL(szLayout,"BIP") )
562 : {
563 0 : nPixelOffset = nItemSize * nBands;
564 0 : nLineOffset = nPixelOffset * nCols;
565 0 : nBandOffset = nItemSize;
566 : }
567 1 : else if( EQUAL(szLayout,"BSQ") )
568 : {
569 1 : nPixelOffset = nItemSize;
570 1 : nLineOffset = nPixelOffset * nCols;
571 1 : nBandOffset = nLineOffset * nRows;
572 : }
573 : else /* assume BIL */
574 : {
575 0 : nPixelOffset = nItemSize;
576 0 : nLineOffset = nItemSize * nBands * nCols;
577 0 : nBandOffset = nItemSize * nCols;
578 : }
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* Create band information objects. */
582 : /* -------------------------------------------------------------------- */
583 : int i;
584 :
585 1 : poDS->nBands = nBands;;
586 2 : for( i = 0; i < poDS->nBands; i++ )
587 : {
588 : RawRasterBand *poBand;
589 :
590 : poBand =
591 : new RawRasterBand( poDS, i+1, poDS->fpImage,
592 : nSkipBytes + nBandOffset * i,
593 : nPixelOffset, nLineOffset, eDataType,
594 : #ifdef CPL_LSB
595 : chByteOrder == 'I' || chByteOrder == 'L',
596 : #else
597 : chByteOrder == 'M',
598 : #endif
599 1 : TRUE );
600 :
601 1 : if( bNoDataSet )
602 1 : poBand->SetNoDataValue( dfNoData );
603 :
604 1 : poDS->SetBand( i+1, poBand );
605 :
606 : // Set offset/scale values at the PAM level.
607 : poBand->SetOffset(
608 1 : CPLAtofM(poDS->GetKeyword("QUBE.CORE_BASE","0.0")));
609 : poBand->SetScale(
610 1 : CPLAtofM(poDS->GetKeyword("QUBE.CORE_MULTIPLIER","1.0")));
611 : }
612 :
613 : /* -------------------------------------------------------------------- */
614 : /* Check for a .prj file. For isis2 I would like to keep this in */
615 : /* -------------------------------------------------------------------- */
616 1 : CPLString osPath, osName;
617 :
618 1 : osPath = CPLGetPath( poOpenInfo->pszFilename );
619 1 : osName = CPLGetBasename(poOpenInfo->pszFilename);
620 1 : const char *pszPrjFile = CPLFormCIFilename( osPath, osName, "prj" );
621 :
622 1 : fp = VSIFOpen( pszPrjFile, "r" );
623 1 : if( fp != NULL )
624 : {
625 : char **papszLines;
626 0 : OGRSpatialReference oSRS;
627 :
628 0 : VSIFClose( fp );
629 :
630 0 : papszLines = CSLLoad( pszPrjFile );
631 :
632 0 : if( oSRS.importFromESRI( papszLines ) == OGRERR_NONE )
633 : {
634 0 : char *pszResult = NULL;
635 0 : oSRS.exportToWkt( &pszResult );
636 0 : poDS->osProjection = pszResult;
637 0 : CPLFree( pszResult );
638 : }
639 :
640 0 : CSLDestroy( papszLines );
641 : }
642 :
643 :
644 1 : if( dfULYMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0 )
645 : {
646 1 : poDS->bGotTransform = TRUE;
647 1 : poDS->adfGeoTransform[0] = dfULXMap;
648 1 : poDS->adfGeoTransform[1] = dfXDim;
649 1 : poDS->adfGeoTransform[2] = 0.0;
650 1 : poDS->adfGeoTransform[3] = dfULYMap;
651 1 : poDS->adfGeoTransform[4] = 0.0;
652 1 : poDS->adfGeoTransform[5] = dfYDim;
653 : }
654 :
655 1 : if( !poDS->bGotTransform )
656 : poDS->bGotTransform =
657 : GDALReadWorldFile( poOpenInfo->pszFilename, "cbw",
658 0 : poDS->adfGeoTransform );
659 :
660 1 : if( !poDS->bGotTransform )
661 : poDS->bGotTransform =
662 : GDALReadWorldFile( poOpenInfo->pszFilename, "wld",
663 0 : poDS->adfGeoTransform );
664 :
665 : /* -------------------------------------------------------------------- */
666 : /* Initialize any PAM information. */
667 : /* -------------------------------------------------------------------- */
668 1 : poDS->SetDescription( poOpenInfo->pszFilename );
669 1 : poDS->TryLoadXML();
670 :
671 : /* -------------------------------------------------------------------- */
672 : /* Check for overviews. */
673 : /* -------------------------------------------------------------------- */
674 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
675 :
676 1 : return( poDS );
677 : }
678 :
679 : /************************************************************************/
680 : /* GetKeyword() */
681 : /************************************************************************/
682 :
683 19 : const char *ISIS2Dataset::GetKeyword( const char *pszPath,
684 : const char *pszDefault )
685 :
686 : {
687 19 : return oKeywords.GetKeyword( pszPath, pszDefault );
688 : }
689 :
690 : /************************************************************************/
691 : /* GetKeywordSub() */
692 : /************************************************************************/
693 :
694 6 : const char *ISIS2Dataset::GetKeywordSub( const char *pszPath,
695 : int iSubscript,
696 : const char *pszDefault )
697 :
698 : {
699 6 : const char *pszResult = oKeywords.GetKeyword( pszPath, NULL );
700 :
701 6 : if( pszResult == NULL )
702 0 : return pszDefault;
703 :
704 6 : if( pszResult[0] != '(' )
705 0 : return pszDefault;
706 :
707 : char **papszTokens = CSLTokenizeString2( pszResult, "(,)",
708 6 : CSLT_HONOURSTRINGS );
709 :
710 6 : if( iSubscript <= CSLCount(papszTokens) )
711 : {
712 6 : oTempResult = papszTokens[iSubscript-1];
713 6 : CSLDestroy( papszTokens );
714 6 : return oTempResult.c_str();
715 : }
716 : else
717 : {
718 0 : CSLDestroy( papszTokens );
719 0 : return pszDefault;
720 : }
721 : }
722 :
723 : /************************************************************************/
724 : /* CleanString() */
725 : /* */
726 : /* Removes single or double quotes, and converts spaces to underscores. */
727 : /* The change is made in-place to CPLString. */
728 : /************************************************************************/
729 :
730 1 : void ISIS2Dataset::CleanString( CPLString &osInput )
731 :
732 : {
733 1 : if( ( osInput.size() < 2 ) ||
734 : ((osInput.at(0) != '"' || osInput.at(osInput.size()-1) != '"' ) &&
735 : ( osInput.at(0) != '\'' || osInput.at(osInput.size()-1) != '\'')) )
736 1 : return;
737 :
738 0 : char *pszWrk = CPLStrdup(osInput.c_str() + 1);
739 : int i;
740 :
741 0 : pszWrk[strlen(pszWrk)-1] = '\0';
742 :
743 0 : for( i = 0; pszWrk[i] != '\0'; i++ )
744 : {
745 0 : if( pszWrk[i] == ' ' )
746 0 : pszWrk[i] = '_';
747 : }
748 :
749 0 : osInput = pszWrk;
750 0 : CPLFree( pszWrk );
751 : }
752 :
753 : /************************************************************************/
754 : /* GDALRegister_ISIS2() */
755 : /************************************************************************/
756 :
757 338 : void GDALRegister_ISIS2()
758 :
759 : {
760 : GDALDriver *poDriver;
761 :
762 338 : if( GDALGetDriverByName( "ISIS2" ) == NULL )
763 : {
764 336 : poDriver = new GDALDriver();
765 :
766 336 : poDriver->SetDescription( "ISIS2" );
767 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
768 336 : "USGS Astrogeology ISIS cube (Version 2)" );
769 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
770 336 : "frmt_various.html#ISIS2" );
771 :
772 336 : poDriver->pfnOpen = ISIS2Dataset::Open;
773 :
774 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
775 : }
776 338 : }
|
