Last active
April 30, 2016 20:34
-
-
Save koloritcm/9d73aa4358f3dafe2fd6808f7d4e9f41 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
- (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