Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
- (WarpResultC*)warpImageWithGeoTransform:(NSArray<NSNumber*>*)geoTransformArray sourceFile:(NSString*)inFilepath destinationFile:(NSString*)outFilepath
{
GDALAllRegister();
GDALDriverH hDriver;
GDALDataType eDT;
GDALDatasetH hDstDS;
GDALDatasetH hSrcDS;
// Open the source file.
hSrcDS = GDALOpen( inFilepath.UTF8String, GA_ReadOnly );
CPLAssert( hSrcDS != NULL );
// Set the GeoTransform on the source image
double geoTransform[] = { geoTransformArray[0].doubleValue, geoTransformArray[1].doubleValue, geoTransformArray[2].doubleValue, geoTransformArray[3].doubleValue, -geoTransformArray[4].doubleValue, -geoTransformArray[5].doubleValue };
GDALSetGeoTransform(hSrcDS, geoTransform);
// Create output with same datatype as first input band.
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));
// Get output driver (GeoTIFF format)
hDriver = GDALGetDriverByName( "GTiff" );
CPLAssert( hDriver != NULL );
// Create a transformer that maps from source pixel/line coordinates
// to destination georeferenced coordinates (not destination
// pixel line). We do that by omitting the destination dataset
// handle (setting it to NULL).
void *hTransformArg = GDALCreateGenImgProjTransformer( hSrcDS, NULL, NULL, NULL, FALSE, 0, 1 );
CPLAssert( hTransformArg != NULL );
// Get approximate output georeferenced bounds and resolution for file.
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
CPLErr eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform, hTransformArg, adfDstGeoTransform, &nPixels, &nLines );
CPLAssert( eErr == CE_None );
GDALDestroyGenImgProjTransformer( hTransformArg );
// Create the output file.
hDstDS = GDALCreate( hDriver, outFilepath.UTF8String, nPixels, nLines, 4, eDT, NULL );
CPLAssert( hDstDS != NULL );
// Write out the projection definition.
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
// Copy the color table, if required.
GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS, 1) );
if( hCT != NULL )
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS, 1), hCT );
// Setup warp options.
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = hDstDS;
/* -------------------------------------------------------------------- */
/* Do we have a source alpha band? */
/* -------------------------------------------------------------------- */
bool enableSrcAlpha = GDALGetRasterColorInterpretation( GDALGetRasterBand(hSrcDS, GDALGetRasterCount(hSrcDS) )) == GCI_AlphaBand;
if(enableSrcAlpha) { printf( "Using band %d of source image as alpha.\n", GDALGetRasterCount(hSrcDS) ); }
/* -------------------------------------------------------------------- */
/* Setup band mapping. */
/* -------------------------------------------------------------------- */
if(enableSrcAlpha)
psWarpOptions->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
else
psWarpOptions->nBandCount = GDALGetRasterCount(hSrcDS);
psWarpOptions->panSrcBands = (int *) CPLMalloc(psWarpOptions->nBandCount*sizeof(int));
psWarpOptions->panDstBands = (int *) CPLMalloc(psWarpOptions->nBandCount*sizeof(int));
for( int i = 0; i < psWarpOptions->nBandCount; i++ )
{
psWarpOptions->panSrcBands[i] = i+1;
psWarpOptions->panDstBands[i] = i+1;
}
/* -------------------------------------------------------------------- */
/* Setup alpha bands used if any. */
/* -------------------------------------------------------------------- */
if( enableSrcAlpha )
psWarpOptions->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
psWarpOptions->nDstAlphaBand = GDALGetRasterCount(hDstDS);
psWarpOptions->pfnProgress = GDALTermProgress;
// Establish reprojection transformer.
psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer( hSrcDS, NULL, hDstDS, NULL, FALSE, 0.0, 1 );
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
// Initialize and execute the warp operation.
GDALWarpOperation oOperation;
oOperation.Initialize( psWarpOptions );
CPLErr warpError = oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) );
CPLAssert( warpError == CE_None );
GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
GDALDestroyWarpOptions( psWarpOptions );
GDALClose( hDstDS );
GDALClose( hSrcDS );
WarpResultC* warpResultC = [WarpResultC new];
warpResultC.geoTransformValues = @[@(adfDstGeoTransform[0]), @(adfDstGeoTransform[1]), @(adfDstGeoTransform[2]), @(adfDstGeoTransform[3]), @(adfDstGeoTransform[4]), @(adfDstGeoTransform[5])];
warpResultC.newX = nPixels;
warpResultC.newY = nLines;
return warpResultC;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.