LCOV - code coverage report
Current view: directory - frmts/netcdf - gmtdataset.cpp (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 225 180 80.0 %
Date: 2013-03-30 Functions: 13 9 69.2 %

       1                 : /******************************************************************************
       2                 :  * $Id: gmtdataset.cpp 25593 2013-02-02 13:03:03Z rouault $
       3                 :  *
       4                 :  * Project:  netCDF read/write Driver
       5                 :  * Purpose:  GDAL bindings over netCDF library for GMT Grids.
       6                 :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       7                 :  *
       8                 :  ******************************************************************************
       9                 :  * Copyright (c) 2004, Frank Warmerdam
      10                 :  *
      11                 :  * Permission is hereby granted, free of charge, to any person obtaining a
      12                 :  * copy of this software and associated documentation files (the "Software"),
      13                 :  * to deal in the Software without restriction, including without limitation
      14                 :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      15                 :  * and/or sell copies of the Software, and to permit persons to whom the
      16                 :  * Software is furnished to do so, subject to the following conditions:
      17                 :  *
      18                 :  * The above copyright notice and this permission notice shall be included
      19                 :  * in all copies or substantial portions of the Software.
      20                 :  *
      21                 :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      22                 :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      23                 :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      24                 :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      25                 :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      26                 :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      27                 :  * DEALINGS IN THE SOFTWARE.
      28                 :  ****************************************************************************/
      29                 : 
      30                 : #include "gdal_pam.h"
      31                 : #include "gdal_frmts.h"
      32                 : #include "netcdf.h"
      33                 : #include "cpl_multiproc.h"
      34                 : 
      35                 : CPL_CVSID("$Id: gmtdataset.cpp 25593 2013-02-02 13:03:03Z rouault $");
      36                 : 
      37                 : extern void *hNCMutex; /* shared with netcdf. See netcdfdataset.cpp */
      38                 : 
      39                 : /************************************************************************/
      40                 : /* ==================================================================== */
      41                 : /*           GMTDataset               */
      42                 : /* ==================================================================== */
      43                 : /************************************************************************/
      44                 : 
      45                 : class GMTRasterBand;
      46                 : 
      47                 : class GMTDataset : public GDALPamDataset
      48               8 : {
      49                 :     int         z_id;
      50                 :     double      adfGeoTransform[6];
      51                 : 
      52                 :   public:
      53                 :     int         cdfid;
      54                 : 
      55                 :     ~GMTDataset();
      56                 :     
      57                 :     static GDALDataset *Open( GDALOpenInfo * );
      58                 : 
      59                 :     CPLErr  GetGeoTransform( double * padfTransform );
      60                 : };
      61                 : 
      62                 : /************************************************************************/
      63                 : /* ==================================================================== */
      64                 : /*                         GMTRasterBand                                */
      65                 : /* ==================================================================== */
      66                 : /************************************************************************/
      67                 : 
      68                 : class GMTRasterBand : public GDALPamRasterBand
      69               8 : {
      70                 :     nc_type nc_datatype;
      71                 :     int         nZId;
      72                 : 
      73                 :   public:
      74                 : 
      75                 :         GMTRasterBand( GMTDataset *poDS, int nZId, int nBand );
      76                 :     
      77                 :     virtual CPLErr IReadBlock( int, int, void * );
      78                 : };
      79                 : 
      80                 : 
      81                 : /************************************************************************/
      82                 : /*                           GMTRasterBand()                            */
      83                 : /************************************************************************/
      84                 : 
      85               8 : GMTRasterBand::GMTRasterBand( GMTDataset *poDS, int nZId, int nBand )
      86                 : 
      87                 : {
      88               8 :     this->poDS = poDS;
      89               8 :     this->nBand = nBand;
      90               8 :     this->nZId = nZId;
      91                 :     
      92               8 :     nBlockXSize = poDS->GetRasterXSize();
      93               8 :     nBlockYSize = 1;
      94                 : 
      95                 : /* -------------------------------------------------------------------- */
      96                 : /*      Get the type of the "z" variable, our target raster array.      */
      97                 : /* -------------------------------------------------------------------- */
      98               8 :     if( nc_inq_var( poDS->cdfid, nZId, NULL, &nc_datatype, NULL, NULL,
      99                 :                     NULL ) != NC_NOERR )
     100                 :     {
     101                 :         CPLError( CE_Failure, CPLE_AppDefined, 
     102               0 :                   "Error in nc_var_inq() on 'z'." );
     103               0 :         return;
     104                 :     }
     105                 : 
     106               8 :     if( nc_datatype == NC_BYTE )
     107               0 :         eDataType = GDT_Byte;
     108               8 :     else if( nc_datatype == NC_SHORT )
     109               4 :         eDataType = GDT_Int16;
     110               4 :     else if( nc_datatype == NC_INT )
     111               1 :         eDataType = GDT_Int32;
     112               3 :     else if( nc_datatype == NC_FLOAT )
     113               2 :         eDataType = GDT_Float32;
     114               1 :     else if( nc_datatype == NC_DOUBLE )
     115               1 :         eDataType = GDT_Float64;
     116                 :     else
     117                 :     {
     118               0 :         if( nBand == 1 )
     119                 :             CPLError( CE_Warning, CPLE_AppDefined, 
     120                 :                       "Unsupported GMT datatype (%d), treat as Float32.", 
     121               0 :                       (int) nc_datatype );
     122               0 :         eDataType = GDT_Float32;
     123                 :     }
     124               0 : }
     125                 : 
     126                 : /************************************************************************/
     127                 : /*                             IReadBlock()                             */
     128                 : /************************************************************************/
     129                 : 
     130              90 : CPLErr GMTRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
     131                 :                                   void * pImage )
     132                 : 
     133                 : {
     134                 :     size_t start[2], edge[2];
     135              90 :     int    nErr = NC_NOERR;
     136              90 :     int    cdfid = ((GMTDataset *) poDS)->cdfid;
     137                 :     
     138              90 :     CPLMutexHolderD(&hNCMutex);
     139                 : 
     140              90 :     start[0] = nBlockYOff * nBlockXSize;
     141              90 :     edge[0] = nBlockXSize;
     142                 : 
     143              90 :     if( eDataType == GDT_Byte )
     144                 :         nErr = nc_get_vara_uchar( cdfid, nZId, start, edge, 
     145               0 :                                   (unsigned char *) pImage );
     146              90 :     else if( eDataType == GDT_Int16 )
     147                 :         nErr = nc_get_vara_short( cdfid, nZId, start, edge, 
     148              40 :                                   (short int *) pImage );
     149              50 :     else if( eDataType == GDT_Int32 )
     150                 :     {
     151                 :         if( sizeof(long) == 4 )
     152                 :             nErr = nc_get_vara_long( cdfid, nZId, start, edge, 
     153                 :                                      (long *) pImage );
     154                 :         else
     155                 :             nErr = nc_get_vara_int( cdfid, nZId, start, edge, 
     156               0 :                                     (int *) pImage );
     157                 :     }
     158              50 :     else if( eDataType == GDT_Float32 )
     159                 :         nErr = nc_get_vara_float( cdfid, nZId, start, edge, 
     160              50 :                                   (float *) pImage );
     161               0 :     else if( eDataType == GDT_Float64 )
     162                 :         nErr = nc_get_vara_double( cdfid, nZId, start, edge, 
     163               0 :                                    (double *) pImage );
     164                 : 
     165              90 :     if( nErr != NC_NOERR )
     166                 :     {
     167                 :         CPLError( CE_Failure, CPLE_AppDefined, 
     168                 :                   "GMT scanline fetch failed: %s", 
     169               0 :                   nc_strerror( nErr ) );
     170               0 :         return CE_Failure;
     171                 :     }
     172                 :     else
     173              90 :         return CE_None;
     174                 : }
     175                 : 
     176                 : /************************************************************************/
     177                 : /* ==================================================================== */
     178                 : /*        GMTDataset        */
     179                 : /* ==================================================================== */
     180                 : /************************************************************************/
     181                 : 
     182                 : /************************************************************************/
     183                 : /*                            ~GMTDataset()                             */
     184                 : /************************************************************************/
     185                 : 
     186               8 : GMTDataset::~GMTDataset()
     187                 : 
     188                 : {
     189               8 :     FlushCache();
     190               8 :     nc_close (cdfid);
     191               8 : }
     192                 : 
     193                 : /************************************************************************/
     194                 : /*                          GetGeoTransform()                           */
     195                 : /************************************************************************/
     196                 : 
     197               7 : CPLErr GMTDataset::GetGeoTransform( double * padfTransform )
     198                 : 
     199                 : {
     200               7 :     memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
     201               7 :     return CE_None;
     202                 : }
     203                 : 
     204                 : /************************************************************************/
     205                 : /*                                Open()                                */
     206                 : /************************************************************************/
     207                 : 
     208           13567 : GDALDataset *GMTDataset::Open( GDALOpenInfo * poOpenInfo )
     209                 : 
     210                 : {
     211                 : /* -------------------------------------------------------------------- */
     212                 : /*      Does this file have the GMT magic number?                    */
     213                 : /* -------------------------------------------------------------------- */
     214           13567 :     if( poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 50 )
     215           12039 :         return NULL;
     216                 : 
     217            2166 :     if( poOpenInfo->pabyHeader[0] != 'C' 
     218             220 :         || poOpenInfo->pabyHeader[1] != 'D' 
     219             209 :         || poOpenInfo->pabyHeader[2] != 'F' 
     220             209 :         || poOpenInfo->pabyHeader[3] != 1 )
     221            1320 :         return NULL;
     222                 :     
     223             208 :     CPLMutexHolderD(&hNCMutex);
     224                 : 
     225                 : /* -------------------------------------------------------------------- */
     226                 : /*      Try opening the dataset.                                        */
     227                 : /* -------------------------------------------------------------------- */
     228                 :     int cdfid, nm_id, dim_count, z_id;
     229                 : 
     230             208 :     if( nc_open( poOpenInfo->pszFilename, NC_NOWRITE, &cdfid ) != NC_NOERR )
     231               0 :         return NULL;
     232                 : 
     233             208 :     if( nc_inq_varid( cdfid, "dimension", &nm_id ) != NC_NOERR 
     234                 :         || nc_inq_varid( cdfid, "z", &z_id ) != NC_NOERR )
     235                 :     {
     236                 : #ifdef notdef
     237                 :         CPLError( CE_Warning, CPLE_AppDefined, 
     238                 :                   "%s is a GMT file, but not in GMT configuration.",
     239                 :                   poOpenInfo->pszFilename );
     240                 : #endif
     241             200 :         nc_close( cdfid );
     242             200 :         return NULL;
     243                 :     }
     244                 : 
     245               8 :     if( nc_inq_ndims( cdfid, &dim_count ) != NC_NOERR || dim_count < 2 )
     246                 :     {
     247               0 :         nc_close( cdfid );
     248               0 :         return NULL;
     249                 :     }
     250                 :     
     251                 : /* -------------------------------------------------------------------- */
     252                 : /*      Confirm the requested access is supported.                      */
     253                 : /* -------------------------------------------------------------------- */
     254               8 :     if( poOpenInfo->eAccess == GA_Update )
     255                 :     {
     256               0 :         nc_close( cdfid );
     257                 :         CPLError( CE_Failure, CPLE_NotSupported, 
     258                 :                   "The GMT driver does not support update access to existing"
     259               0 :                   " datasets.\n" );
     260               0 :         return NULL;
     261                 :     }
     262                 :     
     263                 : /* -------------------------------------------------------------------- */
     264                 : /*      Create a corresponding GDALDataset.                             */
     265                 : /* -------------------------------------------------------------------- */
     266                 :     GMTDataset  *poDS;
     267                 : 
     268               8 :     CPLReleaseMutex(hNCMutex);  // Release mutex otherwise we'll deadlock with GDALDataset own mutex
     269               8 :     poDS = new GMTDataset();
     270               8 :     CPLAcquireMutex(hNCMutex, 1000.0);
     271                 : 
     272               8 :     poDS->cdfid = cdfid;
     273               8 :     poDS->z_id = z_id;
     274                 :     
     275                 : /* -------------------------------------------------------------------- */
     276                 : /*      Get dimensions.  If we can't find this, then this is a          */
     277                 : /*      GMT file, but not a normal grid product.                     */
     278                 : /* -------------------------------------------------------------------- */
     279                 :     size_t start[2], edge[2];
     280                 :     int    nm[2];
     281                 : 
     282               8 :     start[0] = 0;
     283               8 :     edge[0] = 2;
     284                 : 
     285               8 :     nc_get_vara_int(cdfid, nm_id, start, edge, nm);
     286                 :     
     287               8 :     poDS->nRasterXSize = nm[0];
     288               8 :     poDS->nRasterYSize = nm[1];
     289                 : 
     290                 : /* -------------------------------------------------------------------- */
     291                 : /*      Fetch "z" attributes scale_factor, add_offset, and              */
     292                 : /*      node_offset.                                                    */
     293                 : /* -------------------------------------------------------------------- */
     294               8 :     double scale_factor=1.0, add_offset=0.0;
     295               8 :     int node_offset = 1;
     296                 : 
     297               8 :     nc_get_att_double( cdfid, z_id, "scale_factor", &scale_factor );
     298               8 :     nc_get_att_double( cdfid, z_id, "add_offset", &add_offset );
     299               8 :     nc_get_att_int( cdfid, z_id, "node_offset", &node_offset );
     300                 : 
     301                 : /* -------------------------------------------------------------------- */
     302                 : /*      Get x/y range information.                                      */
     303                 : /* -------------------------------------------------------------------- */
     304                 :     int x_range_id, y_range_id;
     305                 : 
     306               8 :     if( nc_inq_varid (cdfid, "x_range", &x_range_id) == NC_NOERR 
     307                 :         && nc_inq_varid (cdfid, "y_range", &y_range_id) == NC_NOERR )
     308                 :     {
     309                 :         double x_range[2], y_range[2];
     310                 : 
     311               8 :         nc_get_vara_double( cdfid, x_range_id, start, edge, x_range );
     312               8 :         nc_get_vara_double( cdfid, y_range_id, start, edge, y_range );
     313                 : 
     314                 :         // Pixel is area
     315               8 :         if( node_offset == 1 )
     316                 :         {
     317               8 :             poDS->adfGeoTransform[0] = x_range[0];
     318                 :             poDS->adfGeoTransform[1] = 
     319               8 :                 (x_range[1] - x_range[0]) / poDS->nRasterXSize;
     320               8 :             poDS->adfGeoTransform[2] = 0.0;
     321               8 :             poDS->adfGeoTransform[3] = y_range[1];
     322               8 :             poDS->adfGeoTransform[4] = 0.0;
     323                 :             poDS->adfGeoTransform[5] = 
     324               8 :                 (y_range[0] - y_range[1]) / poDS->nRasterYSize;
     325                 :         }
     326                 : 
     327                 :         // Pixel is point - offset by half pixel. 
     328                 :         else /* node_offset == 0 */ 
     329                 :         {
     330                 :             poDS->adfGeoTransform[1] = 
     331               0 :                 (x_range[1] - x_range[0]) / (poDS->nRasterXSize-1);
     332                 :             poDS->adfGeoTransform[0] = 
     333               0 :                 x_range[0] - poDS->adfGeoTransform[1]*0.5;
     334               0 :             poDS->adfGeoTransform[2] = 0.0;
     335               0 :             poDS->adfGeoTransform[4] = 0.0;
     336                 :             poDS->adfGeoTransform[5] = 
     337               0 :                 (y_range[0] - y_range[1]) / (poDS->nRasterYSize-1);
     338                 :             poDS->adfGeoTransform[3] = 
     339               0 :                 y_range[1] - poDS->adfGeoTransform[5]*0.5;
     340                 :         }
     341                 :     }
     342                 :     else
     343                 :     {
     344               0 :         poDS->adfGeoTransform[0] = 0.0;
     345               0 :         poDS->adfGeoTransform[1] = 1.0;
     346               0 :         poDS->adfGeoTransform[2] = 0.0;
     347               0 :         poDS->adfGeoTransform[3] = 0.0;
     348               0 :         poDS->adfGeoTransform[4] = 0.0;
     349               0 :         poDS->adfGeoTransform[5] = 1.0;
     350                 :     }
     351                 : 
     352                 : /* -------------------------------------------------------------------- */
     353                 : /*      Create band information objects.                                */
     354                 : /* -------------------------------------------------------------------- */
     355               8 :     poDS->nBands = 1;
     356               8 :     poDS->SetBand( 1, new GMTRasterBand( poDS, z_id, 1 ));
     357                 : 
     358               8 :     if( scale_factor != 1.0 || add_offset != 0.0 )
     359                 :     {
     360               0 :         poDS->GetRasterBand(1)->SetOffset( add_offset );
     361               0 :         poDS->GetRasterBand(1)->SetScale( scale_factor );
     362                 :     }
     363                 : 
     364                 : /* -------------------------------------------------------------------- */
     365                 : /*      Initialize any PAM information.                                 */
     366                 : /* -------------------------------------------------------------------- */
     367               8 :     poDS->SetDescription( poOpenInfo->pszFilename );
     368                 : 
     369               8 :     CPLReleaseMutex(hNCMutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
     370               8 :     poDS->TryLoadXML();
     371                 : 
     372                 : /* -------------------------------------------------------------------- */
     373                 : /*      Check for external overviews.                                   */
     374                 : /* -------------------------------------------------------------------- */
     375               8 :     poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename, poOpenInfo->papszSiblingFiles );
     376               8 :     CPLAcquireMutex(hNCMutex, 1000.0);
     377                 : 
     378               8 :     return( poDS );
     379                 : }
     380                 : 
     381                 : /************************************************************************/
     382                 : /*                           GMTCreateCopy()                            */
     383                 : /*                                                                      */
     384                 : /*      This code mostly cribbed from GMT's "gmt_cdf.c" module.         */
     385                 : /************************************************************************/
     386                 : 
     387                 : static GDALDataset *
     388              19 : GMTCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
     389                 :                   int bStrict, char ** papszOptions, 
     390                 :                   GDALProgressFunc pfnProgress, void * pProgressData )
     391                 : 
     392                 : {
     393                 : /* -------------------------------------------------------------------- */
     394                 : /*      Figure out general characteristics.                             */
     395                 : /* -------------------------------------------------------------------- */
     396                 :     nc_type nc_datatype;
     397                 :     GDALRasterBand *poBand;
     398                 :     int nXSize, nYSize;
     399                 :     
     400              19 :     CPLMutexHolderD(&hNCMutex);
     401                 : 
     402              19 :     if( poSrcDS->GetRasterCount() != 1 )
     403                 :     {
     404                 :         CPLError( CE_Failure, CPLE_AppDefined, 
     405               5 :                   "Currently GMT export only supports 1 band datasets." );
     406               5 :         return NULL;
     407                 :     }
     408                 : 
     409              14 :     poBand = poSrcDS->GetRasterBand(1);
     410                 : 
     411              14 :     nXSize = poSrcDS->GetRasterXSize();
     412              14 :     nYSize = poSrcDS->GetRasterYSize();
     413                 :     
     414              14 :     if( poBand->GetRasterDataType() == GDT_Int16 )
     415               2 :         nc_datatype = NC_SHORT;
     416              12 :     else if( poBand->GetRasterDataType() == GDT_Int32 )
     417               1 :         nc_datatype = NC_INT;
     418              11 :     else if( poBand->GetRasterDataType() == GDT_Float32 )
     419               1 :         nc_datatype = NC_FLOAT;
     420              10 :     else if( poBand->GetRasterDataType() == GDT_Float64 )
     421               1 :         nc_datatype = NC_DOUBLE;
     422               9 :     else if( bStrict )
     423                 :     {
     424                 :         CPLError( CE_Failure, CPLE_AppDefined, 
     425                 :                   "Band data type %s not supported in GMT, giving up.",
     426               8 :                   GDALGetDataTypeName( poBand->GetRasterDataType() ) );
     427               8 :         return NULL;
     428                 :     }
     429               1 :     else if( poBand->GetRasterDataType() == GDT_Byte )
     430               1 :         nc_datatype = NC_SHORT;
     431               0 :     else if( poBand->GetRasterDataType() == GDT_UInt16 )
     432               0 :         nc_datatype = NC_INT;
     433               0 :     else if( poBand->GetRasterDataType() == GDT_UInt32 )
     434               0 :         nc_datatype = NC_INT;
     435                 :     else 
     436               0 :         nc_datatype = NC_FLOAT;
     437                 :     
     438                 : /* -------------------------------------------------------------------- */
     439                 : /*      Establish bounds from geotransform.                             */
     440                 : /* -------------------------------------------------------------------- */
     441                 :     double adfGeoTransform[6];
     442                 :     double dfXMax, dfYMin;
     443                 : 
     444               6 :     poSrcDS->GetGeoTransform( adfGeoTransform );
     445                 :     
     446               6 :     if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 )
     447                 :     {
     448                 :         CPLError( bStrict ? CE_Failure : CE_Warning, CPLE_AppDefined, 
     449               0 :                   "Geotransform has rotational coefficients not supported in GMT." );
     450               0 :         if( bStrict )
     451               0 :             return NULL;
     452                 :     }
     453                 : 
     454               6 :     dfXMax = adfGeoTransform[0] + adfGeoTransform[1] * nXSize;
     455               6 :     dfYMin = adfGeoTransform[3] + adfGeoTransform[5] * nYSize;
     456                 :     
     457                 : /* -------------------------------------------------------------------- */
     458                 : /*      Create base file.                                               */
     459                 : /* -------------------------------------------------------------------- */
     460                 :     int cdfid, err;
     461                 : 
     462               6 :     err = nc_create (pszFilename, NC_CLOBBER,&cdfid);
     463               6 :     if( err != NC_NOERR )
     464                 :     {
     465                 :         CPLError( CE_Failure, CPLE_AppDefined, 
     466                 :                   "nc_create(%s): %s", 
     467               1 :                   pszFilename, nc_strerror( err ) );
     468               1 :         return NULL;
     469                 :     }
     470                 : 
     471                 : /* -------------------------------------------------------------------- */
     472                 : /*      Define the dimensions and so forth.                             */
     473                 : /* -------------------------------------------------------------------- */
     474                 :     int side_dim, xysize_dim, dims[1];
     475                 :     int x_range_id, y_range_id, z_range_id, inc_id, nm_id, z_id;
     476                 : 
     477               5 :     nc_def_dim(cdfid, "side", 2, &side_dim);
     478               5 :     nc_def_dim(cdfid, "xysize", (int) (nXSize * nYSize), &xysize_dim);
     479                 : 
     480               5 :     dims[0]   = side_dim;
     481               5 :     nc_def_var (cdfid, "x_range", NC_DOUBLE, 1, dims, &x_range_id);
     482               5 :     nc_def_var (cdfid, "y_range", NC_DOUBLE, 1, dims, &y_range_id);
     483               5 :     nc_def_var (cdfid, "z_range", NC_DOUBLE, 1, dims, &z_range_id);
     484               5 :     nc_def_var (cdfid, "spacing", NC_DOUBLE, 1, dims, &inc_id);
     485               5 :     nc_def_var (cdfid, "dimension", NC_LONG, 1, dims, &nm_id);
     486                 : 
     487               5 :     dims[0]   = xysize_dim;
     488               5 :     nc_def_var (cdfid, "z", nc_datatype, 1, dims, &z_id);
     489                 : 
     490                 : /* -------------------------------------------------------------------- */
     491                 : /*      Assign attributes.                                              */
     492                 : /* -------------------------------------------------------------------- */
     493               5 :     double default_scale = 1.0;
     494               5 :     double default_offset = 0.0;
     495               5 :     int default_node_offset = 1; // pixel is area
     496                 : 
     497               5 :     nc_put_att_text (cdfid, x_range_id, "units", 7, "meters");
     498               5 :     nc_put_att_text (cdfid, y_range_id, "units", 7, "meters");
     499               5 :     nc_put_att_text (cdfid, z_range_id, "units", 7, "meters");
     500                 : 
     501                 :     nc_put_att_double (cdfid, z_id, "scale_factor", NC_DOUBLE, 1, 
     502               5 :                        &default_scale );
     503                 :     nc_put_att_double (cdfid, z_id, "add_offset", NC_DOUBLE, 1, 
     504               5 :                        &default_offset );
     505                 : 
     506                 :     nc_put_att_int (cdfid, z_id, "node_offset", NC_LONG, 1, 
     507               5 :                     &default_node_offset );
     508               5 :     nc_put_att_text (cdfid, NC_GLOBAL, "title", 1, "");
     509               5 :     nc_put_att_text (cdfid, NC_GLOBAL, "source", 1, "");
     510                 :   
     511                 :     /* leave define mode */
     512               5 :     nc_enddef (cdfid);
     513                 : 
     514                 : /* -------------------------------------------------------------------- */
     515                 : /*      Get raster min/max.                                             */
     516                 : /* -------------------------------------------------------------------- */
     517                 :     double adfMinMax[2];
     518               5 :     GDALComputeRasterMinMax( (GDALRasterBandH) poBand, FALSE, adfMinMax );
     519                 :   
     520                 : /* -------------------------------------------------------------------- */
     521                 : /*      Set range variables.                                            */
     522                 : /* -------------------------------------------------------------------- */
     523                 :     size_t start[2], edge[2];
     524                 :     double dummy[2];
     525                 :     int nm[2];
     526                 :   
     527               5 :     start[0] = 0;
     528               5 :     edge[0] = 2;
     529               5 :     dummy[0] = adfGeoTransform[0];
     530               5 :     dummy[1] = dfXMax;
     531               5 :     nc_put_vara_double(cdfid, x_range_id, start, edge, dummy);
     532                 : 
     533               5 :     dummy[0] = dfYMin;
     534               5 :     dummy[1] = adfGeoTransform[3];
     535               5 :     nc_put_vara_double(cdfid, y_range_id, start, edge, dummy);
     536                 : 
     537               5 :     dummy[0] = adfGeoTransform[1];
     538               5 :     dummy[1] = -adfGeoTransform[5];
     539               5 :     nc_put_vara_double(cdfid, inc_id, start, edge, dummy);
     540                 : 
     541               5 :     nm[0] = nXSize;
     542               5 :     nm[1] = nYSize;
     543               5 :     nc_put_vara_int(cdfid, nm_id, start, edge, nm);
     544                 : 
     545               5 :     nc_put_vara_double(cdfid, z_range_id, start, edge, adfMinMax);
     546                 : 
     547                 : /* -------------------------------------------------------------------- */
     548                 : /*      Write out the image one scanline at a time.                     */
     549                 : /* -------------------------------------------------------------------- */
     550                 :     double *padfData;
     551                 :     int  iLine;
     552                 : 
     553               5 :     padfData = (double *) CPLMalloc( sizeof(double) * nXSize );
     554                 : 
     555               5 :     edge[0] = nXSize;
     556              65 :     for( iLine = 0; iLine < nYSize; iLine++ )
     557                 :     {
     558              60 :         start[0] = iLine * nXSize;
     559                 :         poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, 
     560              60 :                           padfData, nXSize, 1, GDT_Float64, 0, 0 );
     561              60 :         err = nc_put_vara_double( cdfid, z_id, start, edge, padfData );
     562              60 :         if( err != NC_NOERR )
     563                 :         {
     564                 :             CPLError( CE_Failure, CPLE_AppDefined, 
     565                 :                       "nc_put_vara_double(%s): %s", 
     566               0 :                       pszFilename, nc_strerror( err ) );
     567               0 :             nc_close (cdfid);
     568               0 :             return( NULL );
     569                 :         }
     570                 :     }
     571                 :     
     572               5 :     CPLFree( padfData );
     573                 : 
     574                 : /* -------------------------------------------------------------------- */
     575                 : /*      Close file, and reopen.                                         */
     576                 : /* -------------------------------------------------------------------- */
     577               5 :     nc_close (cdfid);
     578                 : 
     579                 : /* -------------------------------------------------------------------- */
     580                 : /*      Re-open dataset, and copy any auxilary pam information.         */
     581                 : /* -------------------------------------------------------------------- */
     582                 :     GDALPamDataset *poDS = (GDALPamDataset *) 
     583               5 :         GDALOpen( pszFilename, GA_ReadOnly );
     584                 : 
     585               5 :     if( poDS )
     586               5 :         poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT );
     587                 : 
     588               5 :     return poDS;
     589                 : }
     590                 : 
     591                 : /************************************************************************/
     592                 : /*                          GDALRegister_GMT()                          */
     593                 : /************************************************************************/
     594                 : 
     595             610 : void GDALRegister_GMT()
     596                 : 
     597                 : {
     598                 :     GDALDriver  *poDriver;
     599                 :     
     600             610 :     if (! GDAL_CHECK_VERSION("GMT driver"))
     601               0 :         return;
     602                 : 
     603             610 :     if( GDALGetDriverByName( "GMT" ) == NULL )
     604                 :     {
     605             588 :         poDriver = new GDALDriver();
     606                 :         
     607             588 :         poDriver->SetDescription( "GMT" );
     608                 :         poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
     609             588 :                                    "GMT NetCDF Grid Format" );
     610                 :         poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, 
     611             588 :                                    "frmt_various.html#GMT" );
     612             588 :         poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nc" );
     613                 : 
     614                 :         poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, 
     615             588 :                                    "Int16 Int32 Float32 Float64" );
     616                 : 
     617             588 :         poDriver->pfnOpen = GMTDataset::Open;
     618             588 :         poDriver->pfnCreateCopy = GMTCreateCopy;
     619                 : 
     620             588 :         GetGDALDriverManager()->RegisterDriver( poDriver );
     621                 :     }
     622                 : }

Generated by: LCOV version 1.7