Unable to warp through a VRTDataset with a destination alpha band

[tcm0116]

I wrote the function below to allow me to warp an image through a VRTDataset in order to capitalize on the advanced capabilities of the VRT filter. However, when attempting to warp a grayscale iamge (with a single band), the CreateCopy fails with an error message similar to "GetBlockRef failed at X block offset 2, Y block offset 0". I haven't tried other source image types, but I suspect the result would be the same.

P.S. This code is almost identical to GDALAutoCreateWarpedVRT with the addition of the alpha band.

GDALDatasetH WarpImage(GDALDatasetH hSrcDS, const char* dstFile,
   const char** dstWKT, GDALResampleAlg alg, CPLStringList options)
{
   GDALWarpOptions *psWO;
   int i;

   /* -------------------------------------------------------------------- */
   /*      Populate the warp options.                                      */
   /* -------------------------------------------------------------------- */
   psWO = GDALCreateWarpOptions();

   psWO->eResampleAlg = alg;

   psWO->hSrcDS = hSrcDS;

   psWO->nBandCount = GDALGetRasterCount( psWO->hSrcDS );
   psWO->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);
   psWO->panDstBands = (int *) CPLMalloc(sizeof(int) * psWO->nBandCount);

   for( i = 0; i < psWO->nBandCount; i++ )
   {
      psWO->panSrcBands[i] = i+1;
      psWO->panDstBands[i] = i+1;
   }

   /* -------------------------------------------------------------------- */
   /*      Create the transformer.                                         */
   /* -------------------------------------------------------------------- */
   psWO->pfnTransformer = GDALGenImgProjTransform;
   psWO->pTransformerArg = 
      GDALCreateGenImgProjTransformer( psWO->hSrcDS, 
      GDALGetProjectionRef( psWO->hSrcDS ), 
      NULL, dstWKT,
      TRUE, 1.0, 0 );

   if( psWO->pTransformerArg == NULL )
   {
      GDALDestroyWarpOptions( psWO );
      return NULL;
   }

   /* -------------------------------------------------------------------- */
   /*      Figure out the desired output bounds and resolution.            */
   /* -------------------------------------------------------------------- */
   double adfDstGeoTransform[6];
   int    nDstPixels, nDstLines;
   CPLErr eErr;

   eErr = 
      GDALSuggestedWarpOutput( psWO->hSrcDS, psWO->pfnTransformer, 
      psWO->pTransformerArg, 
      adfDstGeoTransform, &nDstPixels, &nDstLines );

   /* -------------------------------------------------------------------- */
   /*      Update the transformer to include an output geotransform        */
   /*      back to pixel/line coordinates.                                 */
   /*                                                                      */
   /* -------------------------------------------------------------------- */
   GDALSetGenImgProjTransformerDstGeoTransform( 
      psWO->pTransformerArg, adfDstGeoTransform );

   /* -------------------------------------------------------------------- */
   /*      Create the VRTDataset and populate it with bands.               */
   /* -------------------------------------------------------------------- */
   VRTWarpedDataset *poDS = new VRTWarpedDataset( nDstPixels, nDstLines );

   psWO->hDstDS = (GDALDatasetH) poDS;

   poDS->SetGeoTransform( adfDstGeoTransform );
   poDS->SetProjection( dstWKT );

   for( i = 0; i < psWO->nBandCount; i++ )
   {
      VRTWarpedRasterBand *poBand;
      GDALRasterBand *poSrcBand = (GDALRasterBand *) 
         GDALGetRasterBand( psWO->hSrcDS, i+1 );

      poDS->AddBand( poSrcBand->GetRasterDataType(), NULL );

      poBand = (VRTWarpedRasterBand *) poDS->GetRasterBand( i+1 );
      poBand->CopyCommonInfoFrom( poSrcBand );
   }

   /* -------------------------------------------------------------------- */
   /*      Add an alpha band to the dataset                                */
   /* -------------------------------------------------------------------- */
   poDS->AddBand(GDT_Byte);
   poDS->GetRasterBand(poDS->GetRasterCount())->SetColorInterpretation(GCI_AlphaBand);
   psWO->nDstAlphaBand = poDS->GetRasterCount();

   /* -------------------------------------------------------------------- */
   /*      Initialize the warp on the VRTWarpedDataset.                    */
   /*      NOTE: GDALCreateWarpedVRT doesn't check this return value       */
   /* -------------------------------------------------------------------- */
   eErr = poDS->Initialize( psWO );

   GDALDestroyWarpOptions( psWO );

   if( eErr != CE_None )
   {
      GDALClose( poDS );
      return NULL;
   }

   // Create the output dataset
   GDALDriver* pDriver = pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
   GDALDataset* pDataset = pDriver->CreateCopy(dstFile, poDS, FALSE, options, NULL, NULL);

   if (pDataset == NULL)
      printf("WarpImage: Unable to create destination image - %s\n", CPLGetLastErrorMsg());

   GDALClose( poDS );

   return (GDALDatasetH)pDataset;
}

[Even Rouault]

Hum, this turned to be a subtil problem with how warped VRT works with alpha band, that was not triggered when using gdalwarp utility because it takes care to implicetly set the INIT_DEST warping option. I've pushed a fix in trunk (r24190) to do so.

In the meantime, you can easily fix your code by adding :

psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "0");