<div dir="ltr">I have those nodata-values there since the same code path is used for elevation data and such, but it had no effect to remove it.<div><br></div><div>I have narrowed it down to the case where I use an overview as input. I have some very large images with overviews, so I find the overview closest matching to the output resolution of my image and make a VRT of the overview level like this:</div><div><br></div><div><div><span class="" style="white-space:pre">                                </span>vrtDS = vrtDriver->Create("", ovrW, ovrH, _numBands, _dataType, NULL);</div><div><span class="" style="white-space:pre">                                </span>GDALSetProjection(vrtDS, _srcProjectionWKT.c_str());</div><div><br></div><div><span class="" style="white-space:pre">                              </span>geoTransform[XFM_PIXEL_SIZE_EW] = scaleX*basePixelSizeX;</div><div><span class="" style="white-space:pre">                           </span>geoTransform[XFM_PIXEL_SIZE_NS] = scaleY*basePixelSizeY;</div><div><span class="" style="white-space:pre">                           </span>vrtDS->SetGeoTransform(geoTransform);</div><div><br></div><div><span class="" style="white-space:pre">                          </span>for(int i = 1; i <= _numBands; i++)</div><div><span class="" style="white-space:pre">                             </span>{</div><div><span class="" style="white-space:pre">                                  </span>GDALRasterBand* srcRootBand = _ds->GetRasterBand(i);</div><div><span class="" style="white-space:pre">                                    </span>GDALRasterBand* srcBand = srcRootBand->GetOverview(overviewId);</div><div><span class="" style="white-space:pre">                                 </span>VRTSourcedRasterBand* vrtBand = (VRTSourcedRasterBand*)vrtDS->GetRasterBand(i);</div><div><span class="" style="white-space:pre">                                 </span>vrtBand->AddSimpleSource(srcBand);</div><div><span class="" style="white-space:pre">                              </span>}</div><div><span class="" style="white-space:pre">                          </span>srcDS = vrtDS;</div></div><div><br></div><div>The strange thing is that this used to work when I had the same mistake as you pointed out, where I added the bands twice, i.e the overview had 6 bands, but only three or four of them were used in the warp operation. When fixing the code, I got this problem.</div><div><br></div><div>Is this the best way to do it btw, or is it a better way to make use of overviews when cutting out and reprojecting parts of a big image?</div><div><br></div><div>Thanks again for all your help :-)</div><div><br></div><div>- Thomas</div><div><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jun 12, 2015 at 1:35 PM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@spatialys.com" target="_blank">even.rouault@spatialys.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Selon Thomas Sevaldrud <<a href="mailto:thomas@silentwings.no">thomas@silentwings.no</a>>:<br>
<div><div class="h5"><br>
> Thanks again! The problem was in the warper setup, it used only the first<br>
> band of the image. This also revealed a bug in another part of my code<br>
> related to VRT's and Warp, since I had added double sets of bands there as<br>
> well. However, when fixing this, some old code stopped working :-)<br>
><br>
> Basically, my problem is that I want to reproject an RGB image which may be<br>
> the VRT I made above, or another RGB(A) image, into  an output RGBA. If I<br>
> don't have an input alpha channel, I still want the output to have alpha,<br>
> and I want any borders that result from the reprojection to have alpha = 0.<br>
> What is the correct way of setting up this? Currently I set up the bands<br>
> like this:<br>
><br>
> warpOptions->nBandCount = srcNumBands;<br>
><br>
> warpOptions->panSrcBands =<br>
> (int*)CPLMalloc(warpOptions->nBandCount*sizeof(int));<br>
> warpOptions->panDstBands =<br>
> (int*)CPLMalloc(warpOptions->nBandCount*sizeof(int));<br>
><br>
> warpOptions->padfSrcNoDataReal =<br>
> (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));<br>
> warpOptions->padfSrcNoDataImag =<br>
> (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));<br>
><br>
> warpOptions->padfDstNoDataReal =<br>
> (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));<br>
> warpOptions->padfDstNoDataImag =<br>
> (double*)CPLMalloc(warpOptions->nBandCount*sizeof(double));<br>
><br>
> for (int j = 0; j < warpOptions->nBandCount; j++)<br>
> {<br>
> warpOptions->panSrcBands[j] = j+1;<br>
> warpOptions->panDstBands[j] = j+1;<br>
><br>
> warpOptions->padfSrcNoDataReal[j] = _srcNodata;<br>
> warpOptions->padfSrcNoDataImag[j] = 0.0;<br>
><br>
> warpOptions->padfDstNoDataReal[j] = _nodata;<br>
> warpOptions->padfDstNoDataImag[j] = 0.0;<br>
> }<br>
><br>
> if (srcNumBands == 4)<br>
> warpOptions->nSrcAlphaBand = 4;<br>
><br>
> if (dstNumBands == 4)<br>
> warpOptions->nDstAlphaBand = 4;<br>
><br>
> char **papszWarpOptions = NULL;<br>
> papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "NO_DATA"<br>
> );<br>
><br>
> In some cases this works, but in others I get only a black image after the<br>
> warp. So my question is basically if the above code should have worked. My<br>
> assumption here is that the GDAL Warper will regard all input pixels as<br>
> having full opacity when no alpha channel is specified, and that output<br>
> alpha is initialized to 0. Are these valid assumptions, or do I need to<br>
> create an alpha channel in the input image as well?<br>
<br>
</div></div>I don't see anything fundamentally wrong in the above. If your source image is<br>
fully opaque, you don't need to set padfSrcNoDataReal and padfSrcNoDataImag. And<br>
if you have an output alpha channel, you don't need to set padfDstNoDataReal and<br>
padfDstNoDataImag too. In that cas you would just do papszWarpOptions =<br>
CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0" ); I'm not sure this will<br>
change the end result however.<br>
<div class="HOEnZb"><div class="h5"><br>
><br>
> - Thomas<br>
><br>
> On Thu, Jun 11, 2015 at 5:25 PM, Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>><br>
> wrote:<br>
><br>
> > Selon Thomas Sevaldrud <<a href="mailto:thomas@silentwings.no">thomas@silentwings.no</a>>:<br>
> ><br>
> > > Great, thanks!<br>
> > ><br>
> > > I tried this, but got only a red image as result, so I guess only the<br>
> > first<br>
> > > channel was used.<br>
> > ><br>
> > > This is the relevant code, where _ds is the input paletted data set<br>
> > ><br>
> > > vrtDS = vrtDriver->Create("", origW, origH, 3, GDT_Byte, NULL);<br>
> > > double geoTransform[6];<br>
> > > _ds->GetGeoTransform(geoTransform);<br>
> > > vrtDS->SetGeoTransform(geoTransform);<br>
> > > GDALSetProjection(vrtDS, _srcProjectionWKT.c_str());<br>
> > ><br>
> > > GDALRasterBand* srcBand = _ds->GetRasterBand(1);<br>
> > ><br>
> > > vrtDS->AddBand(GDT_Byte, NULL);<br>
> > > vrtDS->AddBand(GDT_Byte, NULL);<br>
> > > vrtDS->AddBand(GDT_Byte, NULL);<br>
> ><br>
> > You don't need the above 3 AddBand() since you called Create() with 3<br>
> > already<br>
> ><br>
> > > VRTSourcedRasterBand* vrtRBand =<br>
> > > (VRTSourcedRasterBand*)vrtDS->GetRasterBand(1);<br>
> > > VRTSourcedRasterBand* vrtGBand =<br>
> > > (VRTSourcedRasterBand*)vrtDS->GetRasterBand(2);<br>
> > > VRTSourcedRasterBand* vrtBBand =<br>
> > > (VRTSourcedRasterBand*)vrtDS->GetRasterBand(3);<br>
> > > vrtRBand->AddComplexSource(srcBand, -1, -1, -1, -1, -1, -1, -1, -1, 0.0,<br>
> > > 1.0, VRT_NODATA_UNSET, 1);<br>
> > > vrtGBand->AddComplexSource(srcBand, -1, -1, -1, -1, -1, -1, -1, -1, 0.0,<br>
> > > 1.0, VRT_NODATA_UNSET, 2);<br>
> > > vrtBBand->AddComplexSource(srcBand, -1, -1, -1, -1, -1, -1, -1, -1, 0.0,<br>
> > > 1.0, VRT_NODATA_UNSET, 3);<br>
> > ><br>
> > > Is there anything I have missed in the setup here?<br>
> ><br>
> > Looks good otherwise. Did you check that you configured properly the<br>
> > warper to<br>
> > use the 3 bands of the VRT ?<br>
> > ><br>
> > > - Thomas<br>
> > ><br>
> > ><br>
> > ><br>
> > > On Tue, Jun 9, 2015 at 1:44 PM, Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a><br>
> > ><br>
> > > wrote:<br>
> > ><br>
> > > > Thomas,<br>
> > > > ><br>
> > > > > I'm using the GDALWarp api from C++ to reproject and cut various<br>
> > input<br>
> > > > > images. In general this works very well for my purposes, except that<br>
> > for<br>
> > > > > paletted images I have to use NearestNeighbour resampling,<br>
> > > > ><br>
> > > > > I would like to use Bilinear or higher order resampling, and wonder<br>
> > if<br>
> > > > > there is any way to expand a paletted image to true RGB during the<br>
> > warp<br>
> > > > > process, so I can use these resampling methods.<br>
> > > ><br>
> > > > No, you must use an intermediate step before.<br>
> > > ><br>
> > > > ><br>
> > > > > Alternatively, is there a way to do this conversion (from code)<br>
> > before<br>
> > > > the<br>
> > > > > warp without having to run through the entire image in memory?<br>
> > > ><br>
> > > > You can use a in-memory VRT dataset to do on-the-fly expansion to RGB.<br>
> > > ><br>
> > > > With the AddComplexSource() method of VRTSourcedRasterBand class:<br>
> > > ><br>
> > > >     CPLErr         AddComplexSource( GDALRasterBand *poSrcBand,<br>
> > > >                                      int nSrcXOff=-1, int nSrcYOff=-1,<br>
> > > >                                      int nSrcXSize=-1, int<br>
> > nSrcYSize=-1,<br>
> > > >                                      int nDstXOff=-1, int nDstYOff=-1,<br>
> > > >                                      int nDstXSize=-1, int<br>
> > nDstYSize=-1,<br>
> > > >                                      double dfScaleOff=0.0,<br>
> > > >                                      double dfScaleRatio=1.0,<br>
> > > >                                      double dfNoDataValue =<br>
> > > > VRT_NODATA_UNSET,<br>
> > > >                                      int nColorTableComponent = 0);<br>
> > > ><br>
> > > > Set nColorTableComponent to 1,2,3 for respectively R,G,B. This is what<br>
> > > > gdal_translate will do if you use "gdal_translate -of VRT -expand rgb<br>
> > > > pct.tif<br>
> > > > out.vrt"<br>
> > > ><br>
> > > > Even<br>
> > > ><br>
> > > > --<br>
> > > > Spatialys - Geospatial professional services<br>
> > > > <a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
> > > ><br>
> > ><br>
> ><br>
> ><br>
> > --<br>
> > Spatialys - Geospatial professional services<br>
> > <a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
> ><br>
><br>
<br>
<br>
--<br>
Spatialys - Geospatial professional services<br>
<a href="http://www.spatialys.com" rel="noreferrer" target="_blank">http://www.spatialys.com</a><br>
</div></div></blockquote></div><br></div>