<font class="Apple-style-span" face="arial, sans-serif" size="3"><span class="Apple-style-span" style="border-collapse: collapse; -webkit-border-horizontal-spacing: 2px; -webkit-border-vertical-spacing: 2px;"><span class="Apple-style-span" style="font-size: 13px; -webkit-border-horizontal-spacing: 0px; -webkit-border-vertical-spacing: 0px; ">Hello,<div>
<br></div><div> I want to project an input Lat/Long image into UTM. I need to do this automatic, so I have to do it in C++. I have found a source code to implement this issue. I have paste below the source code that I have used.</div>
<div> I When this code is running, appears some times the followin error messages: </div><div><br></div><div> "ERROR 1: no system list, errno: 42" and </div><div> "ERROR 6: SetColorTable() not supported for multi-sample TIFF files."</div>
<div> </div><div> What is it means?</div><div>When I open the output image there are some strange black lines crossing the output image that not appear in the input image. Is there any reason for this behaviour?</div>
<div><br></div><div><br></div><div>Many thanks,</div><div><br></div><div>Jorge </div><div><br></div><div><br></div><div><div>//----------------------------------SOUCE CODE----------------------------------------</div><div>
<br></div><div><br></div><div><div>void Reproject()</div><div>{</div><div> // ------------ Example of initializing an image into the GDALDataset object</div><div><br></div><div> GDALDataset *poDataset;</div><div><br></div>
<div> // -------------- Reprojection example</div><div><br></div><div> GDALDatasetH hSrcDS, hDstDS; //, hRefDS;</div><div><br></div><div> // Open input and output files.</div><div><br></div><div> GDALAllRegister();</div>
<div><br></div><div> hSrcDS = GDALOpen( "/home/jjmf/sample_images/t1.09155.AV_MI.143.250m.tif", GA_ReadOnly );</div><div><br></div><div> poDataset = (GDALDataset *) hSrcDS;</div><div> if( poDataset == NULL )</div>
<div> {</div><div> </div><div> cout<<"The input file does not exist: "<< endl;</div><div><br></div><div> exit(1);</div><div> }</div><div><br></div><div> GDALDriverH hDriver;</div>
<div> GDALDataType eDT;</div><div><br></div><div> // Open the source file.</div><div><br></div><div> CPLAssert( hSrcDS != NULL );</div><div><br></div><div> // Create output with same datatype as first input band.</div>
<div><br></div><div> eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));</div><div><br></div><div> // Get output driver (GeoTIFF format)</div><div> hDriver = GDALGetDriverByName("GTiff");</div><div>
CPLAssert( hDriver != NULL );</div><div><br></div><div> // Get Source coordinate system.</div><div><br></div><div> const char *pszSrcWKT = NULL;</div><div> char *pszDstWKT;</div><div><br></div><div> pszSrcWKT = GDALGetProjectionRef( hSrcDS );</div>
<div> CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );</div><div><br></div><div><br></div><div><br></div><div> OGRSpatialReference oSRS;</div><div><br></div><div> oSRS.SetUTM(12, TRUE);</div><div>
<br></div><div> oSRS.SetWellKnownGeogCS("EPSG:4258");</div><div> //oSRS.importFromProj4("+proj=latlo +datum=WGS84");</div><div><br></div><div> oSRS.exportToWkt( &pszDstWKT );</div><div><br>
</div><div> // Create a transformer that maps from source pixel/line coordinates</div><div> // to destination georeferenced coordinates (not destination</div><div> // pixel line). We do that by omitting the destination dataset</div>
<div> // handle (setting it to NULL).</div><div><br></div><div> void *hTransformArg;</div><div><br></div><div> hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 0, 1 );</div>
<div> CPLAssert( hTransformArg != NULL );</div><div><br></div><div> // Get approximate output georeferenced bounds and resolution for file.</div><div><br></div><div> double adfDstGeoTransform[6];</div><div> int nPixels = 0, nLines = 0;</div>
<div> CPLErr eErr;</div><div><br></div><div> eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines );</div><div> CPLAssert( eErr == CE_None );</div>
<div><br></div><div> GDALDestroyGenImgProjTransformer( hTransformArg );</div><div><br></div><div> // Create the output file.</div><div> char *apszOptions[] = {"INTERLEAVE=PIXEL", NULL};</div><div> </div>
<div> hDstDS = GDALCreate(hDriver, "/home/jjmf/sample_images/out_new.tif", nPixels, nLines, GDALGetRasterCount(hSrcDS), eDT, apszOptions);</div><div><br></div><div> CPLAssert( hDstDS != NULL );</div><div><br>
</div><div> // Write out the projection definition.</div><div><br></div><div> GDALSetProjection(hDstDS, pszDstWKT);</div><div> GDALSetGeoTransform(hDstDS, adfDstGeoTransform);</div><div><br></div><div> // Copy the color table, if required.</div>
<div><br></div><div> GDALColorTableH hCT;</div><div><br></div><div> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,1));</div><div> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,1), hCT);</div><div><br>
</div><div><br></div><div> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,2));</div><div> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,2), hCT);</div><div><br></div><div><br></div><div> hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS,3));</div>
<div> GDALSetRasterColorTable(GDALGetRasterBand(hDstDS,3), hCT);</div><div><br></div><div><br></div><div> cout << "Setup warp options ..."<<endl;</div><div> // Setup warp options.</div><div><br>
</div><div> GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();</div><div><br></div><div> psWarpOptions->hSrcDS = hSrcDS;</div><div> psWarpOptions->hDstDS = hDstDS;</div><div> psWarpOptions->nBandCount = 3;</div>
<div> psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);</div><div> psWarpOptions->panSrcBands[0] = 1;</div><div> psWarpOptions->panSrcBands[1] = 2;</div><div> psWarpOptions->panSrcBands[2] = 3;</div>
<div> psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);</div><div> psWarpOptions->panDstBands[0] = 1;</div><div> psWarpOptions->panDstBands[1] = 2;</div><div> psWarpOptions->panDstBands[2] = 3;</div>
<div><br></div><div> psWarpOptions->pfnProgress = GDALTermProgress;</div><div><br></div><div> cout << "Establishing transformer ..."<<endl;</div><div> // Establish reprojection transformer.</div>
<div><br></div><div> psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(hSrcDS, GDALGetProjectionRef(hSrcDS), hDstDS, GDALGetProjectionRef(hDstDS), FALSE, 0.0, 1);</div><div> psWarpOptions->pfnTransformer = GDALGenImgProjTransform;</div>
<div><br></div><div> // Initialize and execute the warp operation.</div><div><br></div><div> cout << "Executing ..."<<endl;</div><div> GDALWarpOperation oOperation;</div><div><br></div><div> // Initialize transformation operation</div>
<div> oOperation.Initialize(psWarpOptions);</div><div><br></div><div> // Perform the transformation</div><div> oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(hDstDS), GDALGetRasterYSize(hDstDS));</div><div><br></div>
<div> // Clean up the intermediate objects</div><div> GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);</div><div> GDALDestroyWarpOptions(psWarpOptions);</div><div><br></div><div> </div><div> cout <<"Done ..."<<endl;</div>
<div> </div><div> // Close the input and output files</div><div> GDALClose(hDstDS);</div><div> GDALClose(hSrcDS);</div><div><br></div><div><br></div><div><br></div><div>}</div></div><div><br></div><div>//--------------------------------------------------------------------------------------------------------------------</div>
<div><br></div><div><br></div><div><br></div></div></span></span></font>