LCOV - code coverage report
Current view: directory - frmts/raw - ctable2dataset.cpp (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 132 119 90.2 %
Date: 2012-12-26 Functions: 13 9 69.2 %

       1                 : /******************************************************************************
       2                 :  * $Id$
       3                 :  *
       4                 :  * Project:  Horizontal Datum Formats
       5                 :  * Purpose:  Implementation of the CTable2 format, a PROJ.4 specific format
       6                 :  *           that is more compact (due to a lack of unused error band) than NTv2
       7                 :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       8                 :  *
       9                 :  ******************************************************************************
      10                 :  * Copyright (c) 2012, Frank Warmerdam
      11                 :  *
      12                 :  * Permission is hereby granted, free of charge, to any person obtaining a
      13                 :  * copy of this software and associated documentation files (the "Software"),
      14                 :  * to deal in the Software without restriction, including without limitation
      15                 :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16                 :  * and/or sell copies of the Software, and to permit persons to whom the
      17                 :  * Software is furnished to do so, subject to the following conditions:
      18                 :  *
      19                 :  * The above copyright notice and this permission notice shall be included
      20                 :  * in all copies or substantial portions of the Software.
      21                 :  *
      22                 :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      23                 :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      24                 :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      25                 :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      26                 :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      27                 :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      28                 :  * DEALINGS IN THE SOFTWARE.
      29                 :  ****************************************************************************/
      30                 : 
      31                 : #include "rawdataset.h"
      32                 : #include "cpl_string.h"
      33                 : #include "ogr_srs_api.h"
      34                 : 
      35                 : CPL_CVSID("$Id$");
      36                 : 
      37                 : #ifndef M_PI
      38                 : #define M_PI    3.14159265358979323846
      39                 : #endif
      40                 : 
      41                 : /************************************************************************/
      42                 : /* ==================================================================== */
      43                 : /*        CTable2Dataset        */
      44                 : /* ==================================================================== */
      45                 : /************************************************************************/
      46                 : 
      47                 : class CTable2Dataset : public RawDataset
      48                 : {
      49                 :   public:
      50                 :     VSILFILE  *fpImage; // image data file.
      51                 : 
      52                 :     double      adfGeoTransform[6];
      53                 : 
      54                 :   public:
      55                 :         CTable2Dataset();
      56                 :               ~CTable2Dataset();
      57                 :     
      58                 :     virtual CPLErr SetGeoTransform( double * padfTransform );
      59                 :     virtual CPLErr GetGeoTransform( double * padfTransform );
      60                 :     virtual const char *GetProjectionRef();
      61                 :     virtual void   FlushCache(void);
      62                 : 
      63                 :     static GDALDataset *Open( GDALOpenInfo * );
      64                 :     static int          Identify( GDALOpenInfo * );
      65                 :     static GDALDataset *Create( const char * pszFilename,
      66                 :                                 int nXSize, int nYSize, int nBands,
      67                 :                                 GDALDataType eType, char ** papszOptions );
      68                 : };
      69                 : 
      70                 : /************************************************************************/
      71                 : /* ==================================================================== */
      72                 : /*        CTable2Dataset        */
      73                 : /* ==================================================================== */
      74                 : /************************************************************************/
      75                 : 
      76                 : /************************************************************************/
      77                 : /*                             CTable2Dataset()                          */
      78                 : /************************************************************************/
      79                 : 
      80               6 : CTable2Dataset::CTable2Dataset()
      81                 : {
      82               6 :     fpImage = NULL;
      83               6 : }
      84                 : 
      85                 : /************************************************************************/
      86                 : /*                            ~CTable2Dataset()                          */
      87                 : /************************************************************************/
      88                 : 
      89               6 : CTable2Dataset::~CTable2Dataset()
      90                 : 
      91                 : {
      92               6 :     FlushCache();
      93                 : 
      94               6 :     if( fpImage != NULL )
      95               6 :         VSIFCloseL( fpImage );
      96               6 : }
      97                 : 
      98                 : /************************************************************************/
      99                 : /*                             FlushCache()                             */
     100                 : /************************************************************************/
     101                 : 
     102               6 : void CTable2Dataset::FlushCache()
     103                 : 
     104                 : {
     105               6 :     RawDataset::FlushCache();
     106               6 : }
     107                 : 
     108                 : /************************************************************************/
     109                 : /*                              Identify()                              */
     110                 : /************************************************************************/
     111                 : 
     112           11929 : int CTable2Dataset::Identify( GDALOpenInfo *poOpenInfo )
     113                 : 
     114                 : {
     115           11929 :     if( poOpenInfo->nHeaderBytes < 64 )
     116           11360 :         return FALSE;
     117                 : 
     118             569 :     if( !EQUALN((const char *)poOpenInfo->pabyHeader + 0, "CTABLE V2", 9 ) )
     119             563 :         return FALSE;
     120                 : 
     121               6 :     return TRUE;
     122                 : }
     123                 : 
     124                 : /************************************************************************/
     125                 : /*                                Open()                                */
     126                 : /************************************************************************/
     127                 : 
     128            1836 : GDALDataset *CTable2Dataset::Open( GDALOpenInfo * poOpenInfo )
     129                 : 
     130                 : {
     131            1836 :     if( !Identify( poOpenInfo ) )
     132            1830 :         return NULL;
     133                 :         
     134                 : /* -------------------------------------------------------------------- */
     135                 : /*      Create a corresponding GDALDataset.                             */
     136                 : /* -------------------------------------------------------------------- */
     137                 :     CTable2Dataset  *poDS;
     138                 : 
     139               6 :     poDS = new CTable2Dataset();
     140               6 :     poDS->eAccess = poOpenInfo->eAccess;
     141                 : 
     142                 : /* -------------------------------------------------------------------- */
     143                 : /*      Open the file.                                                  */
     144                 : /* -------------------------------------------------------------------- */
     145               6 :     CPLString osFilename;
     146              12 :     osFilename = poOpenInfo->pszFilename;
     147                 :     
     148               6 :     if( poOpenInfo->eAccess == GA_ReadOnly )
     149               3 :         poDS->fpImage = VSIFOpenL( osFilename, "rb" );
     150                 :     else
     151               3 :         poDS->fpImage = VSIFOpenL( osFilename, "rb+" );
     152                 : 
     153               6 :     if( poDS->fpImage == NULL )
     154                 :     {
     155               0 :         delete poDS;
     156               0 :         return NULL;
     157                 :     }
     158                 : 
     159                 : /* -------------------------------------------------------------------- */
     160                 : /*      Read the file header.                                           */
     161                 : /* -------------------------------------------------------------------- */
     162                 :     char  achHeader[160];
     163               6 :     CPLString osDescription;
     164                 : 
     165               6 :     VSIFSeekL( poDS->fpImage, 0, SEEK_SET );
     166               6 :     VSIFReadL( achHeader, 1, 160, poDS->fpImage );
     167                 : 
     168               6 :     achHeader[16+79] = '\0';
     169               6 :     osDescription = (const char *) achHeader+16;
     170               6 :     osDescription.Trim();
     171               6 :     poDS->SetMetadataItem( "DESCRIPTION", osDescription );
     172                 : 
     173                 : /* -------------------------------------------------------------------- */
     174                 : /*      Convert from LSB to local machine byte order.                   */
     175                 : /* -------------------------------------------------------------------- */
     176                 :     CPL_LSBPTR64( achHeader + 96 );
     177                 :     CPL_LSBPTR64( achHeader + 104 );
     178                 :     CPL_LSBPTR64( achHeader + 112 );
     179                 :     CPL_LSBPTR64( achHeader + 120 );
     180                 :     CPL_LSBPTR32( achHeader + 128 );
     181                 :     CPL_LSBPTR32( achHeader + 132 );
     182                 : 
     183                 : /* -------------------------------------------------------------------- */
     184                 : /*      Extract size, and geotransform.                                 */
     185                 : /* -------------------------------------------------------------------- */
     186                 :     GUInt32 nRasterXSize, nRasterYSize;
     187                 :     
     188               6 :     memcpy( &nRasterXSize, achHeader + 128, 4 );
     189               6 :     memcpy( &nRasterYSize, achHeader + 132, 4 );
     190               6 :     poDS->nRasterXSize = nRasterXSize;
     191               6 :     poDS->nRasterYSize = nRasterYSize;
     192                 : 
     193                 :     int i;
     194                 :     double adfValues[4];
     195               6 :     memcpy( adfValues, achHeader + 96, sizeof(double)*4 );
     196                 : 
     197              30 :     for( i = 0; i < 4; i++ )
     198              24 :         adfValues[i] *= 180/M_PI; // Radians to degrees.
     199                 :     
     200                 : 
     201               6 :     poDS->adfGeoTransform[0] = adfValues[0] - adfValues[2]*0.5;
     202               6 :     poDS->adfGeoTransform[1] = adfValues[2];
     203               6 :     poDS->adfGeoTransform[2] = 0.0;
     204               6 :     poDS->adfGeoTransform[3] = adfValues[1] + adfValues[3]*(nRasterYSize-0.5);
     205               6 :     poDS->adfGeoTransform[4] = 0.0;
     206               6 :     poDS->adfGeoTransform[5] = -adfValues[3];
     207                 : 
     208                 : /* -------------------------------------------------------------------- */
     209                 : /*      Setup the bands.                                                */
     210                 : /* -------------------------------------------------------------------- */
     211                 :     RawRasterBand *poBand = 
     212                 :         new RawRasterBand( poDS, 1, poDS->fpImage, 
     213                 :                            160 + 4 + nRasterXSize * (nRasterYSize-1) * 2 * 4,
     214                 :                            8, -8 * nRasterXSize, 
     215               6 :                            GDT_Float32, CPL_IS_LSB, TRUE, FALSE );
     216               6 :     poBand->SetDescription( "Latitude Offset (radians)" );
     217               6 :     poDS->SetBand( 1, poBand );
     218                 :     
     219                 :     poBand = 
     220                 :         new RawRasterBand( poDS, 2, poDS->fpImage, 
     221                 :                            160 + nRasterXSize * (nRasterYSize-1) * 2 * 4,
     222                 :                            8, -8 * nRasterXSize, 
     223               6 :                            GDT_Float32, CPL_IS_LSB, TRUE, FALSE );
     224               6 :     poBand->SetDescription( "Longitude Offset (radians)" );
     225               6 :     poDS->SetBand( 2, poBand );
     226                 : 
     227                 : /* -------------------------------------------------------------------- */
     228                 : /*      Initialize any PAM information.                                 */
     229                 : /* -------------------------------------------------------------------- */
     230               6 :     poDS->SetDescription( poOpenInfo->pszFilename );
     231               6 :     poDS->TryLoadXML();
     232                 : 
     233                 : /* -------------------------------------------------------------------- */
     234                 : /*      Check for overviews.                                            */
     235                 : /* -------------------------------------------------------------------- */
     236               6 :     poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
     237                 : 
     238               6 :     return( poDS );
     239                 : }
     240                 : 
     241                 : /************************************************************************/
     242                 : /*                          GetGeoTransform()                           */
     243                 : /************************************************************************/
     244                 : 
     245               2 : CPLErr CTable2Dataset::GetGeoTransform( double * padfTransform )
     246                 : 
     247                 : {
     248               2 :     memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
     249               2 :     return CE_None;
     250                 : }
     251                 : 
     252                 : /************************************************************************/
     253                 : /*                          SetGeoTransform()                           */
     254                 : /************************************************************************/
     255                 : 
     256               3 : CPLErr CTable2Dataset::SetGeoTransform( double * padfTransform )
     257                 : 
     258                 : {
     259               3 :     if( eAccess == GA_ReadOnly )
     260                 :     {
     261                 :         CPLError( CE_Failure, CPLE_NoWriteAccess,
     262               0 :                   "Unable to update geotransform on readonly file." ); 
     263               0 :         return CE_Failure;
     264                 :     }
     265                 : 
     266               3 :     if( padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
     267                 :     {
     268                 :         CPLError( CE_Failure, CPLE_AppDefined,
     269               0 :                   "Rotated and sheared geotransforms not supported for CTable2."); 
     270               0 :         return CE_Failure;
     271                 :     }
     272                 : 
     273               3 :     memcpy( adfGeoTransform, padfTransform, sizeof(double)*6 );
     274                 : 
     275                 : /* -------------------------------------------------------------------- */
     276                 : /*      Update grid header.                                             */
     277                 : /* -------------------------------------------------------------------- */
     278                 :     double dfValue;
     279                 :     char   achHeader[160];
     280               3 :     double dfDegToRad = M_PI / 180.0;
     281                 : 
     282                 :     // read grid header
     283               3 :     VSIFSeekL( fpImage, 0, SEEK_SET );
     284               3 :     VSIFReadL( achHeader, 1, sizeof(achHeader), fpImage );
     285                 : 
     286                 :     // lower left origin (longitude, center of pixel, radians)
     287               3 :     dfValue = (adfGeoTransform[0] + adfGeoTransform[1]*0.5) * dfDegToRad;
     288                 :     CPL_LSBPTR64( &dfValue );
     289               3 :     memcpy( achHeader + 96, &dfValue, 8 );
     290                 : 
     291                 :     // lower left origin (latitude, center of pixel, radians)
     292               3 :     dfValue = (adfGeoTransform[3] + adfGeoTransform[5] * (nRasterYSize-0.5)) * dfDegToRad;
     293                 :     CPL_LSBPTR64( &dfValue );
     294               3 :     memcpy( achHeader + 104, &dfValue, 8 );
     295                 : 
     296                 :     // pixel width (radians)
     297               3 :     dfValue = adfGeoTransform[1] * dfDegToRad;
     298                 :     CPL_LSBPTR64( &dfValue );
     299               3 :     memcpy( achHeader + 112, &dfValue, 8 );
     300                 : 
     301                 :     // pixel height (radians)
     302               3 :     dfValue = adfGeoTransform[5] * -1 * dfDegToRad;
     303                 :     CPL_LSBPTR64( &dfValue );
     304               3 :     memcpy( achHeader + 120, &dfValue, 8 );
     305                 : 
     306                 :     // write grid header.
     307               3 :     VSIFSeekL( fpImage, 0, SEEK_SET );
     308               3 :     VSIFWriteL( achHeader, 11, 16, fpImage );
     309                 : 
     310               3 :     return CE_None;
     311                 : }
     312                 : 
     313                 : 
     314                 : /************************************************************************/
     315                 : /*                          GetProjectionRef()                          */
     316                 : /************************************************************************/
     317                 : 
     318               0 : const char *CTable2Dataset::GetProjectionRef()
     319                 : 
     320                 : {
     321               0 :     return SRS_WKT_WGS84;
     322                 : }
     323                 : 
     324                 : /************************************************************************/
     325                 : /*                               Create()                               */
     326                 : /************************************************************************/
     327                 : 
     328              43 : GDALDataset *CTable2Dataset::Create( const char * pszFilename,
     329                 :                                      int nXSize, int nYSize, int nBands,
     330                 :                                      GDALDataType eType,
     331                 :                                      char ** papszOptions )
     332                 : 
     333                 : {
     334              43 :     if( eType != GDT_Float32 )
     335                 :     {
     336                 :         CPLError(CE_Failure, CPLE_AppDefined,
     337                 :                  "Attempt to create CTable2 file with unsupported data type '%s'.",
     338              40 :                  GDALGetDataTypeName( eType ) );
     339              40 :         return NULL;
     340                 :     }
     341                 : 
     342                 : /* -------------------------------------------------------------------- */
     343                 : /*      Try to open or create file.                                     */
     344                 : /* -------------------------------------------------------------------- */
     345                 :     VSILFILE  *fp;
     346                 : 
     347               3 :     fp = VSIFOpenL( pszFilename, "wb" );
     348                 :     
     349               3 :     if( fp == NULL )
     350                 :     {
     351                 :         CPLError( CE_Failure, CPLE_OpenFailed,
     352                 :                   "Attempt to create file `%s' failed.\n",
     353               0 :                   pszFilename );
     354               0 :         return NULL;
     355                 :     }
     356                 : 
     357                 : /* -------------------------------------------------------------------- */
     358                 : /*      Create a file header, with a defaulted georeferencing.          */
     359                 : /* -------------------------------------------------------------------- */
     360                 :     char achHeader[160];
     361                 :     int nValue32;
     362                 :     double dfValue;
     363                 : 
     364               3 :     memset( achHeader, 0, sizeof(achHeader));
     365                 : 
     366               3 :     memcpy( achHeader+0, "CTABLE V2.0     ", 16 );
     367                 :     
     368               3 :     if( CSLFetchNameValue( papszOptions, "DESCRIPTION" ) != NULL )
     369                 :         strncpy( achHeader + 16, 
     370                 :                  CSLFetchNameValue( papszOptions, "DESCRIPTION" ), 
     371               0 :                  80 );
     372                 :     
     373                 :     // lower left origin (longitude, center of pixel, radians)
     374               3 :     dfValue = 0;
     375                 :     CPL_LSBPTR64( &dfValue );
     376               3 :     memcpy( achHeader + 96, &dfValue, 8 );
     377                 : 
     378                 :     // lower left origin (latitude, center of pixel, radians)
     379               3 :     dfValue = 0;
     380                 :     CPL_LSBPTR64( &dfValue );
     381               3 :     memcpy( achHeader + 104, &dfValue, 8 );
     382                 : 
     383                 :     // pixel width (radians)
     384               3 :     dfValue = 0.01 * M_PI / 180.0;
     385                 :     CPL_LSBPTR64( &dfValue );
     386               3 :     memcpy( achHeader + 112, &dfValue, 8 );
     387                 : 
     388                 :     // pixel height (radians)
     389               3 :     dfValue = 0.01 * M_PI / 180.0;
     390                 :     CPL_LSBPTR64( &dfValue );
     391               3 :     memcpy( achHeader + 120, &dfValue, 8 );
     392                 : 
     393                 :     // raster width in pixels
     394               3 :     nValue32 = nXSize;
     395                 :     CPL_LSBPTR32( &nValue32 );
     396               3 :     memcpy( achHeader + 128, &nValue32, 4 );
     397                 : 
     398                 :     // raster width in pixels
     399               3 :     nValue32 = nYSize;
     400                 :     CPL_LSBPTR32( &nValue32 );
     401               3 :     memcpy( achHeader + 132, &nValue32, 4 );
     402                 : 
     403               3 :     VSIFWriteL( achHeader, 1, sizeof(achHeader), fp );
     404                 : 
     405                 : /* -------------------------------------------------------------------- */
     406                 : /*      Write zeroed grid data.                                         */
     407                 : /* -------------------------------------------------------------------- */
     408               3 :     float *pafLine = (float *) CPLCalloc(sizeof(float)*2,nXSize);
     409                 :     int i;
     410                 : 
     411             213 :     for( i = 0; i < nYSize; i++ )
     412                 :     {
     413             210 :         if( (int)VSIFWriteL( pafLine, sizeof(float)*2, nXSize, fp ) != nXSize ) 
     414                 :         {
     415                 :             CPLError( CE_Failure, CPLE_FileIO, 
     416                 :                       "Write failed at line %d, perhaps the disk is full?",
     417               0 :                       i );
     418               0 :             return NULL;
     419                 :         }
     420                 :     }
     421                 :     
     422                 : /* -------------------------------------------------------------------- */
     423                 : /*      Cleanup and return.                                             */
     424                 : /* -------------------------------------------------------------------- */
     425               3 :     CPLFree( pafLine );
     426                 : 
     427               3 :     VSIFCloseL( fp );
     428                 : 
     429               3 :     return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
     430                 : }
     431                 : 
     432                 : /************************************************************************/
     433                 : /*                         GDALRegister_CTable2()                          */
     434                 : /************************************************************************/
     435                 : 
     436             582 : void GDALRegister_CTable2()
     437                 : 
     438                 : {
     439                 :     GDALDriver  *poDriver;
     440                 : 
     441             582 :     if( GDALGetDriverByName( "CTable2" ) == NULL )
     442                 :     {
     443             561 :         poDriver = new GDALDriver();
     444                 :         
     445             561 :         poDriver->SetDescription( "CTable2" );
     446                 :         poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
     447             561 :                                    "CTable2 Datum Grid Shift" );
     448             561 :         poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
     449                 : 
     450                 :         poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, 
     451             561 :                                    "Float32" );
     452                 : 
     453             561 :         poDriver->pfnOpen = CTable2Dataset::Open;
     454             561 :         poDriver->pfnIdentify = CTable2Dataset::Identify;
     455             561 :         poDriver->pfnCreate = CTable2Dataset::Create;
     456                 : 
     457             561 :         GetGDALDriverManager()->RegisterDriver( poDriver );
     458                 :     }
     459             582 : }
     460                 : 

Generated by: LCOV version 1.7