[gdal-dev] Setting "INIT_DEST" to nan on windows

[Kyle Shannon]

Thanks Even, I appreciate it.

kss

/**
 *
 * Kyle Shannon
 * ksshannon at gmail.com
 *
 */




On Tue, May 10, 2011 at 15:41, Even Rouault <even.rouault at mines-paris.org>wrote:

> Le mardi 10 mai 2011 17:39:52, Kyle Shannon a écrit :
> > Even, here is the basic flow of my task.  I am opening several
> subdatasets
> > in a for loop, warping them and populating an internal data structure.
> >
> > #include "gdal_priv.h"
> > #include "gdalwarper.h"
> > #include "cpl_conv.h"
> > #include "cpl_port.h"
> > #include <string>
> > #include <vector>
> > int main()
> > {
> > GDALAllRegister();
> > /*
> >  * Variable names for the netcdf forecast file
> >  */
> >
> > std::vector<std::string> varList;
> > varList.push_back( "Temperature_height_above_ground" );
> > varList.push_back( "V-component_of_wind_height_above_ground" );
> > varList.push_back( "U-component_of_wind_height_above_ground" );
> > varList.push_back( "Total_cloud_cover" );
> >
> > std::string temp;
> > GDALDataset* srcDS;
> > GDALDataset* wrpDS;
> > GDALWarpOptions* psWarpOptions;
> > std::string srcWkt, dstWkt;
> > dstWkt = "PROJCS[\"WGS 84 / UTM zone 12N\",GEOGCS[\"WGS
> > 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS
> >
> 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174
> >
> 532925199433]],UNIT[\"metre\",1],PROJECTION[\"Transverse_Mercator\"],PARAME
> >
> TER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-111],PARAMETE
> >
> R[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"
> > false_northing\",0]]";
> >
> > /*
> >  * Loop over variable names and open sub datasets
> >  */
> >
> > for( int i = 0;i < varList.size();i++ ) {
> >
> >     temp = "NETCDF:test.nc:" + varList[i];
> >
> >     srcDS = (GDALDataset*)GDALOpenShared( temp.c_str(), GA_ReadOnly );
> >
> >     srcWkt = srcDS->GetProjectionRef();
> >
> >     /*
> >      * Grab the first band to get the nodata value for the variable,
> >      * assume all bands have the same ndv
> >      */
> >     GDALRasterBand *poBand = srcDS->GetRasterBand( 1 );
> >     int pbSuccess;
> >     double dfNoData = poBand->GetNoDataValue( &pbSuccess );
> >     psWarpOptions = GDALCreateWarpOptions();
> >
> >     int nBandCount = srcDS->GetRasterCount();
> >     psWarpOptions->padfDstNoDataReal =
> >         (double*) CPLMalloc( sizeof( double ) * nBandCount );
> >     psWarpOptions->padfDstNoDataImag =
> >         (double*) CPLMalloc( sizeof( double ) * nBandCount );
> >     for( int b = 0;b < srcDS->GetRasterCount();b++ ) {
> >         psWarpOptions->padfDstNoDataReal[b] = dfNoData;
> >         psWarpOptions->padfDstNoDataImag[b] = dfNoData;
> >     }
> >     if( pbSuccess == false )
> >         dfNoData = -9999.0;
> >     psWarpOptions->papszWarpOptions =
> >     CSLSetNameValue( psWarpOptions->papszWarpOptions,
> >                      "INIT_DEST", "NO_DATA" );
> >
> >     //set the dskWkt with an internal structure
> >
> >     wrpDS = (GDALDataset*) GDALAutoCreateWarpedVRT( srcDS,
> srcWkt.c_str(),
> >                                                     dstWkt.c_str(),
> >                                                     GRA_NearestNeighbour,
> >                                                     1.0, psWarpOptions );
> >
> >     //copy warped data to internal structure
> >
> >     GDALDestroyWarpOptions( psWarpOptions );
> >     GDALClose((GDALDatasetH) srcDS );
> >     GDALClose((GDALDatasetH) wrpDS );
> > }
> > }
> >
>
> Of course, Valgrind has been helpfull once again (too bad you're on Windows
> ;-))
>
> You just forgot : psWarpOptions->nBandCount = nBandCount;
>
> > I am sure I am missing something.
> >
> > kss
> >
> > /**
> >  *
> >  * Kyle Shannon
> >  * ksshannon at gmail.com
> >  *
> >  */
> >
> > On Tue, May 10, 2011 at 03:51, Even Rouault <even.rouault at mines-
> paris.org>wrote:
> > > Selon Kyle Shannon <ksshannon at gmail.com>:
> > > > I did have to set the padfNoDataImag array as well, so I am one step
> > > > closer.  My new issue is ownership of my psWarpOptions.  I assume
> that
> > >
> > > since
> > >
> > > > I create it, and the warp api tutorial implies that is it mine, I
> > > > destroy
> > >
> > > it
> > >
> > > > with GDALDestroyWarpOptions().  My problem is that
> AutoCreateWarpedVRT
> > > > appears to be stepping on my warp options when it destroys it's own
> > > > copy.
> > > >
> > > > I call AutoCreateWarpedVRT() and that function calls
> > > > GDALCloneWarpOptions().  When AutoCreateWarpedVRT destroys the copy
> it
> > >
> > > owns,
> > >
> > > > it appears to be stepping on my copy, and when I call
> > > > GDALDestroyWarpOptions(), I get a seg fault.  Am I doing something
> > > > wrong?
> > >
> > > Difficult to say without seeing your code. Could you post a small
> > > self-sufficient code snippet that demonstrates your issue ?
> > >
> > > By reviewing the code, one potential issue would be if you try to free
> > > the pTransformerArg field by yourself, as it is freed by the destructor
> > > of the VRTWarpedDataset. But if you just use GDALDestroyWarpOptions(),
> > > it should work
> > > OK
> > >
> > > > kss
> > > >
> > > >
> > > >
> > > > /**
> > > >
> > > >  *
> > > >  * Kyle Shannon
> > > >  * ksshannon at gmail.com
> > > >  *
> > > >  */
> > > >
> > > > On Mon, May 9, 2011 at 15:43, Even Rouault
> > > >
> > > > <even.rouault at mines-paris.org>wrote:
> > > > > Le lundi 09 mai 2011 23:30:05, Kyle Shannon a écrit :
> > > > > > Thanks Even,
> > > > > > If that is the case, are there other values in psWarpOptions I
> > > > > > *have*
> > >
> > > to
> > >
> > > > > > set?  I can't find docs discussing that.
> > > > >
> > > > > I don't think so, but you'll just have to try... I'm not  surprised
> > > > > the docs
> > > > > don't talk about NaN. It is a corner case, and likely not really
> well
> > > > > supported in all places.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20110510/0089f4d9/attachment-0001.html