LTP GCOV extension - code coverage report
Current view: directory - apps - gdalbuildvrt.cpp
Test: gdal_filtered.info
Date: 2010-07-12 Instrumented lines: 514
Code covered: 77.6 % Executed lines: 399

       1                 : /******************************************************************************
       2                 :  * $Id: gdalbuildvrt.cpp 19804 2010-06-05 11:33:05Z rouault $
       3                 :  *
       4                 :  * Project:  GDAL Utilities
       5                 :  * Purpose:  Commandline application to build VRT datasets from raster products or content of SHP tile index
       6                 :  * Author:   Even Rouault, even.rouault at mines-paris.org
       7                 :  *
       8                 :  ******************************************************************************
       9                 :  * Copyright (c) 2007, Even Rouault
      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_proxy.h"
      31                 : #include "cpl_string.h"
      32                 : #include "vrt/gdal_vrt.h"
      33                 : 
      34                 : #ifdef OGR_ENABLED
      35                 : #include "ogr_api.h"
      36                 : #endif
      37                 : 
      38                 : CPL_CVSID("$Id: gdalbuildvrt.cpp 19804 2010-06-05 11:33:05Z rouault $");
      39                 : 
      40                 : #define GEOTRSFRM_TOPLEFT_X            0
      41                 : #define GEOTRSFRM_WE_RES               1
      42                 : #define GEOTRSFRM_ROTATION_PARAM1      2
      43                 : #define GEOTRSFRM_TOPLEFT_Y            3
      44                 : #define GEOTRSFRM_ROTATION_PARAM2      4
      45                 : #define GEOTRSFRM_NS_RES               5
      46                 : 
      47                 : typedef enum
      48                 : {
      49                 :     LOWEST_RESOLUTION,
      50                 :     HIGHEST_RESOLUTION,
      51                 :     AVERAGE_RESOLUTION,
      52                 :     USER_RESOLUTION
      53                 : } ResolutionStrategy;
      54                 : 
      55                 : typedef struct
      56                 : {
      57                 :     int    isFileOK;
      58                 :     int    nRasterXSize;
      59                 :     int    nRasterYSize;
      60                 :     double adfGeoTransform[6];
      61                 :     int    nBlockXSize;
      62                 :     int    nBlockYSize;
      63                 :     GDALDataType firstBandType;
      64                 :     int*         panHasNoData;
      65                 :     double*      padfNoDataValues;
      66                 : } DatasetProperty;
      67                 : 
      68                 : typedef struct
      69                 : {
      70                 :     GDALColorInterp        colorInterpretation;
      71                 :     GDALDataType           dataType;
      72                 :     GDALColorTableH        colorTable;
      73                 :     int                    bHasNoData;
      74                 :     double                 noDataValue;
      75                 : } BandProperty;
      76                 : 
      77                 : /************************************************************************/
      78                 : /*                               Usage()                                */
      79                 : /************************************************************************/
      80                 : 
      81               0 : static void Usage()
      82                 : 
      83                 : {
      84                 :     fprintf(stdout, "%s", 
      85                 :             "Usage: gdalbuildvrt [-tileindex field_name] [-resolution {highest|lowest|average|user}]\n"
      86                 :             "                    [-tr xres yres] [-separate] [-allow_projection_difference] [-q]\n"
      87                 :             "                    [-te xmin ymin xmax ymax] [-addalpha] [-hidenodata] \n"
      88                 :             "                    [-srcnodata \"value [value...]\"] [-vrtnodata \"value [value...]\"] \n"
      89                 :             "                    [-input_file_list my_liste.txt] output.vrt [gdalfile]*\n"
      90                 :             "\n"
      91                 :             "eg.\n"
      92                 :             "  % gdalbuildvrt doq_index.vrt doq/*.tif\n"
      93                 :             "  % gdalbuildvrt -input_file_list my_liste.txt doq_index.vrt\n"
      94                 :             "\n"
      95                 :             "NOTES:\n"
      96                 :             "  o With -separate, each files goes into a separate band in the VRT band. Otherwise,\n"
      97                 :             "    the files are considered as tiles of a larger mosaic.\n"
      98                 :             "  o The default tile index field is 'location' unless otherwise specified by -tileindex.\n"
      99                 :             "  o In case the resolution of all input files is not the same, the -resolution flag.\n"
     100                 :             "    enable the user to control the way the output resolution is computed. average is the default.\n"
     101                 :             "  o Input files may be any valid GDAL dataset or a GDAL raster tile index.\n"
     102                 :             "  o For a GDAL raster tile index, all entries will be added to the VRT.\n"
     103                 :             "  o If one GDAL dataset is made of several subdatasets and has 0 raster bands, its\n"
     104                 :             "    datasets will be added to the VRT rather than the dataset itself.\n"
     105                 :             "  o By default, only datasets of same projection and band characteristics may be added to the VRT.\n"
     106               0 :             );
     107               0 :     exit( 1 );
     108                 : }
     109                 : 
     110                 : 
     111                 : /************************************************************************/
     112                 : /*                         GetSrcDstWin()                               */
     113                 : /************************************************************************/
     114                 : 
     115                 : int  GetSrcDstWin(DatasetProperty* psDP,
     116                 :                   double we_res, double ns_res,
     117                 :                   double minX, double minY, double maxX, double maxY,
     118                 :                   int* pnSrcXOff, int* pnSrcYOff, int* pnSrcXSize, int* pnSrcYSize,
     119              38 :                   int* pnDstXOff, int* pnDstYOff, int* pnDstXSize, int* pnDstYSize)
     120                 : {
     121                 :     /* Check that the destination bounding box intersects the source bounding box */
     122              38 :     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] +
     123                 :          psDP->nRasterXSize *
     124                 :          psDP->adfGeoTransform[GEOTRSFRM_WE_RES] < minX )
     125               0 :          return FALSE;
     126              38 :     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] > maxX )
     127               0 :          return FALSE;
     128              38 :     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] +
     129                 :          psDP->nRasterYSize *
     130                 :          psDP->adfGeoTransform[GEOTRSFRM_NS_RES] > maxY )
     131               0 :          return FALSE;
     132              38 :     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] < minY )
     133               0 :          return FALSE;
     134                 : 
     135              38 :     *pnSrcXSize = psDP->nRasterXSize;
     136              38 :     *pnSrcYSize = psDP->nRasterYSize;
     137              38 :     if ( psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] < minX )
     138                 :     {
     139                 :         *pnSrcXOff = (int)((minX - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X]) /
     140               1 :             psDP->adfGeoTransform[GEOTRSFRM_WE_RES] + 0.5);
     141               1 :         *pnDstXOff = 0;
     142                 :     }
     143                 :     else
     144                 :     {
     145              37 :         *pnSrcXOff = 0;
     146                 :         *pnDstXOff = (int)
     147              37 :             (0.5 + (psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_X] - minX) / we_res);
     148                 :     }
     149              38 :     if ( maxY < psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y])
     150                 :     {
     151                 :         *pnSrcYOff = (int)((psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] - maxY) /
     152               1 :             -psDP->adfGeoTransform[GEOTRSFRM_NS_RES] + 0.5);
     153               1 :         *pnDstYOff = 0;
     154                 :     }
     155                 :     else
     156                 :     {
     157              37 :         *pnSrcYOff = 0;
     158                 :         *pnDstYOff = (int)
     159              37 :             (0.5 + (maxY - psDP->adfGeoTransform[GEOTRSFRM_TOPLEFT_Y]) / -ns_res);
     160                 :     }
     161                 :     *pnDstXSize = (int)
     162                 :         (0.5 + psDP->nRasterXSize *
     163              38 :          psDP->adfGeoTransform[GEOTRSFRM_WE_RES] / we_res);
     164                 :     *pnDstYSize = (int)
     165                 :         (0.5 + psDP->nRasterYSize *
     166              38 :          psDP->adfGeoTransform[GEOTRSFRM_NS_RES] / ns_res);
     167                 :          
     168              38 :     return TRUE;
     169                 : }
     170                 : 
     171                 : /************************************************************************/
     172                 : /*                            VRTBuilder                                */
     173                 : /************************************************************************/
     174                 : 
     175                 : class VRTBuilder
     176                 : {
     177                 :     /* Input parameters */
     178                 :     char               *pszOutputFilename;
     179                 :     int                 nInputFiles;
     180                 :     char              **ppszInputFilenames;
     181                 :     ResolutionStrategy  resolutionStrategy;
     182                 :     double              we_res;
     183                 :     double              ns_res;
     184                 :     double              minX;
     185                 :     double              minY;
     186                 :     double              maxX;
     187                 :     double              maxY;
     188                 :     int                 bSeparate;
     189                 :     int                 bAllowProjectionDifference;
     190                 :     int                 bAddAlpha;
     191                 :     int                 bHideNoData;
     192                 :     char               *pszSrcNoData;
     193                 :     char               *pszVRTNoData;
     194                 : 
     195                 :     /* Internal variables */
     196                 :     char               *pszProjectionRef;
     197                 :     int                 nBands;
     198                 :     BandProperty       *pasBandProperties;
     199                 :     int                 bFirst;
     200                 :     int                 bHasGeoTransform;
     201                 :     int                 nRasterXSize;
     202                 :     int                 nRasterYSize;
     203                 :     DatasetProperty    *pasDatasetProperties;
     204                 :     int                 bUserExtent;
     205                 :     int                 bAllowSrcNoData;
     206                 :     double             *padfSrcNoData;
     207                 :     int                 nSrcNoDataCount;
     208                 :     int                 bAllowVRTNoData;
     209                 :     double             *padfVRTNoData;
     210                 :     int                 nVRTNoDataCount;
     211                 :     int                 bHasRunBuild;
     212                 : 
     213                 :     int         AnalyseRaster(GDALDatasetH hDS, const char* dsFileName,
     214                 :                               DatasetProperty* psDatasetProperties);
     215                 : 
     216                 :     void        CreateVRTSeparate(VRTDatasetH hVRTDS);
     217                 :     void        CreateVRTNonSeparate(VRTDatasetH hVRTDS);
     218                 : 
     219                 :     public:
     220                 :                 VRTBuilder(const char* pszOutputFilename,
     221                 :                            int nInputFiles, const char* const * ppszInputFilenames,
     222                 :                            ResolutionStrategy resolutionStrategy,
     223                 :                            double we_res, double ns_res,
     224                 :                            double minX, double minY, double maxX, double maxY,
     225                 :                            int bSeparate, int bAllowProjectionDifference,
     226                 :                            int bAddAlpha, int bHideNoData,
     227                 :                            const char* pszSrcNoData, const char* pszVRTNoData);
     228                 : 
     229                 :                ~VRTBuilder();
     230                 : 
     231                 :         int     Build(GDALProgressFunc pfnProgress, void * pProgressData);
     232                 : };
     233                 : 
     234                 : 
     235                 : /************************************************************************/
     236                 : /*                          VRTBuilder()                                */
     237                 : /************************************************************************/
     238                 : 
     239                 : VRTBuilder::VRTBuilder(const char* pszOutputFilename,
     240                 :                        int nInputFiles, const char* const * ppszInputFilenames,
     241                 :                        ResolutionStrategy resolutionStrategy,
     242                 :                        double we_res, double ns_res,
     243                 :                        double minX, double minY, double maxX, double maxY,
     244                 :                        int bSeparate, int bAllowProjectionDifference,
     245                 :                        int bAddAlpha, int bHideNoData,
     246              13 :                        const char* pszSrcNoData, const char* pszVRTNoData)
     247                 : {
     248              13 :     this->pszOutputFilename = CPLStrdup(pszOutputFilename);
     249              13 :     this->nInputFiles = nInputFiles;
     250                 : 
     251              13 :     this->ppszInputFilenames = (char**) CPLMalloc(nInputFiles * sizeof(char*));
     252                 :     int i;
     253              55 :     for(i=0;i<nInputFiles;i++)
     254                 :     {
     255              42 :         this->ppszInputFilenames[i] = CPLStrdup(ppszInputFilenames[i]);
     256                 :     }
     257                 : 
     258              13 :     this->resolutionStrategy = resolutionStrategy;
     259              13 :     this->we_res = we_res;
     260              13 :     this->ns_res = ns_res;
     261              13 :     this->minX = minX;
     262              13 :     this->minY = minY;
     263              13 :     this->maxX = maxX;
     264              13 :     this->maxY = maxY;
     265              13 :     this->bSeparate = bSeparate;
     266              13 :     this->bAllowProjectionDifference = bAllowProjectionDifference;
     267              13 :     this->bAddAlpha = bAddAlpha;
     268              13 :     this->bHideNoData = bHideNoData;
     269              13 :     this->pszSrcNoData = (pszSrcNoData) ? CPLStrdup(pszSrcNoData) : NULL;
     270              13 :     this->pszVRTNoData = (pszVRTNoData) ? CPLStrdup(pszVRTNoData) : NULL;
     271                 : 
     272              13 :     bUserExtent = FALSE;
     273              13 :     pszProjectionRef = NULL;
     274              13 :     nBands = 0;
     275              13 :     pasBandProperties = NULL;
     276              13 :     bFirst = TRUE;
     277              13 :     bHasGeoTransform = FALSE;
     278              13 :     nRasterXSize = 0;
     279              13 :     nRasterYSize = 0;
     280              13 :     pasDatasetProperties = NULL;
     281              13 :     bAllowSrcNoData = TRUE;
     282              13 :     padfSrcNoData = NULL;
     283              13 :     nSrcNoDataCount = 0;
     284              13 :     bAllowVRTNoData = TRUE;
     285              13 :     padfVRTNoData = NULL;
     286              13 :     nVRTNoDataCount = 0;
     287              13 :     bHasRunBuild = FALSE;
     288              13 : }
     289                 : 
     290                 : /************************************************************************/
     291                 : /*                         ~VRTBuilder()                                */
     292                 : /************************************************************************/
     293                 : 
     294              13 : VRTBuilder::~VRTBuilder()
     295                 : {
     296              13 :     CPLFree(pszOutputFilename);
     297              13 :     CPLFree(pszSrcNoData);
     298              13 :     CPLFree(pszVRTNoData);
     299                 : 
     300                 :     int i;
     301              55 :     for(i=0;i<nInputFiles;i++)
     302                 :     {
     303              42 :         CPLFree(ppszInputFilenames[i]);
     304                 :     }
     305              13 :     CPLFree(ppszInputFilenames);
     306                 : 
     307              13 :     if (pasDatasetProperties != NULL)
     308                 :     {
     309              55 :         for(i=0;i<nInputFiles;i++)
     310                 :         {
     311              42 :             CPLFree(pasDatasetProperties[i].padfNoDataValues);
     312              42 :             CPLFree(pasDatasetProperties[i].panHasNoData);
     313                 :         }
     314                 :     }
     315              13 :     CPLFree(pasDatasetProperties);
     316                 : 
     317              13 :     if (!bSeparate && pasBandProperties != NULL)
     318                 :     {
     319                 :         int j;
     320              24 :         for(j=0;j<nBands;j++)
     321                 :         {
     322              13 :             GDALDestroyColorTable(pasBandProperties[j].colorTable);
     323                 :         }
     324                 :     }
     325              13 :     CPLFree(pasBandProperties);
     326                 : 
     327              13 :     CPLFree(pszProjectionRef);
     328              13 :     CPLFree(padfSrcNoData);
     329              13 :     CPLFree(padfVRTNoData);
     330              13 : }
     331                 : 
     332                 : /************************************************************************/
     333                 : /*                           AnalyseRaster()                            */
     334                 : /************************************************************************/
     335                 : 
     336                 : int VRTBuilder::AnalyseRaster( GDALDatasetH hDS, const char* dsFileName,
     337              42 :                                   DatasetProperty* psDatasetProperties)
     338                 : {
     339              42 :     char** papszMetadata = GDALGetMetadata( hDS, "SUBDATASETS" );
     340              42 :     if( CSLCount(papszMetadata) > 0 && GDALGetRasterCount(hDS) == 0 )
     341                 :     {
     342                 :         pasDatasetProperties =
     343                 :             (DatasetProperty*) CPLRealloc(pasDatasetProperties,
     344               0 :                             (nInputFiles+CSLCount(papszMetadata))*sizeof(DatasetProperty));
     345                 : 
     346                 :         ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
     347               0 :                                 sizeof(char*) * (nInputFiles+CSLCount(papszMetadata)));
     348               0 :         int count = 1;
     349                 :         char subdatasetNameKey[256];
     350               0 :         sprintf(subdatasetNameKey, "SUBDATASET_%d_NAME", count);
     351               0 :         while(*papszMetadata != NULL)
     352                 :         {
     353               0 :             if (EQUALN(*papszMetadata, subdatasetNameKey, strlen(subdatasetNameKey)))
     354                 :             {
     355               0 :                 memset(&pasDatasetProperties[nInputFiles], 0, sizeof(DatasetProperty));
     356                 :                 ppszInputFilenames[nInputFiles++] =
     357               0 :                         CPLStrdup(*papszMetadata+strlen(subdatasetNameKey)+1);
     358               0 :                 count++;
     359               0 :                 sprintf(subdatasetNameKey, "SUBDATASET_%d_NAME", count);
     360                 :             }
     361               0 :             papszMetadata++;
     362                 :         }
     363               0 :         return FALSE;
     364                 :     }
     365                 : 
     366              42 :     const char* proj = GDALGetProjectionRef(hDS);
     367              42 :     double* padfGeoTransform = psDatasetProperties->adfGeoTransform;
     368              42 :     int bGotGeoTransform = GDALGetGeoTransform(hDS, padfGeoTransform) == CE_None;
     369              42 :     if (bSeparate)
     370                 :     {
     371               6 :         if (bFirst)
     372                 :         {
     373               2 :             bHasGeoTransform = bGotGeoTransform;
     374               2 :             if (!bHasGeoTransform)
     375                 :             {
     376               1 :                 if (bUserExtent)
     377                 :                 {
     378                 :                     CPLError(CE_Warning, CPLE_NotSupported,
     379               0 :                         "User extent ignored by gdalbuildvrt -separate with ungeoreferenced images.");
     380                 :                 }
     381               1 :                 if (resolutionStrategy == USER_RESOLUTION)
     382                 :                 {
     383                 :                     CPLError(CE_Warning, CPLE_NotSupported,
     384               0 :                         "User resolution ignored by gdalbuildvrt -separate with ungeoreferenced images.");
     385                 :                 }
     386                 :             }
     387                 :         }
     388               4 :         else if (bHasGeoTransform != bGotGeoTransform)
     389                 :         {
     390                 :             CPLError(CE_Warning, CPLE_NotSupported,
     391                 :                     "gdalbuildvrt -separate cannot stack ungeoreferenced and georeferenced images. Skipping %s",
     392               0 :                     dsFileName);
     393               0 :             return FALSE;
     394                 :         }
     395               4 :         else if (!bHasGeoTransform &&
     396                 :                     (nRasterXSize != GDALGetRasterXSize(hDS) ||
     397                 :                     nRasterYSize != GDALGetRasterYSize(hDS)))
     398                 :         {
     399                 :             CPLError(CE_Warning, CPLE_NotSupported,
     400                 :                     "gdalbuildvrt -separate cannot stack ungeoreferenced images that have not the same dimensions. Skipping %s",
     401               0 :                     dsFileName);
     402               0 :             return FALSE;
     403                 :         }
     404                 :     }
     405                 :     else
     406                 :     {
     407              36 :         if (!bGotGeoTransform)
     408                 :         {
     409                 :             CPLError(CE_Warning, CPLE_NotSupported,
     410                 :                     "gdalbuildvrt does not support ungeoreferenced image. Skipping %s",
     411               0 :                     dsFileName);
     412               0 :             return FALSE;
     413                 :         }
     414              36 :         bHasGeoTransform = TRUE;
     415                 :     }
     416                 : 
     417              42 :     if (bGotGeoTransform)
     418                 :     {
     419              40 :         if (padfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] != 0 ||
     420                 :             padfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] != 0)
     421                 :         {
     422                 :             CPLError(CE_Warning, CPLE_NotSupported,
     423                 :                     "gdalbuildvrt does not support rotated geo transforms. Skipping %s",
     424               0 :                     dsFileName);
     425               0 :             return FALSE;
     426                 :         }
     427              40 :         if (padfGeoTransform[GEOTRSFRM_NS_RES] >= 0)
     428                 :         {
     429                 :             CPLError(CE_Warning, CPLE_NotSupported,
     430                 :                     "gdalbuildvrt does not support positive NS resolution. Skipping %s",
     431               0 :                     dsFileName);
     432               0 :             return FALSE;
     433                 :         }
     434                 :     }
     435                 : 
     436              42 :     psDatasetProperties->nRasterXSize = GDALGetRasterXSize(hDS);
     437              42 :     psDatasetProperties->nRasterYSize = GDALGetRasterYSize(hDS);
     438              42 :     if (bFirst && bSeparate && !bGotGeoTransform)
     439                 :     {
     440               1 :         nRasterXSize = GDALGetRasterXSize(hDS);
     441               1 :         nRasterYSize = GDALGetRasterYSize(hDS);
     442                 :     }
     443                 : 
     444              42 :     double ds_minX = padfGeoTransform[GEOTRSFRM_TOPLEFT_X];
     445              42 :     double ds_maxY = padfGeoTransform[GEOTRSFRM_TOPLEFT_Y];
     446                 :     double ds_maxX = ds_minX +
     447                 :                 GDALGetRasterXSize(hDS) *
     448              42 :                 padfGeoTransform[GEOTRSFRM_WE_RES];
     449                 :     double ds_minY = ds_maxY +
     450                 :                 GDALGetRasterYSize(hDS) *
     451              42 :                 padfGeoTransform[GEOTRSFRM_NS_RES];
     452                 : 
     453                 :     GDALGetBlockSize(GDALGetRasterBand( hDS, 1 ),
     454                 :                         &psDatasetProperties->nBlockXSize,
     455              42 :                         &psDatasetProperties->nBlockYSize);
     456                 : 
     457              42 :     int _nBands = GDALGetRasterCount(hDS);
     458              42 :     if (_nBands == 0)
     459                 :     {
     460                 :         CPLError(CE_Warning, CPLE_AppDefined,
     461               0 :                     "Skipping %s as it has no bands", dsFileName);
     462               0 :         return FALSE;
     463                 :     }
     464              42 :     else if (_nBands > 1 && bSeparate)
     465                 :     {
     466                 :         CPLError(CE_Warning, CPLE_AppDefined, "%s has %d bands. Only the first one will "
     467                 :                     "be taken into account in the -separate case",
     468               0 :                     dsFileName, _nBands);
     469               0 :         _nBands = 1;
     470                 :     }
     471                 : 
     472                 :     /* For the -separate case */
     473              42 :     psDatasetProperties->firstBandType = GDALGetRasterDataType(GDALGetRasterBand(hDS, 1));
     474                 : 
     475              42 :     psDatasetProperties->padfNoDataValues = (double*)CPLCalloc(sizeof(double), _nBands);
     476              42 :     psDatasetProperties->panHasNoData = (int*)CPLCalloc(sizeof(int), _nBands);
     477                 :     int j;
     478              89 :     for(j=0;j<_nBands;j++)
     479                 :     {
     480              47 :         if (nSrcNoDataCount > 0)
     481                 :         {
     482               2 :             psDatasetProperties->panHasNoData[j] = TRUE;
     483               2 :             if (j < nSrcNoDataCount)
     484               2 :                 psDatasetProperties->padfNoDataValues[j] = padfSrcNoData[j];
     485                 :             else
     486               0 :                 psDatasetProperties->padfNoDataValues[j] = padfSrcNoData[nSrcNoDataCount - 1];
     487                 :         }
     488                 :         else
     489                 :         {
     490                 :             psDatasetProperties->padfNoDataValues[j]  =
     491                 :                 GDALGetRasterNoDataValue(GDALGetRasterBand(hDS, j+1),
     492              45 :                                         &psDatasetProperties->panHasNoData[j]);
     493                 :         }
     494                 :     }
     495                 : 
     496              42 :     if (bFirst)
     497                 :     {
     498              13 :         if (proj)
     499              13 :             pszProjectionRef = CPLStrdup(proj);
     500              13 :         if (!bUserExtent)
     501                 :         {
     502              11 :             minX = ds_minX;
     503              11 :             minY = ds_minY;
     504              11 :             maxX = ds_maxX;
     505              11 :             maxY = ds_maxY;
     506                 :         }
     507              13 :         nBands = _nBands;
     508                 : 
     509              13 :         if (!bSeparate)
     510                 :         {
     511              11 :             pasBandProperties = (BandProperty*)CPLMalloc(nBands*sizeof(BandProperty));
     512              24 :             for(j=0;j<nBands;j++)
     513                 :             {
     514              13 :                 GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
     515                 :                 pasBandProperties[j].colorInterpretation =
     516              13 :                         GDALGetRasterColorInterpretation(hRasterBand);
     517              13 :                 pasBandProperties[j].dataType = GDALGetRasterDataType(hRasterBand);
     518              13 :                 if (pasBandProperties[j].colorInterpretation == GCI_PaletteIndex)
     519                 :                 {
     520                 :                     pasBandProperties[j].colorTable =
     521               1 :                             GDALGetRasterColorTable( hRasterBand );
     522               1 :                     if (pasBandProperties[j].colorTable)
     523                 :                     {
     524                 :                         pasBandProperties[j].colorTable =
     525               1 :                                 GDALCloneColorTable(pasBandProperties[j].colorTable);
     526                 :                     }
     527                 :                 }
     528                 :                 else
     529              12 :                     pasBandProperties[j].colorTable = 0;
     530                 : 
     531              13 :                 if (nVRTNoDataCount > 0)
     532                 :                 {
     533               1 :                     pasBandProperties[j].bHasNoData = TRUE;
     534               1 :                     if (j < nVRTNoDataCount)
     535               1 :                         pasBandProperties[j].noDataValue = padfVRTNoData[j];
     536                 :                     else
     537               0 :                         pasBandProperties[j].noDataValue = padfVRTNoData[nVRTNoDataCount - 1];
     538                 :                 }
     539                 :                 else
     540                 :                 {
     541                 :                     pasBandProperties[j].noDataValue =
     542              12 :                             GDALGetRasterNoDataValue(hRasterBand, &pasBandProperties[j].bHasNoData);
     543                 :                 }
     544                 :             }
     545                 :         }
     546                 :     }
     547                 :     else
     548                 :     {
     549              29 :         if ((proj != NULL && pszProjectionRef == NULL) ||
     550                 :             (proj == NULL && pszProjectionRef != NULL) ||
     551                 :             (proj != NULL && pszProjectionRef != NULL && EQUAL(proj, pszProjectionRef) == FALSE))
     552                 :         {
     553               1 :             if (!bAllowProjectionDifference)
     554                 :             {
     555                 :                 CPLError(CE_Warning, CPLE_NotSupported,
     556                 :                             "gdalbuildvrt does not support heterogenous projection. Skipping %s",
     557               1 :                             dsFileName);
     558               1 :                 return FALSE;
     559                 :             }
     560                 :         }
     561              28 :         if (!bSeparate)
     562                 :         {
     563              24 :             if (nBands != _nBands)
     564                 :             {
     565                 :                 CPLError(CE_Warning, CPLE_NotSupported,
     566                 :                             "gdalbuildvrt does not support heterogenous band numbers. Skipping %s",
     567               1 :                         dsFileName);
     568               1 :                 return FALSE;
     569                 :             }
     570              48 :             for(j=0;j<nBands;j++)
     571                 :             {
     572              25 :                 GDALRasterBandH hRasterBand = GDALGetRasterBand( hDS, j+1 );
     573              25 :                 if (pasBandProperties[j].colorInterpretation != GDALGetRasterColorInterpretation(hRasterBand) ||
     574                 :                     pasBandProperties[j].dataType != GDALGetRasterDataType(hRasterBand))
     575                 :                 {
     576                 :                     CPLError(CE_Warning, CPLE_NotSupported,
     577                 :                                 "gdalbuildvrt does not support heterogenous band characteristics. Skipping %s",
     578               0 :                                 dsFileName);
     579               0 :                     return FALSE;
     580                 :                 }
     581              25 :                 if (pasBandProperties[j].colorTable)
     582                 :                 {
     583               1 :                     GDALColorTableH colorTable = GDALGetRasterColorTable( hRasterBand );
     584               1 :                     int nRefColorEntryCount = GDALGetColorEntryCount(pasBandProperties[j].colorTable);
     585                 :                     int i;
     586               1 :                     if (colorTable == NULL ||
     587                 :                         GDALGetColorEntryCount(colorTable) != nRefColorEntryCount)
     588                 :                     {
     589                 :                         CPLError(CE_Warning, CPLE_NotSupported,
     590                 :                                     "gdalbuildvrt does not support rasters with different color tables (different number of color table entries). Skipping %s",
     591               0 :                                 dsFileName);
     592               0 :                         return FALSE;
     593                 :                     }
     594                 : 
     595                 :                     /* Check that the palette are the same too */
     596                 :                     /* We just warn and still process the file. It is not a technical no-go, but the user */
     597                 :                     /* should check that the end result is OK for him. */
     598               3 :                     for(i=0;i<nRefColorEntryCount;i++)
     599                 :                     {
     600               2 :                         const GDALColorEntry* psEntry = GDALGetColorEntry(colorTable, i);
     601               2 :                         const GDALColorEntry* psEntryRef = GDALGetColorEntry(pasBandProperties[j].colorTable, i);
     602               2 :                         if (psEntry->c1 != psEntryRef->c1 || psEntry->c2 != psEntryRef->c2 ||
     603                 :                             psEntry->c3 != psEntryRef->c3 || psEntry->c4 != psEntryRef->c4)
     604                 :                         {
     605                 :                             static int bFirstWarningPCT = TRUE;
     606               0 :                             if (bFirstWarningPCT)
     607                 :                                 CPLError(CE_Warning, CPLE_NotSupported,
     608                 :                                         "%s has different values than the first raster for some entries in the color table.\n"
     609                 :                                         "The end result might produce weird colors.\n"
     610                 :                                         "You're advised to preprocess your rasters with other tools, such as pct2rgb.py or gdal_translate -expand RGB\n"
     611               0 :                                         "to operate gdalbuildvrt on RGB rasters instead", dsFileName);
     612                 :                             else
     613                 :                                 CPLError(CE_Warning, CPLE_NotSupported,
     614                 :                                             "%s has different values than the first raster for some entries in the color table.",
     615               0 :                                             dsFileName);
     616               0 :                             bFirstWarningPCT = FALSE;
     617               0 :                             break;
     618                 :                         }
     619                 :                     }
     620                 :                 }
     621                 :             }
     622                 : 
     623                 :         }
     624              27 :         if (!bUserExtent)
     625                 :         {
     626              24 :             if (ds_minX < minX) minX = ds_minX;
     627              24 :             if (ds_minY < minY) minY = ds_minY;
     628              24 :             if (ds_maxX > maxX) maxX = ds_maxX;
     629              24 :             if (ds_maxY > maxY) maxY = ds_maxY;
     630                 :         }
     631                 :     }
     632                 : 
     633              40 :     if (resolutionStrategy == AVERAGE_RESOLUTION)
     634                 :     {
     635              35 :         we_res += padfGeoTransform[GEOTRSFRM_WE_RES];
     636              35 :         ns_res += padfGeoTransform[GEOTRSFRM_NS_RES];
     637                 :     }
     638               5 :     else if (resolutionStrategy != USER_RESOLUTION)
     639                 :     {
     640               0 :         if (bFirst)
     641                 :         {
     642               0 :             we_res = padfGeoTransform[GEOTRSFRM_WE_RES];
     643               0 :             ns_res = padfGeoTransform[GEOTRSFRM_NS_RES];
     644                 :         }
     645               0 :         else if (resolutionStrategy == HIGHEST_RESOLUTION)
     646                 :         {
     647               0 :             we_res = MIN(we_res, padfGeoTransform[GEOTRSFRM_WE_RES]);
     648                 :             /* Yes : as ns_res is negative, the highest resolution is the max value */
     649               0 :             ns_res = MAX(ns_res, padfGeoTransform[GEOTRSFRM_NS_RES]);
     650                 :         }
     651                 :         else
     652                 :         {
     653               0 :             we_res = MAX(we_res, padfGeoTransform[GEOTRSFRM_WE_RES]);
     654                 :             /* Yes : as ns_res is negative, the lowest resolution is the min value */
     655               0 :             ns_res = MIN(ns_res, padfGeoTransform[GEOTRSFRM_NS_RES]);
     656                 :         }
     657                 :     }
     658                 : 
     659              40 :     return TRUE;
     660                 : }
     661                 : 
     662                 : /************************************************************************/
     663                 : /*                         CreateVRTSeparate()                          */
     664                 : /************************************************************************/
     665                 : 
     666               2 : void VRTBuilder::CreateVRTSeparate(VRTDatasetH hVRTDS)
     667                 : {
     668                 :     int i;
     669               2 :     int iBand = 1;
     670               8 :     for(i=0;i<nInputFiles;i++)
     671                 :     {
     672               6 :         DatasetProperty* psDatasetProperties = &pasDatasetProperties[i];
     673                 : 
     674               6 :         if (psDatasetProperties->isFileOK == FALSE)
     675               0 :             continue;
     676                 : 
     677                 :         int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     678                 :             nDstXOff, nDstYOff, nDstXSize, nDstYSize;
     679               6 :         if (bHasGeoTransform)
     680                 :         {
     681               4 :             if ( ! GetSrcDstWin(psDatasetProperties,
     682                 :                         we_res, ns_res, minX, minY, maxX, maxY,
     683                 :                         &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
     684                 :                         &nDstXOff, &nDstYOff, &nDstXSize, &nDstYSize) )
     685               0 :                 continue;
     686                 :         }
     687                 :         else
     688                 :         {
     689               2 :             nSrcXOff = nSrcYOff = nDstXOff = nDstYOff = 0;
     690               2 :             nSrcXSize = nDstXSize = nRasterXSize;
     691               2 :             nSrcYSize = nDstYSize = nRasterYSize;
     692                 :         }
     693                 : 
     694               6 :         const char* dsFileName = ppszInputFilenames[i];
     695                 : 
     696               6 :         GDALAddBand(hVRTDS, psDatasetProperties->firstBandType, NULL);
     697                 : 
     698                 :         GDALProxyPoolDatasetH hProxyDS =
     699                 :             GDALProxyPoolDatasetCreate(dsFileName,
     700                 :                                         psDatasetProperties->nRasterXSize,
     701                 :                                         psDatasetProperties->nRasterYSize,
     702                 :                                         GA_ReadOnly, TRUE, pszProjectionRef,
     703               6 :                                         psDatasetProperties->adfGeoTransform);
     704                 :         GDALProxyPoolDatasetAddSrcBandDescription(hProxyDS,
     705                 :                                             psDatasetProperties->firstBandType,
     706                 :                                             psDatasetProperties->nBlockXSize,
     707               6 :                                             psDatasetProperties->nBlockYSize);
     708                 : 
     709                 :         VRTSourcedRasterBandH hVRTBand =
     710               6 :                 (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, iBand);
     711                 : 
     712               6 :         if (bHideNoData)
     713               0 :             GDALSetMetadataItem(hVRTBand,"HideNoDataValue","1",NULL);
     714                 : 
     715               6 :         if (bAllowSrcNoData && psDatasetProperties->panHasNoData[0])
     716                 :         {
     717               0 :             GDALSetRasterNoDataValue(hVRTBand, psDatasetProperties->padfNoDataValues[0]);
     718                 :             VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
     719                 :                             nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     720                 :                             nDstXOff, nDstYOff, nDstXSize, nDstYSize,
     721               0 :                             0, 1, psDatasetProperties->padfNoDataValues[0]);
     722                 :         }
     723                 :         else
     724                 :             /* Place the raster band at the right position in the VRT */
     725                 :             VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
     726                 :                             nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     727                 :                             nDstXOff, nDstYOff, nDstXSize, nDstYSize,
     728               6 :                             "near", VRT_NODATA_UNSET);
     729                 : 
     730               6 :         GDALDereferenceDataset(hProxyDS);
     731                 : 
     732               6 :         iBand ++;
     733                 :     }
     734               2 : }
     735                 : 
     736                 : /************************************************************************/
     737                 : /*                       CreateVRTNonSeparate()                         */
     738                 : /************************************************************************/
     739                 : 
     740              11 : void VRTBuilder::CreateVRTNonSeparate(VRTDatasetH hVRTDS)
     741                 : {
     742                 :     int i, j;
     743                 : 
     744              24 :     for(j=0;j<nBands;j++)
     745                 :     {
     746                 :         GDALRasterBandH hBand;
     747              13 :         GDALAddBand(hVRTDS, pasBandProperties[j].dataType, NULL);
     748              13 :         hBand = GDALGetRasterBand(hVRTDS, j+1);
     749              13 :         GDALSetRasterColorInterpretation(hBand, pasBandProperties[j].colorInterpretation);
     750              13 :         if (pasBandProperties[j].colorInterpretation == GCI_PaletteIndex)
     751                 :         {
     752               1 :             GDALSetRasterColorTable(hBand, pasBandProperties[j].colorTable);
     753                 :         }
     754              13 :         if (bAllowVRTNoData && pasBandProperties[j].bHasNoData)
     755               4 :             GDALSetRasterNoDataValue(hBand, pasBandProperties[j].noDataValue);
     756              13 :         if ( bHideNoData )
     757               0 :             GDALSetMetadataItem(hBand,"HideNoDataValue","1",NULL);
     758                 :     }
     759                 : 
     760              11 :     if (bAddAlpha)
     761                 :     {
     762                 :         GDALRasterBandH hBand;
     763               0 :         GDALAddBand(hVRTDS, GDT_Byte, NULL);
     764               0 :         hBand = GDALGetRasterBand(hVRTDS, nBands + 1);
     765               0 :         GDALSetRasterColorInterpretation(hBand, GCI_AlphaBand);
     766                 :     }
     767                 : 
     768              47 :     for(i=0;i<nInputFiles;i++)
     769                 :     {
     770              36 :         DatasetProperty* psDatasetProperties = &pasDatasetProperties[i];
     771                 : 
     772              36 :         if (psDatasetProperties->isFileOK == FALSE)
     773               2 :             continue;
     774                 : 
     775                 :         int nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     776                 :             nDstXOff, nDstYOff, nDstXSize, nDstYSize;
     777              34 :         if ( ! GetSrcDstWin(psDatasetProperties,
     778                 :                         we_res, ns_res, minX, minY, maxX, maxY,
     779                 :                         &nSrcXOff, &nSrcYOff, &nSrcXSize, &nSrcYSize,
     780                 :                         &nDstXOff, &nDstYOff, &nDstXSize, &nDstYSize) )
     781               0 :             continue;
     782                 : 
     783              34 :         const char* dsFileName = ppszInputFilenames[i];
     784                 : 
     785                 :         GDALProxyPoolDatasetH hProxyDS =
     786                 :             GDALProxyPoolDatasetCreate(dsFileName,
     787                 :                                         psDatasetProperties->nRasterXSize,
     788                 :                                         psDatasetProperties->nRasterYSize,
     789                 :                                         GA_ReadOnly, TRUE, pszProjectionRef,
     790              34 :                                         psDatasetProperties->adfGeoTransform);
     791                 : 
     792              72 :         for(j=0;j<nBands;j++)
     793                 :         {
     794                 :             GDALProxyPoolDatasetAddSrcBandDescription(hProxyDS,
     795                 :                                             pasBandProperties[j].dataType,
     796                 :                                             psDatasetProperties->nBlockXSize,
     797              38 :                                             psDatasetProperties->nBlockYSize);
     798                 :         }
     799                 : 
     800              72 :         for(j=0;j<nBands;j++)
     801                 :         {
     802                 :             VRTSourcedRasterBandH hVRTBand =
     803              38 :                     (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, j + 1);
     804                 : 
     805                 :             /* Place the raster band at the right position in the VRT */
     806              46 :             if (bAllowSrcNoData && psDatasetProperties->panHasNoData[j])
     807                 :                 VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, j + 1),
     808                 :                                 nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     809                 :                                 nDstXOff, nDstYOff, nDstXSize, nDstYSize,
     810               8 :                                 0, 1, psDatasetProperties->padfNoDataValues[j]);
     811                 :             else
     812                 :                 VRTAddSimpleSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, j + 1),
     813                 :                                 nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     814                 :                                 nDstXOff, nDstYOff, nDstXSize, nDstYSize,
     815              30 :                                 "near", VRT_NODATA_UNSET);
     816                 :         }
     817                 : 
     818              34 :         if (bAddAlpha)
     819                 :         {
     820                 :             VRTSourcedRasterBandH hVRTBand =
     821               0 :                     (VRTSourcedRasterBandH)GDALGetRasterBand(hVRTDS, nBands + 1);
     822                 :             /* Little trick : we use an offset of 255 and a scaling of 0, so that in areas covered */
     823                 :             /* by the source, the value of the alpha band will be 255, otherwise it will be 0 */
     824                 :             VRTAddComplexSource(hVRTBand, GDALGetRasterBand((GDALDatasetH)hProxyDS, 1),
     825                 :                                 nSrcXOff, nSrcYOff, nSrcXSize, nSrcYSize,
     826                 :                                 nDstXOff, nDstYOff, nDstXSize, nDstYSize,
     827               0 :                                 255, 0, VRT_NODATA_UNSET);
     828                 :         }
     829                 : 
     830              34 :         GDALDereferenceDataset(hProxyDS);
     831                 :     }
     832              11 : }
     833                 : 
     834                 : /************************************************************************/
     835                 : /*                             Build()                                  */
     836                 : /************************************************************************/
     837                 : 
     838              13 : int VRTBuilder::Build(GDALProgressFunc pfnProgress, void * pProgressData)
     839                 : {
     840                 :     int i;
     841                 : 
     842              13 :     if (bHasRunBuild)
     843               0 :         return CE_Failure;
     844              13 :     bHasRunBuild = TRUE;
     845                 : 
     846              13 :     if( pfnProgress == NULL )
     847               0 :         pfnProgress = GDALDummyProgress;
     848                 : 
     849              13 :     bUserExtent = (minX != 0 || minY != 0 || maxX != 0 || maxY != 0);
     850              13 :     if (bUserExtent)
     851                 :     {
     852               2 :         if (minX >= maxX || minY >= maxY )
     853                 :         {
     854               0 :             CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user extent");
     855               0 :             return CE_Failure;
     856                 :         }
     857                 :     }
     858                 : 
     859              13 :     if (resolutionStrategy == USER_RESOLUTION)
     860                 :     {
     861               2 :         if (we_res <= 0 || ns_res <= 0)
     862                 :         {
     863               0 :             CPLError(CE_Failure, CPLE_IllegalArg, "Invalid user resolution");
     864               0 :             return CE_Failure;
     865                 :         }
     866                 : 
     867                 :         /* We work with negative north-south resolution in all the following code */
     868               2 :         ns_res = -ns_res;
     869                 :     }
     870                 :     else
     871                 :     {
     872              11 :         we_res = ns_res = 0;
     873                 :     }
     874                 : 
     875                 :     pasDatasetProperties =
     876              13 :             (DatasetProperty*) CPLCalloc(nInputFiles, sizeof(DatasetProperty));
     877                 : 
     878              13 :     if (pszSrcNoData != NULL)
     879                 :     {
     880               1 :         if (EQUAL(pszSrcNoData, "none"))
     881                 :         {
     882               0 :             bAllowSrcNoData = FALSE;
     883                 :         }
     884                 :         else
     885                 :         {
     886               1 :             char **papszTokens = CSLTokenizeString( pszSrcNoData );
     887               1 :             nSrcNoDataCount = CSLCount(papszTokens);
     888               1 :             padfSrcNoData = (double *) CPLMalloc(sizeof(double) * nSrcNoDataCount);
     889               2 :             for(i=0;i<nSrcNoDataCount;i++)
     890               1 :                 padfSrcNoData[i] = CPLAtofM(papszTokens[i]);
     891               1 :             CSLDestroy(papszTokens);
     892                 :         }
     893                 :     }
     894                 : 
     895              13 :     if (pszVRTNoData != NULL)
     896                 :     {
     897               1 :         if (EQUAL(pszVRTNoData, "none"))
     898                 :         {
     899               0 :             bAllowVRTNoData = FALSE;
     900                 :         }
     901                 :         else
     902                 :         {
     903               1 :             char **papszTokens = CSLTokenizeString( pszVRTNoData );
     904               1 :             nVRTNoDataCount = CSLCount(papszTokens);
     905               1 :             padfVRTNoData = (double *) CPLMalloc(sizeof(double) * nVRTNoDataCount);
     906               2 :             for(i=0;i<nVRTNoDataCount;i++)
     907               1 :                 padfVRTNoData[i] = CPLAtofM(papszTokens[i]);
     908               1 :             CSLDestroy(papszTokens);
     909                 :         }
     910                 :     }
     911                 : 
     912              13 :     int nCountValid = 0;
     913              55 :     for(i=0;i<nInputFiles;i++)
     914                 :     {
     915              42 :         const char* dsFileName = ppszInputFilenames[i];
     916                 : 
     917              42 :         if (!pfnProgress( 1.0 * (i+1) / nInputFiles, NULL, pProgressData))
     918                 :         {
     919               0 :             return CE_Failure;
     920                 :         }
     921                 : 
     922              42 :         GDALDatasetH hDS = GDALOpen(ppszInputFilenames[i], GA_ReadOnly );
     923              42 :         pasDatasetProperties[i].isFileOK = FALSE;
     924                 : 
     925              42 :         if (hDS)
     926                 :         {
     927              42 :             if (AnalyseRaster( hDS, dsFileName, &pasDatasetProperties[i] ))
     928                 :             {
     929              40 :                 pasDatasetProperties[i].isFileOK = TRUE;
     930              40 :                 nCountValid ++;
     931              40 :                 bFirst = FALSE;
     932                 :             }
     933              42 :             GDALClose(hDS);
     934                 :         }
     935                 :         else
     936                 :         {
     937                 :             CPLError(CE_Warning, CPLE_AppDefined, 
     938               0 :                      "Can't open %s. Skipping it", dsFileName);
     939                 :         }
     940                 :     }
     941                 : 
     942              13 :     if (nCountValid == 0)
     943               0 :         return CE_None;
     944                 : 
     945              13 :     if (bHasGeoTransform)
     946                 :     {
     947              12 :         if (resolutionStrategy == AVERAGE_RESOLUTION)
     948                 :         {
     949              10 :             we_res /= nCountValid;
     950              10 :             ns_res /= nCountValid;
     951                 :         }
     952                 : 
     953              12 :         nRasterXSize = (int)(0.5 + (maxX - minX) / we_res);
     954              12 :         nRasterYSize = (int)(0.5 + (maxY - minY) / -ns_res);
     955                 :     }
     956                 : 
     957              13 :     if (nRasterXSize == 0 || nRasterYSize == 0)
     958                 :     {
     959                 :         CPLError(CE_Failure, CPLE_AppDefined, 
     960               0 :                   "Computed VRT dimension is invalid. You've probably specified unappropriate resolution.");
     961               0 :         return CE_Failure;
     962                 :     }
     963                 : 
     964              13 :     VRTDatasetH hVRTDS = VRTCreate(nRasterXSize, nRasterYSize);
     965              13 :     GDALSetDescription(hVRTDS, pszOutputFilename);
     966                 : 
     967              13 :     if (pszProjectionRef)
     968                 :     {
     969              13 :         GDALSetProjection(hVRTDS, pszProjectionRef);
     970                 :     }
     971                 : 
     972              13 :     if (bHasGeoTransform)
     973                 :     {
     974                 :         double adfGeoTransform[6];
     975              12 :         adfGeoTransform[GEOTRSFRM_TOPLEFT_X] = minX;
     976              12 :         adfGeoTransform[GEOTRSFRM_WE_RES] = we_res;
     977              12 :         adfGeoTransform[GEOTRSFRM_ROTATION_PARAM1] = 0;
     978              12 :         adfGeoTransform[GEOTRSFRM_TOPLEFT_Y] = maxY;
     979              12 :         adfGeoTransform[GEOTRSFRM_ROTATION_PARAM2] = 0;
     980              12 :         adfGeoTransform[GEOTRSFRM_NS_RES] = ns_res;
     981              12 :         GDALSetGeoTransform(hVRTDS, adfGeoTransform);
     982                 :     }
     983                 : 
     984              13 :     if (bSeparate)
     985                 :     {
     986               2 :         CreateVRTSeparate(hVRTDS);
     987                 :     }
     988                 :     else
     989                 :     {
     990              11 :         CreateVRTNonSeparate(hVRTDS);
     991                 :     }
     992                 : 
     993              13 :     GDALClose(hVRTDS);
     994                 : 
     995              13 :     return CE_None;
     996                 : }
     997                 : 
     998                 : /************************************************************************/
     999                 : /*                        add_file_to_list()                            */
    1000                 : /************************************************************************/
    1001                 : 
    1002                 : static void add_file_to_list(const char* filename, const char* tile_index,
    1003              39 :                              int* pnInputFiles, char*** pppszInputFilenames)
    1004                 : {
    1005                 :    
    1006              39 :     int nInputFiles = *pnInputFiles;
    1007              39 :     char** ppszInputFilenames = *pppszInputFilenames;
    1008                 :     
    1009              39 :     if (EQUAL(CPLGetExtension(filename), "SHP"))
    1010                 :     {
    1011                 : #ifndef OGR_ENABLED
    1012                 :         CPLError(CE_Failure, CPLE_AppDefined, "OGR support needed to read tileindex");
    1013                 :         *pnInputFiles = 0;
    1014                 :         *pppszInputFilenames = NULL;
    1015                 : #else
    1016                 :         OGRDataSourceH hDS;
    1017                 :         OGRLayerH      hLayer;
    1018                 :         OGRFeatureDefnH hFDefn;
    1019                 :         int j, ti_field;
    1020                 : 
    1021               1 :         OGRRegisterAll();
    1022                 :         
    1023                 :         /* Handle GDALTIndex Shapefile as a special case */
    1024               1 :         hDS = OGROpen( filename, FALSE, NULL );
    1025               1 :         if( hDS  == NULL )
    1026                 :         {
    1027                 :             fprintf( stderr, "Unable to open shapefile `%s'.\n", 
    1028               0 :                     filename );
    1029               0 :             exit(2);
    1030                 :         }
    1031                 :         
    1032               1 :         hLayer = OGR_DS_GetLayer(hDS, 0);
    1033                 : 
    1034               1 :         hFDefn = OGR_L_GetLayerDefn(hLayer);
    1035                 : 
    1036               1 :         for( ti_field = 0; ti_field < OGR_FD_GetFieldCount(hFDefn); ti_field++ )
    1037                 :         {
    1038               1 :             OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, ti_field );
    1039               1 :             const char* pszName = OGR_Fld_GetNameRef(hFieldDefn);
    1040                 : 
    1041               1 :             if (strcmp(pszName, "LOCATION") == 0 && strcmp("LOCATION", tile_index) != 0 )
    1042                 :             {
    1043                 :                 fprintf( stderr, "This shapefile seems to be a tile index of "
    1044               0 :                                 "OGR features and not GDAL products.\n");
    1045                 :             }
    1046               1 :             if( strcmp(pszName, tile_index) == 0 )
    1047               1 :                 break;
    1048                 :         }
    1049                 :     
    1050               1 :         if( ti_field == OGR_FD_GetFieldCount(hFDefn) )
    1051                 :         {
    1052                 :             fprintf( stderr, "Unable to find field `%s' in DBF file `%s'.\n", 
    1053               0 :                     tile_index, filename );
    1054               0 :             return;
    1055                 :         }
    1056                 :     
    1057                 :         /* Load in memory existing file names in SHP */
    1058               1 :         int nTileIndexFiles = OGR_L_GetFeatureCount(hLayer, TRUE);
    1059               1 :         if (nTileIndexFiles == 0)
    1060                 :         {
    1061               0 :             fprintf( stderr, "Tile index %s is empty. Skipping it.\n", filename);
    1062               0 :             return;
    1063                 :         }
    1064                 :         
    1065                 :         ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
    1066               1 :                               sizeof(char*) * (nInputFiles+nTileIndexFiles));
    1067               5 :         for(j=0;j<nTileIndexFiles;j++)
    1068                 :         {
    1069               4 :             OGRFeatureH hFeat = OGR_L_GetNextFeature(hLayer);
    1070                 :             ppszInputFilenames[nInputFiles++] =
    1071               4 :                     CPLStrdup(OGR_F_GetFieldAsString(hFeat, ti_field ));
    1072               4 :             OGR_F_Destroy(hFeat);
    1073                 :         }
    1074                 : 
    1075               1 :         OGR_DS_Destroy( hDS );
    1076                 : #endif
    1077                 :     }
    1078                 :     else
    1079                 :     {
    1080                 :         ppszInputFilenames = (char**)CPLRealloc(ppszInputFilenames,
    1081              38 :                                                  sizeof(char*) * (nInputFiles+1));
    1082              38 :         ppszInputFilenames[nInputFiles++] = CPLStrdup(filename);
    1083                 :     }
    1084                 : 
    1085              39 :     *pnInputFiles = nInputFiles;
    1086              39 :     *pppszInputFilenames = ppszInputFilenames;
    1087                 : }
    1088                 : 
    1089                 : /************************************************************************/
    1090                 : /*                                main()                                */
    1091                 : /************************************************************************/
    1092                 : 
    1093              14 : int main( int nArgc, char ** papszArgv )
    1094                 : 
    1095                 : {
    1096              14 :     const char *tile_index = "location";
    1097              14 :     const char *resolution = NULL;
    1098              14 :     int nInputFiles = 0;
    1099              14 :     char ** ppszInputFilenames = NULL;
    1100              14 :     const char * pszOutputFilename = NULL;
    1101                 :     int i, iArg;
    1102              14 :     int bSeparate = FALSE;
    1103              14 :     int bAllowProjectionDifference = FALSE;
    1104              14 :     int bQuiet = FALSE;
    1105              14 :     GDALProgressFunc pfnProgress = NULL;
    1106              14 :     double we_res = 0, ns_res = 0;
    1107              14 :     double xmin = 0, ymin = 0, xmax = 0, ymax = 0;
    1108              14 :     int bAddAlpha = FALSE;
    1109              14 :     int bForceOverwrite = FALSE;
    1110              14 :     int bHideNoData = FALSE;
    1111              14 :     const char* pszSrcNoData = NULL;
    1112              14 :     const char* pszVRTNoData = NULL;
    1113                 : 
    1114              14 :     GDALAllRegister();
    1115                 : 
    1116              14 :     nArgc = GDALGeneralCmdLineProcessor( nArgc, &papszArgv, 0 );
    1117              14 :     if( nArgc < 1 )
    1118               0 :         exit( -nArgc );
    1119                 : 
    1120                 : /* -------------------------------------------------------------------- */
    1121                 : /*      Parse commandline.                                              */
    1122                 : /* -------------------------------------------------------------------- */
    1123              70 :     for( iArg = 1; iArg < nArgc; iArg++ )
    1124                 :     {
    1125              57 :         if( EQUAL(papszArgv[iArg], "--utility_version") )
    1126                 :         {
    1127                 :             printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
    1128               1 :                    papszArgv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
    1129               1 :             return 0;
    1130                 :         }
    1131              56 :         else if( EQUAL(papszArgv[iArg],"-tileindex") &&
    1132                 :                  iArg + 1 < nArgc)
    1133                 :         {
    1134               0 :             tile_index = papszArgv[++iArg];
    1135                 :         }
    1136              56 :         else if( EQUAL(papszArgv[iArg],"-resolution") &&
    1137                 :                  iArg + 1 < nArgc)
    1138                 :         {
    1139               0 :             resolution = papszArgv[++iArg];
    1140                 :         }
    1141              57 :         else if( EQUAL(papszArgv[iArg],"-input_file_list") &&
    1142                 :                  iArg + 1 < nArgc)
    1143                 :         {
    1144               1 :             const char* input_file_list = papszArgv[++iArg];
    1145               1 :             FILE* f = VSIFOpen(input_file_list, "r");
    1146               1 :             if (f)
    1147                 :             {
    1148               4 :                 while(1)
    1149                 :                 {
    1150               5 :                     const char* filename = CPLReadLine(f);
    1151               5 :                     if (filename == NULL)
    1152               1 :                         break;
    1153                 :                     add_file_to_list(filename, tile_index,
    1154               4 :                                      &nInputFiles, &ppszInputFilenames);
    1155                 :                 }
    1156               1 :                 VSIFClose(f);
    1157                 :             }
    1158                 :         }
    1159              55 :         else if ( EQUAL(papszArgv[iArg],"-separate") )
    1160                 :         {
    1161               2 :             bSeparate = TRUE;
    1162                 :         }
    1163              53 :         else if ( EQUAL(papszArgv[iArg],"-allow_projection_difference") )
    1164                 :         {
    1165               0 :             bAllowProjectionDifference = TRUE;
    1166                 :         }
    1167                 :         /* Alternate syntax for output file */
    1168              53 :         else if( EQUAL(papszArgv[iArg],"-o")  &&
    1169                 :                  iArg + 1 < nArgc)
    1170                 :         {
    1171               0 :             pszOutputFilename = papszArgv[++iArg];
    1172                 :         }
    1173              53 :         else if ( EQUAL(papszArgv[iArg],"-q") || EQUAL(papszArgv[iArg],"-quiet") )
    1174                 :         {
    1175               0 :             bQuiet = TRUE;
    1176                 :         }
    1177              55 :         else if ( EQUAL(papszArgv[iArg],"-tr") && iArg + 2 < nArgc)
    1178                 :         {
    1179               2 :             we_res = CPLAtofM(papszArgv[++iArg]);
    1180               2 :             ns_res = CPLAtofM(papszArgv[++iArg]);
    1181                 :         }
    1182              53 :         else if ( EQUAL(papszArgv[iArg],"-te") && iArg + 4 < nArgc)
    1183                 :         {
    1184               2 :             xmin = CPLAtofM(papszArgv[++iArg]);
    1185               2 :             ymin = CPLAtofM(papszArgv[++iArg]);
    1186               2 :             xmax = CPLAtofM(papszArgv[++iArg]);
    1187               2 :             ymax = CPLAtofM(papszArgv[++iArg]);
    1188                 :         }
    1189              49 :         else if ( EQUAL(papszArgv[iArg],"-addalpha") )
    1190                 :         {
    1191               0 :             bAddAlpha = TRUE;
    1192                 :         }
    1193              49 :         else if ( EQUAL(papszArgv[iArg],"-hidenodata") )
    1194                 :         {
    1195               0 :             bHideNoData = TRUE;
    1196                 :         }
    1197              49 :         else if ( EQUAL(papszArgv[iArg],"-overwrite") )
    1198                 :         {
    1199               0 :             bForceOverwrite = TRUE;
    1200                 :         }
    1201              50 :         else if ( EQUAL(papszArgv[iArg],"-srcnodata") && iArg + 1 < nArgc)
    1202                 :         {
    1203               1 :             pszSrcNoData = papszArgv[++iArg];
    1204                 :         }
    1205              48 :         else if ( EQUAL(papszArgv[iArg],"-vrtnodata") && iArg + 1 < nArgc)
    1206                 :         {
    1207               0 :             pszVRTNoData = papszArgv[++iArg];
    1208                 :         }
    1209              48 :         else if ( papszArgv[iArg][0] == '-' )
    1210                 :         {
    1211               0 :             printf("Unrecognized option : %s\n", papszArgv[iArg]);
    1212               0 :             exit(-1);
    1213                 :         }
    1214              48 :         else if( pszOutputFilename == NULL )
    1215                 :         {
    1216              13 :             pszOutputFilename = papszArgv[iArg];
    1217                 :         }
    1218                 :         else
    1219                 :         {
    1220                 :             add_file_to_list(papszArgv[iArg], tile_index,
    1221              35 :                              &nInputFiles, &ppszInputFilenames);
    1222                 :         }
    1223                 :     }
    1224                 : 
    1225              13 :     if( pszOutputFilename == NULL || nInputFiles == 0 )
    1226               0 :         Usage();
    1227                 : 
    1228              13 :     if (!bQuiet)
    1229              13 :         pfnProgress = GDALTermProgress;
    1230                 :        
    1231                 :     /* Avoid overwriting a non VRT dataset if the user did not put the */
    1232                 :     /* filenames in the right order */
    1233                 :     VSIStatBuf sBuf;
    1234              13 :     if (!bForceOverwrite)
    1235                 :     {
    1236              13 :         int bExists = (VSIStat(pszOutputFilename, &sBuf) == 0);
    1237              13 :         if (bExists)
    1238                 :         {
    1239               6 :             GDALDriverH hDriver = GDALIdentifyDriver( pszOutputFilename, NULL );
    1240               6 :             if (hDriver && !EQUAL(GDALGetDriverShortName(hDriver), "VRT"))
    1241                 :             {
    1242                 :                 fprintf(stderr,
    1243                 :                         "'%s' is an existing GDAL dataset managed by %s driver.\n"
    1244                 :                         "There is an high chance you did not put filenames in the right order.\n"
    1245                 :                         "If you want to overwrite %s, add -overwrite option to the command line.\n\n",
    1246               0 :                         pszOutputFilename, GDALGetDriverShortName(hDriver), pszOutputFilename);
    1247               0 :                 Usage();
    1248                 :             }
    1249                 :         }
    1250                 :     }
    1251                 :     
    1252              13 :     if (we_res != 0 && ns_res != 0 &&
    1253                 :         resolution != NULL && !EQUAL(resolution, "user"))
    1254                 :     {
    1255               0 :         fprintf(stderr, "-tr option is not compatible with -resolution %s\n", resolution);
    1256               0 :         Usage();
    1257                 :     }
    1258                 :     
    1259              13 :     if (bAddAlpha && bSeparate)
    1260                 :     {
    1261               0 :         fprintf(stderr, "-addalpha option is not compatible with -separate\n");
    1262               0 :         Usage();
    1263                 :     }
    1264                 :         
    1265              13 :     ResolutionStrategy eStrategy = AVERAGE_RESOLUTION;
    1266              26 :     if ( resolution == NULL || EQUAL(resolution, "user") )
    1267                 :     {
    1268              15 :         if ( we_res != 0 || ns_res != 0)
    1269               2 :             eStrategy = USER_RESOLUTION;
    1270              11 :         else if ( resolution != NULL && EQUAL(resolution, "user") )
    1271                 :         {
    1272               0 :             fprintf(stderr, "-tr option must be used with -resolution user\n");
    1273               0 :             Usage();
    1274                 :         }
    1275                 :     }
    1276               0 :     else if ( EQUAL(resolution, "average") )
    1277               0 :         eStrategy = AVERAGE_RESOLUTION;
    1278               0 :     else if ( EQUAL(resolution, "highest") )
    1279               0 :         eStrategy = HIGHEST_RESOLUTION;
    1280               0 :     else if ( EQUAL(resolution, "lowest") )
    1281               0 :         eStrategy = LOWEST_RESOLUTION;
    1282                 :     else
    1283                 :     {
    1284               0 :         fprintf(stderr, "invalid value (%s) for -resolution\n", resolution);
    1285               0 :         Usage();
    1286                 :     }
    1287                 :     
    1288                 :     /* If -srcnodata is specified, use it as the -vrtnodata if the latter is not */
    1289                 :     /* specified */
    1290              13 :     if (pszSrcNoData != NULL && pszVRTNoData == NULL)
    1291               1 :         pszVRTNoData = pszSrcNoData;
    1292                 : 
    1293                 :     VRTBuilder oBuilder(pszOutputFilename, nInputFiles, ppszInputFilenames,
    1294                 :                         eStrategy, we_res, ns_res, xmin, ymin, xmax, ymax,
    1295                 :                         bSeparate, bAllowProjectionDifference, bAddAlpha, bHideNoData,
    1296              13 :                         pszSrcNoData, pszVRTNoData);
    1297                 : 
    1298              13 :     oBuilder.Build(pfnProgress, NULL);
    1299                 :     
    1300              55 :     for(i=0;i<nInputFiles;i++)
    1301                 :     {
    1302              42 :         CPLFree(ppszInputFilenames[i]);
    1303                 :     }
    1304              13 :     CPLFree(ppszInputFilenames);
    1305                 : 
    1306                 : 
    1307              13 :     CSLDestroy( papszArgv );
    1308              13 :     GDALDumpOpenDatasets( stderr );
    1309              13 :     GDALDestroyDriverManager();
    1310                 : #ifdef OGR_ENABLED
    1311              13 :     OGRCleanupAll();
    1312                 : #endif
    1313                 : 
    1314              13 :     return 0;
    1315                 : }

Generated by: LTP GCOV extension version 1.5