LCOV - code coverage report
Current view: directory - frmts/pds - pdsdataset.cpp (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 447 351 78.5 %
Date: 2011-12-18 Functions: 25 19 76.0 %

       1                 : /******************************************************************************
       2                 :  * $Id: pdsdataset.cpp 23471 2011-12-05 20:37:37Z rouault $
       3                 :  *
       4                 :  * Project:  PDS Driver; Planetary Data System Format
       5                 :  * Purpose:  Implementation of PDSDataset
       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) 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                 : // Set up PDS NULL values
      36                 : #define NULL1 0
      37                 : #define NULL2 -32768
      38                 : //#define NULL3 -0.3402822655089E+39
      39                 : //Same as ESRI_GRID_FLOAT_NO_DATA
      40                 : //#define NULL3 -340282346638528859811704183484516925440.0
      41                 : #define NULL3 -3.4028226550889044521e+38
      42                 : 
      43                 : #include "rawdataset.h"
      44                 : #include "gdal_proxy.h"
      45                 : #include "ogr_spatialref.h"
      46                 : #include "cpl_string.h" 
      47                 : #include "nasakeywordhandler.h"
      48                 : 
      49                 : CPL_CVSID("$Id: pdsdataset.cpp 23471 2011-12-05 20:37:37Z rouault $");
      50                 : 
      51                 : CPL_C_START
      52                 : void  GDALRegister_PDS(void);
      53                 : CPL_C_END
      54                 : 
      55                 : /************************************************************************/
      56                 : /* ==================================================================== */
      57                 : /*             PDSDataset                         */
      58                 : /* ==================================================================== */
      59                 : /************************************************************************/
      60                 : 
      61                 : class PDSDataset : public RawDataset
      62                 : {
      63                 :     VSILFILE  *fpImage; // image data file.
      64                 :     GDALDataset *poCompressedDS;
      65                 : 
      66                 :     NASAKeywordHandler  oKeywords;
      67                 :   
      68                 :     int         bGotTransform;
      69                 :     double      adfGeoTransform[6];
      70                 :   
      71                 :     CPLString   osProjection;
      72                 : 
      73                 :     CPLString   osTempResult;
      74                 : 
      75                 :     void        ParseSRS();
      76                 :     int         ParseCompressedImage();
      77                 :     int         ParseImage( CPLString osPrefix = "" );
      78                 :     void        CleanString( CPLString &osInput );
      79                 : 
      80                 :     const char *GetKeyword( std::string osPath,
      81                 :                             const char *pszDefault = "");
      82                 :     const char *GetKeywordSub( std::string osPath,
      83                 :                                int iSubscript, 
      84                 :                                const char *pszDefault = "");
      85                 :     const char *GetKeywordUnit( const char *pszPath, 
      86                 :                                int iSubscript, 
      87                 :                                const char *pszDefault = "");
      88                 : 
      89                 :   protected:
      90                 :     virtual int         CloseDependentDatasets();
      91                 : 
      92                 : public:
      93                 :     PDSDataset();
      94                 :     ~PDSDataset();
      95                 :   
      96                 :     virtual CPLErr GetGeoTransform( double * padfTransform );
      97                 :     virtual const char *GetProjectionRef(void);
      98                 :   
      99                 :     virtual char      **GetFileList(void);
     100                 : 
     101                 :     virtual CPLErr IBuildOverviews( const char *, int, int *,
     102                 :                                     int, int *, GDALProgressFunc, void * );
     103                 :     
     104                 :     virtual CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
     105                 :                               void *, int, int, GDALDataType,
     106                 :                               int, int *, int, int, int );
     107                 : 
     108                 :     static int          Identify( GDALOpenInfo * );
     109                 :     static GDALDataset *Open( GDALOpenInfo * );
     110                 :     static GDALDataset *Create( const char * pszFilename,
     111                 :                                 int nXSize, int nYSize, int nBands,
     112                 :                                 GDALDataType eType, char ** papszParmList );
     113                 : };
     114                 : 
     115                 : /************************************************************************/
     116                 : /*                            PDSDataset()                            */
     117                 : /************************************************************************/
     118                 : 
     119              14 : PDSDataset::PDSDataset()
     120                 : {
     121              14 :     fpImage = NULL;
     122              14 :     bGotTransform = FALSE;
     123              14 :     adfGeoTransform[0] = 0.0;
     124              14 :     adfGeoTransform[1] = 1.0;
     125              14 :     adfGeoTransform[2] = 0.0;
     126              14 :     adfGeoTransform[3] = 0.0;
     127              14 :     adfGeoTransform[4] = 0.0;
     128              14 :     adfGeoTransform[5] = 1.0;
     129              14 :     poCompressedDS = NULL;
     130              14 : }
     131                 : 
     132                 : /************************************************************************/
     133                 : /*                            ~PDSDataset()                            */
     134                 : /************************************************************************/
     135                 : 
     136              14 : PDSDataset::~PDSDataset()
     137                 : 
     138                 : {
     139              14 :     FlushCache();
     140              14 :     if( fpImage != NULL )
     141              10 :         VSIFCloseL( fpImage );
     142                 : 
     143              14 :     CloseDependentDatasets();
     144              14 : }
     145                 : 
     146                 : /************************************************************************/
     147                 : /*                        CloseDependentDatasets()                      */
     148                 : /************************************************************************/
     149                 : 
     150              14 : int PDSDataset::CloseDependentDatasets()
     151                 : {
     152              14 :     int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
     153                 : 
     154              14 :     if( poCompressedDS )
     155                 :     {
     156               2 :         bHasDroppedRef = TRUE;
     157               2 :         delete poCompressedDS;
     158               2 :         poCompressedDS = NULL;
     159                 :     }
     160                 : 
     161              26 :     for( int iBand = 0; iBand < nBands; iBand++ )
     162                 :     {
     163              12 :        delete papoBands[iBand];
     164                 :     }
     165              14 :     nBands = 0;
     166                 : 
     167              14 :     return bHasDroppedRef;
     168                 : }
     169                 : 
     170                 : /************************************************************************/
     171                 : /*                            GetFileList()                             */
     172                 : /************************************************************************/
     173                 : 
     174               1 : char **PDSDataset::GetFileList()
     175                 : 
     176                 : {
     177               1 :     char **papszFileList = RawDataset::GetFileList();
     178                 : 
     179               1 :     if( poCompressedDS != NULL )
     180                 :     {
     181               1 :         char **papszCFileList = poCompressedDS->GetFileList();
     182                 : 
     183                 :         papszFileList = CSLInsertStrings( papszFileList, -1, 
     184               1 :                                           papszCFileList );
     185               1 :         CSLDestroy( papszCFileList );
     186                 :     }
     187                 : 
     188               1 :     return papszFileList;
     189                 : }
     190                 : 
     191                 : /************************************************************************/
     192                 : /*                          IBuildOverviews()                           */
     193                 : /************************************************************************/
     194                 : 
     195               0 : CPLErr PDSDataset::IBuildOverviews( const char *pszResampling, 
     196                 :                                     int nOverviews, int *panOverviewList, 
     197                 :                                     int nListBands, int *panBandList,
     198                 :                                     GDALProgressFunc pfnProgress, 
     199                 :                                     void * pProgressData )
     200                 : {
     201               0 :     if( poCompressedDS != NULL )
     202                 :         return poCompressedDS->BuildOverviews( pszResampling, 
     203                 :                                                nOverviews, panOverviewList,
     204                 :                                                nListBands, panBandList, 
     205               0 :                                                pfnProgress, pProgressData );
     206                 :     else
     207                 :         return RawDataset::IBuildOverviews( pszResampling, 
     208                 :                                             nOverviews, panOverviewList,
     209                 :                                             nListBands, panBandList, 
     210               0 :                                             pfnProgress, pProgressData );
     211                 : }
     212                 :         
     213                 : /************************************************************************/
     214                 : /*                             IRasterIO()                              */
     215                 : /************************************************************************/
     216                 : 
     217               0 : CPLErr PDSDataset::IRasterIO( GDALRWFlag eRWFlag,
     218                 :                               int nXOff, int nYOff, int nXSize, int nYSize,
     219                 :                               void * pData, int nBufXSize, int nBufYSize,
     220                 :                               GDALDataType eBufType, 
     221                 :                               int nBandCount, int *panBandMap,
     222                 :                               int nPixelSpace, int nLineSpace, int nBandSpace)
     223                 : 
     224                 : {
     225               0 :     if( poCompressedDS != NULL )
     226                 :         return poCompressedDS->RasterIO( eRWFlag, 
     227                 :                                          nXOff, nYOff, nXSize, nYSize, 
     228                 :                                          pData, nBufXSize, nBufYSize, 
     229                 :                                          eBufType, nBandCount, panBandMap,
     230               0 :                                          nPixelSpace, nLineSpace, nBandSpace );
     231                 :     else
     232                 :         return RawDataset::IRasterIO( eRWFlag, 
     233                 :                                       nXOff, nYOff, nXSize, nYSize, 
     234                 :                                       pData, nBufXSize, nBufYSize, 
     235                 :                                       eBufType, nBandCount, panBandMap,
     236               0 :                                       nPixelSpace, nLineSpace, nBandSpace );
     237                 : }
     238                 : 
     239                 : /************************************************************************/
     240                 : /*                          GetProjectionRef()                          */
     241                 : /************************************************************************/
     242                 : 
     243               4 : const char *PDSDataset::GetProjectionRef()
     244                 : 
     245                 : {
     246               4 :     if( strlen(osProjection) > 0 )
     247               4 :         return osProjection;
     248                 :     else
     249               0 :         return GDALPamDataset::GetProjectionRef();
     250                 : }
     251                 : 
     252                 : /************************************************************************/
     253                 : /*                          GetGeoTransform()                           */
     254                 : /************************************************************************/
     255                 : 
     256               7 : CPLErr PDSDataset::GetGeoTransform( double * padfTransform )
     257                 : 
     258                 : {
     259               7 :     if( bGotTransform )
     260                 :     {
     261               6 :         memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
     262               6 :         return CE_None;
     263                 :     }
     264                 :     else
     265                 :     {
     266               1 :         return GDALPamDataset::GetGeoTransform( padfTransform );
     267                 :     }
     268                 : }
     269                 : 
     270                 : /************************************************************************/
     271                 : /*                              ParseSRS()                              */
     272                 : /************************************************************************/
     273                 : 
     274              12 : void PDSDataset::ParseSRS()
     275                 : 
     276                 : {
     277              12 :     const char *pszFilename = GetDescription();
     278                 : 
     279                 : /* ==================================================================== */
     280                 : /*      Get the geotransform.                                           */
     281                 : /* ==================================================================== */
     282                 :     /***********   Grab Cellsize ************/
     283                 :     //example:
     284                 :     //MAP_SCALE   = 14.818 <KM/PIXEL>
     285                 :     //added search for unit (only checks for CM, KM - defaults to Meters)
     286                 :     const char *value;
     287                 :     //Georef parameters
     288              12 :     double dfULXMap=0.5;
     289              12 :     double dfULYMap = 0.5;
     290              12 :     double dfXDim = 1.0;
     291              12 :     double dfYDim = 1.0;
     292              12 :     double xulcenter = 0.0;
     293              12 :     double yulcenter = 0.0;
     294                 : 
     295              12 :     value = GetKeyword("IMAGE_MAP_PROJECTION.MAP_SCALE");
     296              12 :     if (strlen(value) > 0 ) {
     297              10 :         dfXDim = atof(value);
     298              10 :         dfYDim = atof(value) * -1;
     299                 :         
     300              10 :         CPLString unit = GetKeywordUnit("IMAGE_MAP_PROJECTION.MAP_SCALE",2); //KM
     301                 :         //value = GetKeywordUnit("IMAGE_MAP_PROJECTION.MAP_SCALE",3); //PIXEL
     302              10 :         if((EQUAL(unit,"M"))  || (EQUAL(unit,"METER")) || (EQUAL(unit,"METERS"))) {
     303                 :             // do nothing
     304                 :         }
     305               6 :         else if (EQUAL(unit,"CM")) {
     306                 :             // convert from cm to m
     307               0 :             dfXDim = dfXDim / 100.0;
     308               0 :             dfYDim = dfYDim / 100.0;
     309                 :         } else {
     310                 :             //defaults to convert km to m
     311               6 :             dfXDim = dfXDim * 1000.0;
     312               6 :             dfYDim = dfYDim * 1000.0;
     313              10 :         }            
     314                 :     }
     315                 :     
     316                 : /* -------------------------------------------------------------------- */
     317                 : /*      Calculate upper left corner of pixel in meters from the         */
     318                 : /*      upper  left center pixel sample/line offsets.  It doesn't       */
     319                 : /*      mean the defaults will work for every PDS image, as these       */
     320                 : /*      values are used inconsistantly.  Thus we have included          */
     321                 : /*      conversion options to allow the user to override the            */
     322                 : /*      documented PDS3 default. Jan. 2011, for known mapping issues    */
     323                 : /*      see GDAL PDS page or mapping within ISIS3 source (USGS)         */
     324                 : /*      $ISIS3DATA/base/translations/pdsProjectionLineSampToXY.def      */
     325                 : /* -------------------------------------------------------------------- */
     326                 :    
     327                 :     // defaults should be correct for what is documented in the PDS3 standard
     328                 :     double   dfSampleOffset_Shift;
     329                 :     double   dfLineOffset_Shift;
     330                 :     double   dfSampleOffset_Mult;
     331                 :     double   dfLineOffset_Mult;
     332                 : 
     333                 :     dfSampleOffset_Shift = 
     334              12 :         atof(CPLGetConfigOption( "PDS_SampleProjOffset_Shift", "-0.5" ));
     335                 :     
     336                 :     dfLineOffset_Shift = 
     337              12 :         atof(CPLGetConfigOption( "PDS_LineProjOffset_Shift", "-0.5" ));
     338                 : 
     339                 :     dfSampleOffset_Mult =
     340              12 :         atof(CPLGetConfigOption( "PDS_SampleProjOffset_Mult", "-1.0") );
     341                 : 
     342                 :     dfLineOffset_Mult = 
     343              12 :         atof( CPLGetConfigOption( "PDS_LineProjOffset_Mult", "1.0") );
     344                 : 
     345                 :     /***********   Grab LINE_PROJECTION_OFFSET ************/
     346              12 :     value = GetKeyword("IMAGE_MAP_PROJECTION.LINE_PROJECTION_OFFSET");
     347              12 :     if (strlen(value) > 0) {
     348              10 :         yulcenter = atof(value);
     349              10 :         dfULYMap = ((yulcenter + dfLineOffset_Shift) * -dfYDim * dfLineOffset_Mult);
     350                 :         //notice dfYDim is negative here which is why it is again negated here
     351                 :     }
     352                 :     /***********   Grab SAMPLE_PROJECTION_OFFSET ************/
     353              12 :     value = GetKeyword("IMAGE_MAP_PROJECTION.SAMPLE_PROJECTION_OFFSET");
     354              12 :     if( strlen(value) > 0 ) {
     355              10 :         xulcenter = atof(value);
     356              10 :         dfULXMap = ((xulcenter + dfSampleOffset_Shift) * dfXDim * dfSampleOffset_Mult);
     357                 :     }
     358                 : 
     359                 : /* ==================================================================== */
     360                 : /*      Get the coordinate system.                                      */
     361                 : /* ==================================================================== */
     362              12 :     int bProjectionSet = TRUE;
     363              12 :     double semi_major = 0.0;
     364              12 :     double semi_minor = 0.0;
     365              12 :     double iflattening = 0.0;
     366              12 :     double center_lat = 0.0;
     367              12 :     double center_lon = 0.0;
     368              12 :     double first_std_parallel = 0.0;
     369              12 :     double second_std_parallel = 0.0;
     370              12 :     OGRSpatialReference oSRS;
     371                 : 
     372                 :     /***********  Grab TARGET_NAME  ************/
     373                 :     /**** This is the planets name i.e. MARS ***/
     374              12 :     CPLString target_name = GetKeyword("TARGET_NAME");
     375              12 :     CleanString( target_name );
     376                 :      
     377                 :     /**********   Grab MAP_PROJECTION_TYPE *****/
     378                 :     CPLString map_proj_name = 
     379              12 :         GetKeyword( "IMAGE_MAP_PROJECTION.MAP_PROJECTION_TYPE");
     380              12 :     CleanString( map_proj_name );
     381                 :      
     382                 :     /******  Grab semi_major & convert to KM ******/
     383                 :     semi_major = 
     384              12 :         atof(GetKeyword( "IMAGE_MAP_PROJECTION.A_AXIS_RADIUS")) * 1000.0;
     385                 :     
     386                 :     /******  Grab semi-minor & convert to KM ******/
     387                 :     semi_minor = 
     388              24 :         atof(GetKeyword( "IMAGE_MAP_PROJECTION.C_AXIS_RADIUS")) * 1000.0;
     389                 : 
     390                 :     /***********   Grab CENTER_LAT ************/
     391                 :     center_lat = 
     392              24 :         atof(GetKeyword( "IMAGE_MAP_PROJECTION.CENTER_LATITUDE"));
     393                 : 
     394                 :     /***********   Grab CENTER_LON ************/
     395                 :     center_lon = 
     396              24 :         atof(GetKeyword( "IMAGE_MAP_PROJECTION.CENTER_LONGITUDE"));
     397                 : 
     398                 :     /**********   Grab 1st std parallel *******/
     399                 :     first_std_parallel = 
     400              24 :         atof(GetKeyword( "IMAGE_MAP_PROJECTION.FIRST_STANDARD_PARALLEL"));
     401                 : 
     402                 :     /**********   Grab 2nd std parallel *******/
     403                 :     second_std_parallel = 
     404              24 :         atof(GetKeyword( "IMAGE_MAP_PROJECTION.SECOND_STANDARD_PARALLEL"));
     405                 :      
     406                 :     /*** grab  PROJECTION_LATITUDE_TYPE = "PLANETOCENTRIC" ****/
     407                 :     // Need to further study how ocentric/ographic will effect the gdal library.
     408                 :     // So far we will use this fact to define a sphere or ellipse for some projections
     409                 :     // Frank - may need to talk this over
     410              12 :     char bIsGeographic = TRUE;
     411              24 :     value = GetKeyword("IMAGE_MAP_PROJECTION.COORDINATE_SYSTEM_NAME");
     412              12 :     if (EQUAL( value, "PLANETOCENTRIC" ))
     413               4 :         bIsGeographic = FALSE; 
     414                 : 
     415                 : /**   Set oSRS projection and parameters --- all PDS supported types added if apparently supported in oSRS
     416                 :       "AITOFF",  ** Not supported in GDAL??
     417                 :       "ALBERS", 
     418                 :       "BONNE",
     419                 :       "BRIESEMEISTER",   ** Not supported in GDAL??
     420                 :       "CYLINDRICAL EQUAL AREA",
     421                 :       "EQUIDISTANT",
     422                 :       "EQUIRECTANGULAR",
     423                 :       "GNOMONIC",
     424                 :       "HAMMER",    ** Not supported in GDAL??
     425                 :       "HENDU",     ** Not supported in GDAL??
     426                 :       "LAMBERT AZIMUTHAL EQUAL AREA",
     427                 :       "LAMBERT CONFORMAL",
     428                 :       "MERCATOR",
     429                 :       "MOLLWEIDE",
     430                 :       "OBLIQUE CYLINDRICAL",
     431                 :       "ORTHOGRAPHIC",
     432                 :       "SIMPLE CYLINDRICAL",
     433                 :       "SINUSOIDAL",
     434                 :       "STEREOGRAPHIC",
     435                 :       "TRANSVERSE MERCATOR",
     436                 :       "VAN DER GRINTEN",     ** Not supported in GDAL??
     437                 :       "WERNER"     ** Not supported in GDAL?? 
     438                 : **/ 
     439              12 :     CPLDebug( "PDS","using projection %s\n\n", map_proj_name.c_str());
     440                 : 
     441              12 :     if ((EQUAL( map_proj_name, "EQUIRECTANGULAR" )) ||
     442                 :         (EQUAL( map_proj_name, "SIMPLE_CYLINDRICAL" )) ||
     443                 :         (EQUAL( map_proj_name, "EQUIDISTANT" )) )  {
     444               7 :         oSRS.SetEquirectangular2 ( 0.0, center_lon, center_lat, 0, 0 );
     445               5 :     } else if (EQUAL( map_proj_name, "ORTHOGRAPHIC" )) {
     446               0 :         oSRS.SetOrthographic ( center_lat, center_lon, 0, 0 );
     447               5 :     } else if (EQUAL( map_proj_name, "SINUSOIDAL" )) {
     448               3 :         oSRS.SetSinusoidal ( center_lon, 0, 0 );
     449               2 :     } else if (EQUAL( map_proj_name, "MERCATOR" )) {
     450               0 :         oSRS.SetMercator ( center_lat, center_lon, 1, 0, 0 );
     451               2 :     } else if (EQUAL( map_proj_name, "STEREOGRAPHIC" )) {
     452               0 :         oSRS.SetStereographic ( center_lat, center_lon, 1, 0, 0 );
     453               2 :     } else if (EQUAL( map_proj_name, "POLAR_STEREOGRAPHIC")) {
     454               0 :         oSRS.SetPS ( center_lat, center_lon, 1, 0, 0 );
     455               2 :     } else if (EQUAL( map_proj_name, "TRANSVERSE_MERCATOR" )) {
     456               0 :         oSRS.SetTM ( center_lat, center_lon, 1, 0, 0 );
     457               2 :     } else if (EQUAL( map_proj_name, "LAMBERT_CONFORMAL_CONIC" )) {
     458                 :         oSRS.SetLCC ( first_std_parallel, second_std_parallel, 
     459               0 :                       center_lat, center_lon, 0, 0 );
     460               2 :     } else if (EQUAL( map_proj_name, "LAMBERT_AZIMUTHAL_EQUAL_AREA" )) {
     461               0 :         oSRS.SetLAEA( center_lat, center_lon, 0, 0 );
     462               2 :     } else if (EQUAL( map_proj_name, "CYLINDRICAL_EQUAL_AREA" )) {
     463               0 :         oSRS.SetCEA  ( first_std_parallel, center_lon, 0, 0 );
     464               2 :     } else if (EQUAL( map_proj_name, "MOLLWEIDE" )) {
     465               0 :         oSRS.SetMollweide ( center_lon, 0, 0 );
     466               2 :     } else if (EQUAL( map_proj_name, "ALBERS" )) {
     467                 :         oSRS.SetACEA ( first_std_parallel, second_std_parallel, 
     468               0 :                        center_lat, center_lon, 0, 0 );
     469               2 :     } else if (EQUAL( map_proj_name, "BONNE" )) {
     470               0 :         oSRS.SetBonne ( first_std_parallel, center_lon, 0, 0 );
     471               2 :     } else if (EQUAL( map_proj_name, "GNOMONIC" )) {
     472               0 :         oSRS.SetGnomonic ( center_lat, center_lon, 0, 0 );
     473               2 :     } else if (EQUAL( map_proj_name, "OBLIQUE_CYLINDRICAL" )) { 
     474                 :         // hope Swiss Oblique Cylindrical is the same
     475               0 :         oSRS.SetSOC ( center_lat, center_lon, 0, 0 );
     476                 :     } else {
     477                 :         CPLDebug( "PDS",
     478                 :                   "Dataset projection %s is not supported. Continuing...",
     479               2 :                   map_proj_name.c_str() );
     480               2 :         bProjectionSet = FALSE;
     481                 :     }
     482                 : 
     483              12 :     if (bProjectionSet) {
     484                 :         //Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
     485              10 :         CPLString proj_target_name = map_proj_name + " " + target_name;
     486              10 :         oSRS.SetProjCS(proj_target_name); //set ProjCS keyword
     487                 :      
     488                 :         //The geographic/geocentric name will be the same basic name as the body name
     489                 :         //'GCS' = Geographic/Geocentric Coordinate System
     490              10 :         CPLString geog_name = "GCS_" + target_name;
     491                 :         
     492                 :         //The datum and sphere names will be the same basic name aas the planet
     493              10 :         CPLString datum_name = "D_" + target_name;
     494              10 :         CPLString sphere_name = target_name; // + "_IAU_IAG");  //Might not be IAU defined so don't add
     495                 :           
     496                 :         //calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
     497              10 :         if ((semi_major - semi_minor) < 0.0000001) 
     498               8 :             iflattening = 0;
     499                 :         else
     500               2 :             iflattening = semi_major / (semi_major - semi_minor);
     501                 :      
     502                 :         //Set the body size but take into consideration which proj is being used to help w/ compatibility
     503                 :         //Notice that most PDS projections are spherical based on the fact that ISIS/PICS are spherical 
     504                 :         //Set the body size but take into consideration which proj is being used to help w/ proj4 compatibility
     505                 :         //The use of a Sphere, polar radius or ellipse here is based on how ISIS does it internally
     506              10 :         if ( ( (EQUAL( map_proj_name, "STEREOGRAPHIC" ) && (fabs(center_lat) == 90)) ) || 
     507                 :              (EQUAL( map_proj_name, "POLAR_STEREOGRAPHIC" )))  
     508                 :         {
     509               0 :             if (bIsGeographic) { 
     510                 :                 //Geograpraphic, so set an ellipse
     511                 :                 oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
     512                 :                                 semi_major, iflattening, 
     513               0 :                                 "Reference_Meridian", 0.0 );
     514                 :             } else {
     515                 :                 //Geocentric, so force a sphere using the semi-minor axis. I hope... 
     516               0 :                 sphere_name += "_polarRadius";
     517                 :                 oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
     518                 :                                 semi_minor, 0.0, 
     519               0 :                                 "Reference_Meridian", 0.0 );
     520                 :             }
     521                 :         }
     522              10 :         else if ( (EQUAL( map_proj_name, "SIMPLE_CYLINDRICAL" )) || 
     523                 :                   (EQUAL( map_proj_name, "EQUIDISTANT" )) || 
     524                 :                   (EQUAL( map_proj_name, "ORTHOGRAPHIC" )) || 
     525                 :                   (EQUAL( map_proj_name, "STEREOGRAPHIC" )) || 
     526                 :                   (EQUAL( map_proj_name, "SINUSOIDAL" )) ) {
     527                 :             //isis uses the spherical equation for these projections so force a sphere
     528                 :             oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
     529                 :                             semi_major, 0.0, 
     530               6 :                             "Reference_Meridian", 0.0 );
     531                 :         } 
     532               4 :         else if (EQUAL( map_proj_name, "EQUIRECTANGULAR" )) { 
     533                 :             //isis uses local radius as a sphere, which is pre-calculated in the PDS label as the semi-major
     534               4 :             sphere_name += "_localRadius";
     535                 :             oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
     536                 :                             semi_major, 0.0, 
     537               4 :                             "Reference_Meridian", 0.0 );
     538                 :         } 
     539                 :         else { 
     540                 :             //All other projections: Mercator, Transverse Mercator, Lambert Conformal, etc.
     541                 :             //Geographic, so set an ellipse
     542               0 :             if (bIsGeographic) {
     543                 :                 oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
     544                 :                                 semi_major, iflattening, 
     545               0 :                                 "Reference_Meridian", 0.0 );
     546                 :             } else { 
     547                 :                 //Geocentric, so force a sphere. I hope... 
     548                 :                 oSRS.SetGeogCS( geog_name, datum_name, sphere_name,
     549                 :                                 semi_major, 0.0, 
     550               0 :                                 "Reference_Meridian", 0.0 );
     551                 :             }
     552                 :         }
     553                 : 
     554                 :         // translate back into a projection string.
     555              10 :         char *pszResult = NULL;
     556              10 :         oSRS.exportToWkt( &pszResult );
     557              10 :         osProjection = pszResult;
     558              10 :         CPLFree( pszResult );
     559                 :     }
     560                 : 
     561                 : /* ==================================================================== */
     562                 : /*      Check for a .prj and world file to override the georeferencing. */
     563                 : /* ==================================================================== */
     564                 :     {
     565              12 :         CPLString osPath, osName;
     566                 :         VSILFILE *fp;
     567                 : 
     568              12 :         osPath = CPLGetPath( pszFilename );
     569              12 :         osName = CPLGetBasename(pszFilename);
     570              12 :         const char  *pszPrjFile = CPLFormCIFilename( osPath, osName, "prj" );
     571                 : 
     572              12 :         fp = VSIFOpenL( pszPrjFile, "r" );
     573              12 :         if( fp != NULL )
     574                 :         {
     575                 :             char  **papszLines;
     576               0 :             OGRSpatialReference oSRS;
     577                 : 
     578               0 :             VSIFCloseL( fp );
     579                 :         
     580               0 :             papszLines = CSLLoad( pszPrjFile );
     581                 : 
     582               0 :             if( oSRS.importFromESRI( papszLines ) == OGRERR_NONE )
     583                 :             {
     584               0 :                 char *pszResult = NULL;
     585               0 :                 oSRS.exportToWkt( &pszResult );
     586               0 :                 osProjection = pszResult;
     587               0 :                 CPLFree( pszResult );
     588                 :             }
     589                 : 
     590               0 :             CSLDestroy( papszLines );
     591              12 :         }
     592                 :     }
     593                 :     
     594              12 :     if( dfULYMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0 )
     595                 :     {
     596              10 :         bGotTransform = TRUE;
     597              10 :         adfGeoTransform[0] = dfULXMap;
     598              10 :         adfGeoTransform[1] = dfXDim;
     599              10 :         adfGeoTransform[2] = 0.0;
     600              10 :         adfGeoTransform[3] = dfULYMap;
     601              10 :         adfGeoTransform[4] = 0.0;
     602              10 :         adfGeoTransform[5] = dfYDim;
     603                 :     }
     604                 :     
     605              12 :     if( !bGotTransform )
     606                 :         bGotTransform = 
     607                 :             GDALReadWorldFile( pszFilename, "psw", 
     608               2 :                                adfGeoTransform );
     609                 : 
     610              12 :     if( !bGotTransform )
     611                 :         bGotTransform = 
     612                 :             GDALReadWorldFile( pszFilename, "wld", 
     613               2 :                                adfGeoTransform );
     614                 : 
     615              12 : }
     616                 : 
     617                 : /************************************************************************/
     618                 : /*                             ParseImage()                             */
     619                 : /************************************************************************/
     620                 : 
     621              12 : int PDSDataset::ParseImage( CPLString osPrefix )
     622                 : {
     623                 : /* ------------------------------------------------------------------- */
     624                 : /*  We assume the user is pointing to the label (ie. .lbl) file.       */
     625                 : /* ------------------------------------------------------------------- */
     626                 :     // IMAGE can be inline or detached and point to an image name
     627                 :     // ^IMAGE = 3
     628                 :     // ^IMAGE             = "GLOBAL_ALBEDO_8PPD.IMG"
     629                 :     // ^IMAGE             = "MEGT90N000CB.IMG"
     630                 :     // ^IMAGE     = ("BLAH.IMG",1)   -- start at record 1 (1 based)
     631                 :     // ^IMAGE     = ("BLAH.IMG")   -- still start at record 1 (equiv of "BLAH.IMG")
     632                 :     // ^IMAGE     = ("BLAH.IMG", 5 <BYTES>) -- start at byte 5 (the fifth byte in the file)
     633                 :     // ^IMAGE             = 10851 <BYTES>
     634                 :     // ^SPECTRAL_QUBE = 5  for multi-band images
     635                 : 
     636              12 :     CPLString osImageKeyword = osPrefix + "^IMAGE";
     637              12 :     CPLString osQube = GetKeyword( osImageKeyword, "" );
     638              12 :     CPLString osTargetFile = GetDescription();
     639                 : 
     640              12 :     if (EQUAL(osQube,"")) {
     641               2 :         osImageKeyword = "^SPECTRAL_QUBE";
     642               2 :         osQube = GetKeyword( osImageKeyword );
     643                 :     }
     644                 : 
     645              12 :     int nQube = atoi(osQube);
     646              12 :     int nDetachedOffset = 0;
     647              12 :     int bDetachedOffsetInBytes = FALSE;
     648                 : 
     649              12 :     if( osQube[0] == '(' )
     650                 :     {
     651               2 :         osQube = "\"";
     652               2 :         osQube += GetKeywordSub( osImageKeyword, 1 );
     653               2 :         osQube +=  "\"";
     654               2 :         nDetachedOffset = atoi(GetKeywordSub( osImageKeyword, 2, "1")) - 1;
     655                 : 
     656                 :         // If this is not explicitly in bytes, then it is assumed to be in
     657                 :         // records, and we need to translate to bytes.
     658               2 :         if (strstr(GetKeywordSub(osImageKeyword,2),"<BYTES>") != NULL)
     659               1 :             bDetachedOffsetInBytes = TRUE;
     660                 :     }
     661                 : 
     662              12 :     if( osQube[0] == '"' )
     663                 :     {
     664               3 :         CPLString osTPath = CPLGetPath(GetDescription());
     665               3 :         CPLString osFilename = osQube;
     666               3 :         CleanString( osFilename );
     667               3 :         osTargetFile = CPLFormCIFilename( osTPath, osFilename, NULL );
     668                 :     }
     669                 : 
     670              12 :     GDALDataType eDataType = GDT_Byte;
     671                 : 
     672                 :     
     673                 :     //image parameters
     674              12 :     int nRows, nCols, nBands = 1;
     675              12 :     int nSkipBytes = 0;
     676                 :     int itype;
     677                 :     int record_bytes;
     678              12 :     char chByteOrder = 'M';  //default to MSB
     679              12 :     double dfNoData = 0.0;
     680                 :  
     681                 :     /* -------------------------------------------------------------------- */
     682                 :     /*      Checks to see if this is raw PDS image not compressed image     */
     683                 :     /*      so ENCODING_TYPE either does not exist or it equals "N/A".      */
     684                 :     /*      Compressed types will not be supported in this routine          */
     685                 :     /* -------------------------------------------------------------------- */
     686                 :     const char *value;
     687                 : 
     688              12 :     CPLString osEncodingType = GetKeyword(osPrefix+"IMAGE.ENCODING_TYPE","N/A");
     689              12 :     CleanString(osEncodingType);
     690              12 :     if ( !EQUAL(osEncodingType.c_str(),"N/A") )
     691                 :     {
     692                 :         CPLError( CE_Failure, CPLE_OpenFailed, 
     693                 :                   "*** PDS image file has an ENCODING_TYPE parameter:\n"
     694                 :                   "*** gdal pds driver does not support compressed image types\n"
     695               0 :                   "found: (%s)\n\n", osEncodingType.c_str() );
     696               0 :         return FALSE;
     697                 :     } 
     698                 :     /**************** end ENCODING_TYPE check ***********************/
     699                 :     
     700                 :     
     701                 :     /***********   Grab layout type (BSQ, BIP, BIL) ************/
     702                 :     //  AXIS_NAME = (SAMPLE,LINE,BAND)
     703                 :     /***********   Grab samples lines band        **************/
     704                 :     /** if AXIS_NAME = "" then Bands=1 and Sample and Lines   **/
     705                 :     /** are there own keywords  "LINES" and "LINE_SAMPLES"    **/
     706                 :     /** if not NULL then CORE_ITEMS keyword i.e. (234,322,2)  **/
     707                 :     /***********************************************************/
     708              12 :     char szLayout[10] = "BSQ"; //default to band seq.
     709              12 :     value = GetKeyword( osPrefix+"IMAGE.AXIS_NAME", "" );
     710              12 :     if (EQUAL(value,"(SAMPLE,LINE,BAND)") ) {
     711               0 :         strcpy(szLayout,"BSQ");
     712               0 :         nCols = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",1));
     713               0 :         nRows = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",2));
     714               0 :         nBands = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",3));
     715                 :     }
     716              12 :     else if (EQUAL(value,"(BAND,LINE,SAMPLE)") ) {
     717               0 :         strcpy(szLayout,"BIP");
     718               0 :         nBands = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",1));
     719               0 :         nRows = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",2));
     720               0 :         nCols = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",3));
     721                 :     }
     722              12 :     else if (EQUAL(value,"(SAMPLE,BAND,LINE)") ) {
     723               0 :         strcpy(szLayout,"BIL");
     724               0 :         nCols = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",1));
     725               0 :         nBands = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",2));
     726               0 :         nRows = atoi(GetKeywordSub(osPrefix+"IMAGE.CORE_ITEMS",3));
     727                 :     }
     728              12 :     else if ( EQUAL(value,"") ) {
     729              12 :         strcpy(szLayout,"BSQ");
     730              12 :         nCols = atoi(GetKeyword(osPrefix+"IMAGE.LINE_SAMPLES",""));
     731              12 :         nRows = atoi(GetKeyword(osPrefix+"IMAGE.LINES",""));
     732              12 :         nBands = atoi(GetKeyword(osPrefix+"IMAGE.BANDS","1"));
     733                 :     }
     734                 :     else {
     735                 :         CPLError( CE_Failure, CPLE_OpenFailed, 
     736               0 :                   "%s layout not supported. Abort\n\n", value);
     737               0 :         return FALSE;
     738                 :     }
     739                 :     
     740                 :     /***********   Grab Qube record bytes  **********/
     741              12 :     record_bytes = atoi(GetKeyword(osPrefix+"IMAGE.RECORD_BYTES"));
     742              12 :     if (record_bytes == 0)
     743              12 :         record_bytes = atoi(GetKeyword(osPrefix+"RECORD_BYTES"));
     744                 : 
     745                 :     // this can happen with "record_type = undefined". 
     746              12 :     if( record_bytes == 0 )
     747               0 :         record_bytes = 1;
     748                 : 
     749              12 :     if( nQube >0 && osQube.find("<BYTES>") != CPLString::npos )
     750               0 :         nSkipBytes = nQube - 1;
     751              12 :     else if (nQube > 0 )
     752               7 :         nSkipBytes = (nQube - 1) * record_bytes;
     753               5 :     else if( nDetachedOffset > 0 )
     754                 :     {
     755               1 :         if (bDetachedOffsetInBytes)
     756               1 :             nSkipBytes = nDetachedOffset;
     757                 :         else
     758               0 :             nSkipBytes = nDetachedOffset * record_bytes;
     759                 :     }
     760                 :     else
     761               4 :         nSkipBytes = 0;     
     762                 : 
     763              12 :     nSkipBytes += atoi(GetKeyword(osPrefix+"IMAGE.LINE_PREFIX_BYTES",""));
     764                 :     
     765                 :     /***********   Grab SAMPLE_TYPE *****************/
     766                 :     /** if keyword not found leave as "M" or "MSB" **/
     767              12 :     CPLString osST = GetKeyword( osPrefix+"IMAGE.SAMPLE_TYPE" );
     768              12 :     if( osST.size() >= 2 && osST[0] == '"' && osST[osST.size()-1] == '"' )
     769               0 :         osST = osST.substr( 1, osST.size() - 2 );
     770                 : 
     771              12 :     if( (EQUAL(osST,"LSB_INTEGER")) || 
     772                 :         (EQUAL(osST,"LSB")) || // just incase
     773                 :         (EQUAL(osST,"LSB_UNSIGNED_INTEGER")) || 
     774                 :         (EQUAL(osST,"LSB_SIGNED_INTEGER")) || 
     775                 :         (EQUAL(osST,"UNSIGNED_INTEGER")) || 
     776                 :         (EQUAL(osST,"VAX_REAL")) || 
     777                 :         (EQUAL(osST,"VAX_INTEGER")) || 
     778                 :         (EQUAL(osST,"PC_INTEGER")) ||  //just incase 
     779                 :         (EQUAL(osST,"PC_REAL")) ) {
     780               8 :         chByteOrder = 'I';
     781                 :     }
     782                 : 
     783                 :     /**** Grab format type - pds supports 1,2,4,8,16,32,64 (in theory) **/
     784                 :     /**** I have only seen 8, 16, 32 (float) in released datasets      **/
     785              12 :     itype = atoi(GetKeyword(osPrefix+"IMAGE.SAMPLE_BITS",""));
     786              12 :     switch(itype) {
     787                 :       case 8 :
     788               7 :         eDataType = GDT_Byte;
     789               7 :         dfNoData = NULL1;
     790               7 :         break;
     791                 :       case 16 :
     792               3 :         if( strstr(osST,"UNSIGNED") != NULL )
     793               2 :             eDataType = GDT_UInt16;
     794                 :         else
     795               1 :             eDataType = GDT_Int16;
     796               3 :         dfNoData = NULL2;
     797               3 :         break;
     798                 :       case 32 :
     799               0 :         eDataType = GDT_Float32;
     800               0 :         dfNoData = NULL3;
     801               0 :         break;
     802                 :       case 64 :
     803               0 :         eDataType = GDT_Float64;
     804               0 :         dfNoData = NULL3;
     805               0 :         break;
     806                 :       default :
     807                 :         CPLError( CE_Failure, CPLE_AppDefined,
     808                 :                   "Sample_bits of %d is not supported in this gdal PDS reader.",
     809               2 :                   itype); 
     810               2 :         return FALSE;
     811                 :     }
     812                 : 
     813                 : /* -------------------------------------------------------------------- */
     814                 : /*      Is there a specific nodata value in the file? Either the        */
     815                 : /*      MISSING or MISSING_CONSTANT keywords are nodata.                */
     816                 : /* -------------------------------------------------------------------- */
     817              10 :     if( GetKeyword( osPrefix+"IMAGE.MISSING", NULL ) != NULL )
     818               2 :         dfNoData = CPLAtofM( GetKeyword( osPrefix+"IMAGE.MISSING", "" ) );
     819                 : 
     820              10 :     if( GetKeyword( osPrefix+"IMAGE.MISSING_CONSTANT", NULL ) != NULL )
     821               1 :         dfNoData = CPLAtofM( GetKeyword( osPrefix+"IMAGE.MISSING_CONSTANT",""));
     822                 : 
     823                 : /* -------------------------------------------------------------------- */
     824                 : /*      Did we get the required keywords?  If not we return with        */
     825                 : /*      this never having been considered to be a match. This isn't     */
     826                 : /*      an error!                                                       */
     827                 : /* -------------------------------------------------------------------- */
     828              10 :     if( nRows < 1 || nCols < 1 || nBands < 1 )
     829                 :     {
     830                 :         CPLError( CE_Failure, CPLE_AppDefined,
     831                 :                   "File %s appears to be a PDS file, but failed to find some required keywords.", 
     832               0 :                   GetDescription() );
     833               0 :         return FALSE;
     834                 :     }
     835                 : 
     836                 : /* -------------------------------------------------------------------- */
     837                 : /*      Capture some information from the file that is of interest.     */
     838                 : /* -------------------------------------------------------------------- */
     839              10 :     nRasterXSize = nCols;
     840              10 :     nRasterYSize = nRows;
     841                 : 
     842                 : /* -------------------------------------------------------------------- */
     843                 : /*      Open target binary file.                                        */
     844                 : /* -------------------------------------------------------------------- */
     845                 :     
     846              10 :     if( eAccess == GA_ReadOnly )
     847                 :     {
     848              10 :         fpImage = VSIFOpenL( osTargetFile, "rb" );
     849              10 :         if( fpImage == NULL )
     850                 :         {
     851                 :             CPLError( CE_Failure, CPLE_OpenFailed, 
     852                 :                     "Failed to open %s.\n%s", 
     853                 :                     osTargetFile.c_str(),
     854               0 :                     VSIStrerror( errno ) );
     855               0 :             return FALSE;
     856                 :         }
     857                 :     }
     858                 :     else
     859                 :     {
     860               0 :         fpImage = VSIFOpenL( osTargetFile, "r+b" );
     861               0 :         if( fpImage == NULL )
     862                 :         {
     863                 :             CPLError( CE_Failure, CPLE_OpenFailed, 
     864                 :                     "Failed to open %s with write permission.\n%s", 
     865                 :                     osTargetFile.c_str(),
     866               0 :                     VSIStrerror( errno ) );
     867               0 :             return FALSE;
     868                 :         }
     869                 :     }
     870                 : 
     871                 : /* -------------------------------------------------------------------- */
     872                 : /*      Compute the line offset.                                        */
     873                 : /* -------------------------------------------------------------------- */
     874              10 :     int     nItemSize = GDALGetDataTypeSize(eDataType)/8;
     875              10 :     int     nLineOffset = record_bytes;
     876                 :     int     nPixelOffset, nBandOffset;
     877                 : 
     878              10 :     if( EQUAL(szLayout,"BIP") )
     879                 :     {
     880               0 :         nPixelOffset = nItemSize * nBands;
     881               0 :         nBandOffset = nItemSize;
     882                 :         nLineOffset = ((nPixelOffset * nCols + record_bytes - 1)/record_bytes)
     883               0 :             * record_bytes;
     884                 :     }
     885              10 :     else if( EQUAL(szLayout,"BSQ") )
     886                 :     {
     887              10 :         nPixelOffset = nItemSize;
     888                 :         nLineOffset = ((nPixelOffset * nCols + record_bytes - 1)/record_bytes)
     889              10 :             * record_bytes;
     890              10 :         nBandOffset = nLineOffset * nRows;
     891                 :     }
     892                 :     else /* assume BIL */
     893                 :     {
     894               0 :         nPixelOffset = nItemSize;
     895               0 :         nBandOffset = nItemSize * nCols;
     896                 :         nLineOffset = ((nBandOffset * nCols + record_bytes - 1)/record_bytes)
     897               0 :             * record_bytes;
     898                 :     }
     899                 :     
     900                 : /* -------------------------------------------------------------------- */
     901                 : /*      Create band information objects.                                */
     902                 : /* -------------------------------------------------------------------- */
     903                 :     int i;
     904                 : 
     905              20 :     for( i = 0; i < nBands; i++ )
     906                 :     {
     907                 :         RawRasterBand *poBand;
     908                 : 
     909                 :         poBand = 
     910                 :             new RawRasterBand( this, i+1, fpImage,
     911                 :                                nSkipBytes + nBandOffset * i, 
     912                 :                                nPixelOffset, nLineOffset, eDataType,
     913                 : #ifdef CPL_LSB                               
     914                 :                                chByteOrder == 'I' || chByteOrder == 'L',
     915                 : #else
     916                 :                                chByteOrder == 'M',
     917                 : #endif        
     918              10 :                                TRUE );
     919                 : 
     920              10 :         if( nBands == 1 )
     921                 :         {
     922              10 :             const char* pszMin = GetKeyword(osPrefix+"IMAGE.MINIMUM", NULL);
     923              10 :             const char* pszMax = GetKeyword(osPrefix+"IMAGE.MAXIMUM", NULL);
     924              10 :             const char* pszMean = GetKeyword(osPrefix+"IMAGE.MEAN", NULL);
     925              10 :             const char* pszStdDev= GetKeyword(osPrefix+"IMAGE.STANDARD_DEVIATION", NULL);
     926              10 :             if (pszMin != NULL && pszMax != NULL &&
     927                 :                 pszMean != NULL && pszStdDev != NULL)
     928                 :             {
     929                 :                 poBand->SetStatistics( CPLAtofM(pszMin),
     930                 :                                        CPLAtofM(pszMax),
     931                 :                                        CPLAtofM(pszMean),
     932               0 :                                        CPLAtofM(pszStdDev));
     933                 :             }
     934                 :         }
     935                 :         
     936              10 :         poBand->SetNoDataValue( dfNoData );
     937                 : 
     938              10 :         SetBand( i+1, poBand );
     939                 : 
     940                 :         // Set offset/scale values at the PAM level.
     941                 :         poBand->SetOffset( 
     942              10 :             CPLAtofM(GetKeyword(osPrefix+"IMAGE.OFFSET","0.0")));
     943                 :         poBand->SetScale( 
     944              10 :             CPLAtofM(GetKeyword(osPrefix+"IMAGE.SCALING_FACTOR","1.0")));
     945                 :     }
     946                 : 
     947              10 :     return TRUE;
     948                 : }
     949                 : 
     950                 : /************************************************************************/
     951                 : /* ==================================================================== */
     952                 : /*                         PDSWrapperRasterBand                         */
     953                 : /*                                                                      */
     954                 : /*      proxy for the jp2 or other compressed bands.                    */
     955                 : /* ==================================================================== */
     956                 : /************************************************************************/
     957                 : class PDSWrapperRasterBand : public GDALProxyRasterBand
     958                 : {
     959                 :   GDALRasterBand* poBaseBand;
     960                 : 
     961                 :   protected:
     962              22 :     virtual GDALRasterBand* RefUnderlyingRasterBand() { return poBaseBand; }
     963                 : 
     964                 :   public:
     965               2 :     PDSWrapperRasterBand( GDALRasterBand* poBaseBand ) 
     966               2 :         {
     967               2 :             this->poBaseBand = poBaseBand;
     968               2 :             eDataType = poBaseBand->GetRasterDataType();
     969               2 :             poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
     970               2 :         }
     971               2 :     ~PDSWrapperRasterBand() {}
     972                 : };
     973                 : 
     974                 : /************************************************************************/
     975                 : /*                       ParseCompressedImage()                         */
     976                 : /************************************************************************/
     977                 : 
     978               2 : int PDSDataset::ParseCompressedImage()
     979                 : 
     980                 : {
     981               2 :     CPLString osFileName = GetKeyword( "COMPRESSED_FILE.FILE_NAME", "" );
     982               2 :     CleanString( osFileName );
     983                 : 
     984               2 :     CPLString osPath = CPLGetPath(GetDescription());
     985               2 :     CPLString osFullFileName = CPLFormFilename( osPath, osFileName, NULL );
     986                 :     int iBand;
     987                 : 
     988               2 :     poCompressedDS = (GDALDataset*) GDALOpen( osFullFileName, GA_ReadOnly );
     989                 :     
     990               2 :     if( poCompressedDS == NULL )
     991               0 :         return FALSE;
     992                 : 
     993               2 :     nRasterXSize = poCompressedDS->GetRasterXSize();
     994               2 :     nRasterYSize = poCompressedDS->GetRasterYSize();
     995                 : 
     996               4 :     for( iBand = 0; iBand < poCompressedDS->GetRasterCount(); iBand++ )
     997                 :     {
     998               2 :         SetBand( iBand+1, new PDSWrapperRasterBand( poCompressedDS->GetRasterBand( iBand+1 ) ) );
     999                 :     }
    1000                 :     
    1001               2 :     return TRUE;
    1002                 : }
    1003                 : 
    1004                 : /************************************************************************/
    1005                 : /*                              Identify()                              */
    1006                 : /************************************************************************/
    1007                 : 
    1008           11498 : int PDSDataset::Identify( GDALOpenInfo * poOpenInfo )
    1009                 : 
    1010                 : {
    1011           11498 :     if( poOpenInfo->pabyHeader == NULL )
    1012           10392 :         return FALSE;
    1013                 : 
    1014            1106 :     return strstr((char*)poOpenInfo->pabyHeader,"PDS_VERSION_ID") != NULL;
    1015                 : }
    1016                 : 
    1017                 : /************************************************************************/
    1018                 : /*                                Open()                                */
    1019                 : /************************************************************************/
    1020                 : 
    1021            2253 : GDALDataset *PDSDataset::Open( GDALOpenInfo * poOpenInfo )
    1022                 : {
    1023            2253 :     if( !Identify( poOpenInfo ) )
    1024            2239 :         return NULL;
    1025                 : 
    1026              14 :     if( strstr((const char *)poOpenInfo->pabyHeader,"PDS3") == NULL )
    1027                 :     {
    1028                 :         CPLError( CE_Failure, CPLE_OpenFailed,
    1029               0 :                   "It appears this is an older PDS image type.  Only PDS_VERSION_ID = PDS3 are currently supported by this gdal PDS reader.");
    1030               0 :         return NULL;
    1031                 :     }
    1032                 : 
    1033                 : /* -------------------------------------------------------------------- */
    1034                 : /*      Open and parse the keyword header.  Sometimes there is stuff    */
    1035                 : /*      before the PDS_VERSION_ID, which we want to ignore.             */
    1036                 : /* -------------------------------------------------------------------- */
    1037              14 :     VSILFILE *fpQube = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
    1038                 : 
    1039              14 :     if( fpQube == NULL )
    1040               0 :         return NULL;
    1041                 : 
    1042                 :     PDSDataset  *poDS;
    1043                 : 
    1044              14 :     poDS = new PDSDataset();
    1045              14 :     poDS->SetDescription( poOpenInfo->pszFilename );
    1046              14 :     poDS->eAccess = poOpenInfo->eAccess;
    1047                 :     
    1048              14 :     const char* pszPDSVersionID = strstr((const char *)poOpenInfo->pabyHeader,"PDS_VERSION_ID");
    1049              14 :     int nOffset = 0;
    1050              14 :     if (pszPDSVersionID)
    1051              14 :         nOffset = pszPDSVersionID - (const char *)poOpenInfo->pabyHeader;
    1052                 : 
    1053              14 :     if( ! poDS->oKeywords.Ingest( fpQube, nOffset ) )
    1054                 :     {
    1055               0 :         delete poDS;
    1056               0 :         VSIFCloseL( fpQube );
    1057               0 :         return NULL;
    1058                 :     }
    1059              14 :     VSIFCloseL( fpQube );
    1060                 : 
    1061                 : /* -------------------------------------------------------------------- */
    1062                 : /*      Is this a compressed image with COMPRESSED_FILE subdomain?      */
    1063                 : /*                                                                      */
    1064                 : /*      The corresponding parse operations will read keywords,          */
    1065                 : /*      establish bands and raster size.                                */
    1066                 : /* -------------------------------------------------------------------- */
    1067              14 :     CPLString osEncodingType = poDS->GetKeyword( "COMPRESSED_FILE.ENCODING_TYPE", "" );
    1068                 : 
    1069              28 :     if( osEncodingType.size() != 0 )
    1070                 :     {
    1071               2 :         if( !poDS->ParseCompressedImage() )
    1072                 :         {
    1073               0 :             delete poDS;
    1074               0 :             return NULL;
    1075                 :         }
    1076                 :     }
    1077                 :     else
    1078                 :     {
    1079              12 :         CPLString osPrefix;
    1080              12 :       CPLString osObject = poDS->GetKeyword( "UNCOMPRESSED_FILE.IMAGE.NAME", "");
    1081                 : 
    1082              24 :         if( osObject != "" )
    1083               1 :             osPrefix = "UNCOMPRESSED_FILE.";
    1084                 :         
    1085              12 :         if( !poDS->ParseImage(osPrefix) )
    1086                 :         {
    1087               2 :             delete poDS;
    1088               2 :             return NULL;
    1089               0 :         }
    1090                 :     }
    1091                 : 
    1092                 : /* -------------------------------------------------------------------- */
    1093                 : /*      Set the coordinate system and geotransform.                     */
    1094                 : /* -------------------------------------------------------------------- */
    1095              12 :     poDS->ParseSRS();
    1096                 : 
    1097                 : /* -------------------------------------------------------------------- */
    1098                 : /*      Transfer a few interesting keywords as metadata.                */
    1099                 : /* -------------------------------------------------------------------- */
    1100                 :     int i;
    1101                 :     static const char *apszKeywords[] = 
    1102                 :         { "FILTER_NAME", "DATA_SET_ID", "PRODUCT_ID", 
    1103                 :           "PRODUCER_INSTITUTION_NAME", "PRODUCT_TYPE", "MISSION_NAME",
    1104                 :           "SPACECRAFT_NAME", "INSTRUMENT_NAME", "INSTRUMENT_ID", 
    1105                 :           "TARGET_NAME", "CENTER_FILTER_WAVELENGTH", "BANDWIDTH",
    1106                 :           "PRODUCT_CREATION_TIME", "NOTE",
    1107                 :           NULL };
    1108                 :     
    1109             180 :     for( i = 0; apszKeywords[i] != NULL; i++ )
    1110                 :     {
    1111             168 :         const char *pszKeywordValue = poDS->GetKeyword( apszKeywords[i] );
    1112                 : 
    1113             168 :         if( pszKeywordValue != NULL )
    1114             168 :             poDS->SetMetadataItem( apszKeywords[i], pszKeywordValue );
    1115                 :     }
    1116                 : 
    1117                 : /* -------------------------------------------------------------------- */
    1118                 : /*      Initialize any PAM information.                                 */
    1119                 : /* -------------------------------------------------------------------- */
    1120              12 :     poDS->TryLoadXML();
    1121                 : 
    1122                 : /* -------------------------------------------------------------------- */
    1123                 : /*      Check for overviews.                                            */
    1124                 : /* -------------------------------------------------------------------- */
    1125              12 :     poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
    1126                 : 
    1127              12 :     return( poDS );
    1128                 : }
    1129                 : 
    1130                 : /************************************************************************/
    1131                 : /*                             GetKeyword()                             */
    1132                 : /************************************************************************/
    1133                 : 
    1134             557 : const char *PDSDataset::GetKeyword( std::string osPath,
    1135                 :                                     const char *pszDefault )
    1136                 : 
    1137                 : {
    1138             557 :     return oKeywords.GetKeyword( osPath.c_str(), pszDefault );
    1139                 : }
    1140                 : 
    1141                 : /************************************************************************/
    1142                 : /*                            GetKeywordSub()                           */
    1143                 : /************************************************************************/
    1144                 : 
    1145               6 : const char *PDSDataset::GetKeywordSub( std::string osPath,
    1146                 :                                        int iSubscript,
    1147                 :                                        const char *pszDefault )
    1148                 : 
    1149                 : {
    1150               6 :     const char *pszResult = oKeywords.GetKeyword( osPath.c_str(), NULL );
    1151                 :     
    1152               6 :     if( pszResult == NULL )
    1153               0 :         return pszDefault;
    1154                 : 
    1155               6 :     if( pszResult[0] != '(' ) 
    1156               0 :         return pszDefault;
    1157                 : 
    1158                 :     char **papszTokens = CSLTokenizeString2( pszResult, "(,)", 
    1159               6 :                                              CSLT_HONOURSTRINGS );
    1160                 : 
    1161               6 :     if( iSubscript <= CSLCount(papszTokens) )
    1162                 :     {
    1163               6 :         osTempResult = papszTokens[iSubscript-1];
    1164               6 :         CSLDestroy( papszTokens );
    1165               6 :         return osTempResult.c_str();
    1166                 :     }
    1167                 :     else
    1168                 :     {
    1169               0 :         CSLDestroy( papszTokens );
    1170               0 :         return pszDefault;
    1171                 :     }
    1172                 : }
    1173                 : 
    1174                 : /************************************************************************/
    1175                 : /*                            GetKeywordUnit()                          */
    1176                 : /************************************************************************/
    1177                 : 
    1178              10 : const char *PDSDataset::GetKeywordUnit( const char *pszPath, 
    1179                 :                                          int iSubscript,
    1180                 :                                          const char *pszDefault )
    1181                 : 
    1182                 : {
    1183              10 :     const char *pszResult = oKeywords.GetKeyword( pszPath, NULL );
    1184                 :     
    1185              10 :     if( pszResult == NULL )
    1186               0 :         return pszDefault;
    1187                 :  
    1188                 :     char **papszTokens = CSLTokenizeString2( pszResult, "</>", 
    1189              10 :                                              CSLT_HONOURSTRINGS );
    1190                 : 
    1191              10 :     if( iSubscript <= CSLCount(papszTokens) )
    1192                 :     {
    1193               8 :         osTempResult = papszTokens[iSubscript-1];
    1194               8 :         CSLDestroy( papszTokens );
    1195               8 :         return osTempResult.c_str();
    1196                 :     }
    1197                 :     else
    1198                 :     {
    1199               2 :         CSLDestroy( papszTokens );
    1200               2 :         return pszDefault;
    1201                 :     }
    1202                 : }
    1203                 : 
    1204                 : /************************************************************************/
    1205                 : /*                            CleanString()                             */
    1206                 : /*                                                                      */
    1207                 : /* Removes single or double quotes, and converts spaces to underscores. */
    1208                 : /* The change is made in-place to CPLString.                            */
    1209                 : /************************************************************************/
    1210                 : 
    1211              41 : void PDSDataset::CleanString( CPLString &osInput )
    1212                 : 
    1213                 : {
    1214              41 :    if(  ( osInput.size() < 2 ) ||
    1215                 :         ((osInput.at(0) != '"'   || osInput.at(osInput.size()-1) != '"' ) &&
    1216                 :         ( osInput.at(0) != '\'' || osInput.at(osInput.size()-1) != '\'')) )
    1217              23 :         return;
    1218                 : 
    1219              18 :     char *pszWrk = CPLStrdup(osInput.c_str() + 1);
    1220                 :     int i;
    1221                 : 
    1222              18 :     pszWrk[strlen(pszWrk)-1] = '\0';
    1223                 :     
    1224             180 :     for( i = 0; pszWrk[i] != '\0'; i++ )
    1225                 :     {
    1226             162 :         if( pszWrk[i] == ' ' )
    1227               3 :             pszWrk[i] = '_';
    1228                 :     }
    1229                 : 
    1230              18 :     osInput = pszWrk;
    1231              18 :     CPLFree( pszWrk );
    1232                 : }
    1233                 : /************************************************************************/
    1234                 : /*                         GDALRegister_PDS()                         */
    1235                 : /************************************************************************/
    1236                 : 
    1237             558 : void GDALRegister_PDS()
    1238                 : 
    1239                 : {
    1240                 :     GDALDriver  *poDriver;
    1241                 : 
    1242             558 :     if( GDALGetDriverByName( "PDS" ) == NULL )
    1243                 :     {
    1244             537 :         poDriver = new GDALDriver();
    1245                 :         
    1246             537 :         poDriver->SetDescription( "PDS" );
    1247                 :         poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
    1248             537 :                                    "NASA Planetary Data System" );
    1249                 :         poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
    1250             537 :                                    "frmt_pds.html" );
    1251             537 :         poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
    1252                 : 
    1253             537 :         poDriver->pfnOpen = PDSDataset::Open;
    1254             537 :         poDriver->pfnIdentify = PDSDataset::Identify;
    1255                 : 
    1256             537 :         GetGDALDriverManager()->RegisterDriver( poDriver );
    1257                 :     }
    1258             558 : }
    1259                 : 

Generated by: LCOV version 1.7