Thanks Even, I appreciate it. kss
/** * * Kyle Shannon * [email protected] * */ On Tue, May 10, 2011 at 15:41, Even Rouault <[email protected]>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 > > * [email protected] > > * > > */ > > > > On Tue, May 10, 2011 at 03:51, Even Rouault <even.rouault@mines- > paris.org>wrote: > > > Selon Kyle Shannon <[email protected]>: > > > > 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 > > > > * [email protected] > > > > * > > > > */ > > > > > > > > On Mon, May 9, 2011 at 15:43, Even Rouault > > > > > > > > <[email protected]>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. >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
