/******************************************************************************
 * $Id: gdalwarp.cpp,v 1.9 2005/02/18 21:22:53 jonathan Exp $
 *
 * Project:  High Performance Image Reprojector
 * Purpose:  Test program for high performance warper API.
 * Author:   Frank Warmerdam <warmerdam@pobox.com>
 *
 ******************************************************************************
 * Copyright (c) 2002, i3 - information integration and imaging 
 *                          Fort Collin, CO
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ******************************************************************************
 *
 * $Log: gdalwarp.cpp,v $
 * Revision 1.9  2005/02/18 21:22:53  jonathan
 * Optimize the loop in gdalwarp which builds a mask. Handle the majority of an
 * image in a loop, creating 8 bits at a time. Later, handle the edge case where
 * less than 8 bits are packed.
 *
 * Revision 1.8  2005/02/18 14:54:17  jonathan
 * Refactored baserenderer.py and renderer.py to remove baserenderer.py's
 * dependencies on wxPython. Added a new method projected_raster_layer()
 * that returns a raster layer image in projected space. This must be
 * implemented in classes derived from BaseRenderer. This also eliminates
 * the dependency on gdal in baserenderer.py.
 *
 * Revision 1.7  2005/02/16 21:14:47  jonathan
 * Further wxPython 2.5 changes using patches from Daniel Calvelo Aros
 * so that that wxproj doesn't crash. Added GUI support for selecting
 * alpha channel (opacity can't be selected yet).
 *
 * Revision 1.6  2005/02/07 19:51:13  jonathan
 * Removed unnecessary/unused macros.
 *
 * Revision 1.5  2005/01/28 19:10:07  jonathan
 * Recoded how the mask is packed into the bit array.
 *
 * Revision 1.4  2005/01/28 15:54:00  jonathan
 * Make layer's use_mask flag default to true. Support a bit array describing
 * the mask to use. Improve error handling in ProjectRasterFile (also addresses
 * RT #2947).
 *
 * Revision 1.3  2005/01/27 14:17:01  jonathan
 * Replace the old gdalwarp.cpp code with the non-simple version supplied with
 * gdal. This allows added features such as creating an alpha band.
 *
 * Revision 1.13  2004/11/14 04:57:04  fwarmerdam
 * added -srcalpha switch, and automatic alpha detection
 *
 * Revision 1.12  2004/11/05 06:15:08  fwarmerdam
 * Don't double free the warpoptions array.
 *
 * Revision 1.11  2004/11/05 05:53:43  fwarmerdam
 * Avoid various memory leaks.
 *
 * Revision 1.10  2004/10/07 15:53:42  fwarmerdam
 * added preliminary alpha band support
 *
 * Revision 1.9  2004/08/31 19:58:57  warmerda
 * Added error check if dst srs given but no source srs available.
 * http://208.24.120.44/show_bug.cgi?id=603
 *
 * Revision 1.8  2004/08/11 21:10:29  warmerda
 * Removed extra dumpopendatasets call.
 *
 * Revision 1.7  2004/08/11 20:11:24  warmerda
 * Added special VRT mode
 *
 * Revision 1.6  2004/07/28 17:56:00  warmerda
 * use return instead of exit() to avoid lame warnings on windows
 *
 * Revision 1.5  2004/04/02 17:33:22  warmerda
 * added GDALGeneralCmdLineProcessor()
 *
 * Revision 1.4  2004/04/01 19:51:18  warmerda
 * Added the -dstnodata commandline switch.
 *
 * Revision 1.3  2004/03/17 05:49:26  warmerda
 * Fixed assert check in GDALWarpCreateOutput().
 *
 * Revision 1.2  2003/09/19 17:52:21  warmerda
 * removed planned -gcp option
 *
 * Revision 1.1  2003/09/19 17:40:46  warmerda
 * Renamed from gdalwarptest.cpp to gdalwarp.cpp, replacing old gdalwarp. 
 *
 * Revision 1.12  2003/07/04 11:53:14  dron
 * Added `-rcs' option to select bicubic B-spline resampling.
 *
 * Revision 1.11  2003/06/05 16:55:42  warmerda
 * enable INIT_DEST=0 by default when creating new files
 *
 * Revision 1.10  2003/05/28 18:18:43  warmerda
 * added -q (quiet) flag
 *
 * Revision 1.9  2003/05/21 14:40:40  warmerda
 * fix error message
 *
 * Revision 1.8  2003/05/20 18:35:38  warmerda
 * added error reporting if SRS import fails
 *
 * Revision 1.7  2003/05/06 18:11:29  warmerda
 * added -multi in usage
 *
 * Revision 1.6  2003/04/23 05:18:02  warmerda
 * added -multi switch
 *
 * Revision 1.5  2003/04/21 17:21:04  warmerda
 * Fixed -wm switch.
 *
 * Revision 1.4  2003/03/28 17:43:04  warmerda
 * added -wm option to control warp memory
 *
 * Revision 1.3  2003/03/18 17:37:44  warmerda
 * add color table copying
 *
 * Revision 1.2  2003/03/02 05:24:02  warmerda
 * added -srcnodata option
 *
 * Revision 1.1  2003/02/22 02:03:41  warmerda
 * New
 *
 */

#include <Python.h>

#include "gdal.h"
#include "gdal_alg.h"
#include "gdal_priv.h"
#include "gdalwarper.h"
#include "cpl_string.h"
#include "ogr_srs_api.h"

#define PYTHON_CPL_ERR(x) \
    {const char *str = CPLGetLastErrorMsg(); \
        str != NULL ? PyErr_SetString(x, str) : PyErr_SetString(x, "");}

#define LEAVE_NOW(e) { err = e; goto getOut; }

#define OPTS_MASK  1
#define OPTS_ALPHA 2
#define OPTS_INVERT_MASK_BITS 4

CPL_CVSID("$Id: gdalwarp.cpp,v 1.9 2005/02/18 21:22:53 jonathan Exp $");

static GDALDatasetH 
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename, 
                      const char *pszFormat, const char *pszSourceSRS, 
                      const char *pszTargetSRS, int nOrder, 
                      char **papszCreateOptions, GDALDataType eDT );
static PyObject* get_gdal_version(PyObject *, PyObject * args);
static PyObject* ProjectRasterFile(PyObject *, PyObject * args);

static double	       dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
static double	       dfXRes=0.0, dfYRes=0.0;
static int             nForcePixels=0, nForceLines=0;
static int             bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
static int             bMakeMask, bMakeAlpha, bInvertMask;
const char             *pszSrcFilename = NULL, *pszDstFilename = "MEM:::";

static PyMethodDef gdalwarp_methods[] = {
    { "ProjectRasterFile", ProjectRasterFile, METH_VARARGS },
    { "get_gdal_version", get_gdal_version, METH_VARARGS },
    {NULL, NULL}
};

/************************************************************************/
/*                          get_gdal_version                            */
/************************************************************************/
static PyObject*
get_gdal_version(PyObject *, PyObject * args)
{
    PyObject *version = PyString_FromString(GDALVersionInfo("RELEASE_NAME"));

    Py_XINCREF(version);

    return version;
}

/************************************************************************/
/*                            initgdalwarp                              */
/************************************************************************/
#ifdef __cplusplus
extern "C"
#endif
void initgdalwarp(void)
{
    PyObject * shapelib = NULL;
    PyObject * c_api_func = NULL;
    PyObject * cobj = NULL;
                                                                                
    Py_InitModule("gdalwarp", gdalwarp_methods);
                                                                                
    Py_XDECREF(cobj);
    Py_XDECREF(c_api_func);
    Py_XDECREF(shapelib);
}

/************************************************************************/
/*                             SanitizeSRS                              */
/************************************************************************/

char *SanitizeSRS( const char *pszUserInput )

{
    OGRSpatialReferenceH hSRS;
    char *pszResult = NULL;

    CPLErrorReset();

    hSRS = OSRNewSpatialReference( NULL );
    if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
        OSRExportToWkt( hSRS, &pszResult );
#if 0
    else
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Translating source or target SRS failed:\n%s",
                  pszUserInput );
        exit( 1 );
    }
#endif
    
    OSRDestroySpatialReference( hSRS );

    return pszResult;
}

/************************************************************************/
/*                             GetImageData                             */
/*                                                                      */ 
/* Extract the image data from the GDALDataset and pack it into a single*/
/* array of RGB triples. Use the ColorTable to determine the RGB        */
/* values. We extract the data from the GDALDataset rather than create  */
/* our own driver because the data needs to be translated from          */
/* 4 byte pixel information into 3 byte RGB information. This could be  */
/* done as the data is written to the data set or afterwards, as it is  */
/* done here. Any minor savings from our own driver are outweighed by   */
/* the high development/maintenance costs.                              */
/*                                                                      */ 
/************************************************************************/

static CPLErr GetImageData(GDALDataset *ds, 
                           unsigned char **imgbuf, 
                           unsigned int   *imglen,
                           unsigned char **maskbuf, 
                           unsigned int   *masklen) 
{
    CPLErr ret = CE_None;

    GDALColorTable *pal = NULL;


    ds->FlushCache();

    int rasterCount  = ds->GetRasterCount();
    int nRasterXSize = ds->GetRasterXSize();
    int nRasterYSize = ds->GetRasterYSize();
    if ( ! (nRasterXSize > 0 && nRasterYSize > 0 ))
    {
        PyErr_Format(PyExc_ValueError, 
                "The dimensions (%ix%i) are invalid in %s", 
                nRasterXSize, nRasterYSize, pszSrcFilename);
        return CE_Failure;
    }

    //
    // create the new image array for RGBRGB... values
    //
    *imglen = 3 * nRasterXSize * nRasterYSize;
    *imgbuf = (unsigned char*)CPLMalloc(*imglen);
    if ( *imgbuf == NULL ) 
    {
        PyErr_Format(PyExc_MemoryError, 
                "The system does not have enough memory to project %s", 
                pszSrcFilename);
        return CE_Failure;
    }

    //
    // if there are three or more bands assume that the first three
    // are for RGB, unless told otherwise 
    //
    if (rasterCount >= 3)
    {
        for (int i=1; i <= 3; i++)
        {
            int offs = 0;
            GDALRasterBand *band = ds->GetRasterBand(i);

            switch (band->GetColorInterpretation())
            {
                case GCI_Undefined: offs = i-1; break;
                case GCI_RedBand:   offs = 1; break;
                case GCI_GreenBand: offs = 2; break;
                case GCI_BlueBand:  offs = 3; break;
                default:            offs = -1; break;
            }

            //
            // copy the image into the buffer using the proper offset
            // so we first copy over all Red values, then all Green
            // values, and then all Blue values
            //

            if (0 <= offs && offs < 3)
            {
                ret = band->RasterIO(GF_Read, 0, 0, 
                                     nRasterXSize, nRasterYSize,
                                     *imgbuf+offs, nRasterXSize, nRasterYSize, 
                                     GDT_Byte, 3, 0);
                if (ret == CE_Failure)
                {
                    PyErr_Format(PyExc_IOError, 
                        "An unknown error occured while reading band %i in %s",
                        i, pszSrcFilename);
                    break;
                }
            }
        }
    }
    else if (rasterCount >= 1)
    {
        //
        // one band is either a palette based image, or greyscale
        //

        GDALRasterBand *band = ds->GetRasterBand(1);

        switch (band->GetColorInterpretation())
        {
            case GCI_PaletteIndex:

                pal = band->GetColorTable();
                
                if (pal == NULL)
                {
                    PyErr_Format(PyExc_IOError, 
                        "Couldn't find a palette for palette-based image %s", 
                        pszSrcFilename);
                    ret = CE_Failure;
                }
                else
                {
                    GDALPaletteInterp pal_interp 
                        = pal->GetPaletteInterpretation();

                    //
                    // copy over all the palette indices and then
                    // loop through the buffer replacing the values
                    // with the correct RGB triples. 
                    //
                    ret = band->RasterIO(GF_Read, 0, 0, 
                                         nRasterXSize, nRasterYSize,
                                         *imgbuf, nRasterXSize, nRasterYSize, 
                                         GDT_UInt16, 3, 0);
                            
                    if (ret == CE_Failure) 
                    {
                        PyErr_Format(PyExc_IOError, 
                            "An unknown error occured while reading band 1 in %s", pszSrcFilename);
                        break;
                    }

                    for (unsigned char *data = *imgbuf;
                         data != (*imgbuf+*imglen);
                         data += 3)
                    {

                        unsigned short int val = *((unsigned short int *)data);

                        const GDALColorEntry *color = pal->GetColorEntry(val);

                        if (pal_interp == GPI_Gray)
                        {
                            *(data + 0) = color->c1;
                            *(data + 1) = color->c1;
                            *(data + 2) = color->c1;
                        }
                        else
                        {
                            *(data + 0) = color->c1;
                            *(data + 1) = color->c2;
                            *(data + 2) = color->c3;
                        }
                    }
                }
                break;

            case GCI_Undefined: // can we try to make a greyscale image?
            case GCI_GrayIndex:

                //
                // copy over all the palette indices and then
                // loop through the buffer replacing the values
                // with the correct RGB triples. 
                //
                ret = band->RasterIO(GF_Read, 0, 0, 
                                     nRasterXSize, nRasterYSize,
                                     *imgbuf, nRasterXSize, nRasterYSize, 
                                     GDT_Byte, 3, 0);
                        
                if (ret == CE_Failure) 
                {
                    PyErr_Format(PyExc_IOError, 
                        "An unknown error occured while reading band 1 in %s", 
                        pszSrcFilename);
                    break;
                }

                for (unsigned char *data = *imgbuf;
                     data != (*imgbuf+*imglen);
                     data += 3)
                {
                    //pal->GetColorEntry(*data, &color);

                    //*(data + 0) = *data; // already correct
                    *(data + 1) = *data;
                    *(data + 2) = *data;
                }
                break;

            default:
                PyErr_Format(PyExc_ValueError, 
                    "Unsupported color interpretation '%s' in image %s", 
                    GDALGetColorInterpretationName(
                        band->GetColorInterpretation()),
                    pszSrcFilename);

                ret = CE_Failure;
                break;
        }
    }
    else
    {
        PyErr_Format(PyExc_ValueError, 
                  "Unsupported number of raster bands (%i) in image %s\n",
                  rasterCount, pszSrcFilename);

        ret = CE_Failure;
    }

    if (ret == CE_None && bEnableDstAlpha && rasterCount > 1)
    {
        if (bMakeMask)
        {
            //
            // The mask is really an XBM image. In other words, each
            // pixel is represented by one bit in a byte array.
            //
            // First read the alpha band, and then convert it to
            // a bit array by thresholding each pixel value at 128.
            //
            
            *masklen = ((nRasterXSize + 7) / 8) * nRasterYSize;
            *maskbuf = (unsigned char *)CPLMalloc(*masklen);

            if ( *maskbuf != NULL )
            {
                unsigned char *tmp 
                    = (unsigned char *)CPLMalloc(nRasterXSize * nRasterYSize);

                if ( tmp == NULL )
                {
                    CPLFree(*maskbuf);
                    *maskbuf = NULL;
                }
                else
                {
                    GDALRasterBand *band = ds->GetRasterBand(rasterCount);

                    ret = band->RasterIO(GF_Read, 0, 0, 
                                         nRasterXSize, nRasterYSize,
                                         tmp, nRasterXSize, nRasterYSize, 
                                         GDT_Byte, 0, 0);

                    if (ret != CE_Failure)
                    {
                        int i, j, b=1, c=0;
                        unsigned char *ptr = *maskbuf;
                        unsigned char *tptr = tmp;

                        //unsigned int empty_count=0;

                        for (i=0; i < nRasterYSize; i++)
                        {
                            for (j=nRasterXSize; j >= 8; j -= 8)
                            {
                                c=0; b=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                if (*tptr++ >= 128) {c|=b;} b<<=1;
                                
                                if (bInvertMask) 
                                    *(ptr++) = ~c;
                                else
                                    *(ptr++) = c;
                            }

                            c=0; b=1;
                            switch (nRasterXSize & 7)
                            {
                                case 7: if (*tptr++ >= 128) {c|=b;} b<<=1;
                                case 6: if (*tptr++ >= 128) {c|=b;} b<<=1;
                                case 5: if (*tptr++ >= 128) {c|=b;} b<<=1;
                                case 4: if (*tptr++ >= 128) {c|=b;} b<<=1;
                                case 3: if (*tptr++ >= 128) {c|=b;} b<<=1;
                                case 2: if (*tptr++ >= 128) {c|=b;} b<<=1;
                                case 1: 
                                    if (*tptr++ >= 128) {c|=b;} 
                                    b<<=1;

                                    //
                                    // byte should be padded with 0's so
                                    // it's not a simple inversion
                                    //
                                    if (bInvertMask) 
                                        *(ptr++) = ~c & (b-1);
                                    else
                                        *(ptr++) = c;

                                default: break;
                            }
                        }

#if 0
                        if (empty_count == *masklen)
                        {
                            fprintf(stderr, "mask not used\n");

                            CPLFree(*maskbuf);
                            *maskbuf = NULL;
                        }
#endif
                    }

                    CPLFree(tmp);
                    tmp = NULL;
                }

            }
        }
        else if (bMakeAlpha)
        {
            //
            // This is the simple case. The array we get back from RasterIO
            // is already in the correct format.
            //
            
            *masklen = nRasterXSize * nRasterYSize;
            *maskbuf = (unsigned char *)CPLMalloc(*masklen);

            if ( *maskbuf != NULL )
            {
                GDALRasterBand *band = ds->GetRasterBand(rasterCount);

                ret = band->RasterIO(GF_Read, 0, 0, 
                                     nRasterXSize, nRasterYSize,
                                     *maskbuf, nRasterXSize, nRasterYSize, 
                                     GDT_Byte, 0, 0);

#if 0
                if (ret == CE_Failure)
                {
                    CPLFree(*maskbuf);
                    *maskbuf = NULL;
                }
#endif
            }
        }
    }

    if (ret != CE_None) 
    {
        if (*imgbuf  != NULL) { CPLFree(*imgbuf);  *imgbuf  = NULL; }
        if (*maskbuf != NULL) { CPLFree(*maskbuf); *maskbuf = NULL; }
    }

    return ret;
}

/************************************************************************/
/*                          ProjectRasterFile                           */
/************************************************************************/

static PyObject*
ProjectRasterFile(PyObject *, PyObject * args)
{
    GDALDatasetH	hSrcDS=NULL, hDstDS=NULL;
    const char         *pszFormat = "MEM";
    char               *pszTargetSRS = NULL;
    char               *pszSourceSRS = NULL;
    int                 bCreateOutput = FALSE, i, nOrder = 0;
    void               *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL;
    char               **papszWarpOptions = NULL;
    double             dfErrorThreshold = 0.125;
    double             dfWarpMemoryLimit = 0.0;
    GDALTransformerFunc pfnTransformer = NULL;
    char                **papszCreateOptions = NULL;
    GDALDataType        eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown; 
    GDALResampleAlg     eResampleAlg = GRA_NearestNeighbour;
    int                 bMulti = FALSE;
    int                 err = 0;
    unsigned char      *imgbuf = NULL;
    unsigned int        imglen = 0;
    unsigned char      *maskbuf = NULL;
    unsigned int        masklen = 0;
    int                 options = 0;

    GDALWarpOperation oWO;
    GDALWarpOptions *psWO = NULL;

    PyObject * pyImageData = NULL;
    PyObject * pyMaskData = NULL;
    PyObject * filename;
    PyObject * srcImageArgs;
    PyObject * dstImageArgs;
    PyObject * extents;
    PyObject * resolution;
    PyObject * imageRes;
    PyObject * opts;
    PyObject * pyReturnData = NULL;


    if (!PyArg_ParseTuple(args, "OOOOOOO", &filename, &srcImageArgs, 
                &dstImageArgs, 
                &extents, 
                &resolution, &imageRes,
                &opts))
    {
        return NULL;
    }

    dfXRes=0.0; dfYRes=0.0;

    pszSrcFilename = PyString_AsString( filename );

    pszSourceSRS = SanitizeSRS( PyString_AsString( srcImageArgs ) );
    pszTargetSRS = SanitizeSRS( PyString_AsString( dstImageArgs ) );

    dfMinX = PyFloat_AsDouble( PyTuple_GetItem( extents, 0 ) );
    dfMinY = PyFloat_AsDouble( PyTuple_GetItem( extents, 1 ) );
    dfMaxX = PyFloat_AsDouble( PyTuple_GetItem( extents, 2 ) );
    dfMaxY = PyFloat_AsDouble( PyTuple_GetItem( extents, 3 ) );

    nForcePixels = ( int )PyInt_AsLong( PyTuple_GetItem( imageRes, 0 ) );
    nForceLines  = ( int )PyInt_AsLong( PyTuple_GetItem( imageRes, 1 ) );

    options  = ( int )PyInt_AsLong( opts );
    bMakeMask  = (options & OPTS_MASK) == OPTS_MASK;
    bMakeAlpha = (options & OPTS_ALPHA) == OPTS_ALPHA;
    bInvertMask = (options & OPTS_INVERT_MASK_BITS) == OPTS_INVERT_MASK_BITS;

    // FIXME: error if bMakeMask == bMaskAlpha

    bEnableDstAlpha = bMakeMask || bMakeAlpha;

    GDALAllRegister();

/* -------------------------------------------------------------------- */
/*      Open source dataset.                                            */
/* -------------------------------------------------------------------- */
    hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
    
    if( hSrcDS == NULL )
    {
        PYTHON_CPL_ERR( PyExc_IOError );
        LEAVE_NOW( CPLGetLastErrorNo() );
    }

    if( pszSourceSRS == NULL )
    {
        if( GDALGetProjectionRef( hSrcDS ) != NULL 
            && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
            pszSourceSRS = CPLStrdup(GDALGetProjectionRef( hSrcDS ));

        else if( GDALGetGCPProjection( hSrcDS ) != NULL
                 && strlen(GDALGetGCPProjection(hSrcDS)) > 0 
                 && GDALGetGCPCount( hSrcDS ) > 1 )
            pszSourceSRS = CPLStrdup(GDALGetGCPProjection( hSrcDS ));
        else
            pszSourceSRS = CPLStrdup("");
    }

    if( pszTargetSRS != NULL && strlen(pszSourceSRS) == 0 )
    {
        PyErr_Format(PyExc_ValueError, 
            "A target projection was specified, "
            "but there is no source projection in %s", pszSrcFilename );
        LEAVE_NOW( 1 );
    }

    if( pszTargetSRS == NULL )
        pszTargetSRS = CPLStrdup(pszSourceSRS);

    if( GDALGetRasterColorInterpretation( 
            GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) ) 
        == GCI_AlphaBand 
        && !bEnableSrcAlpha )
    {
        bEnableSrcAlpha = TRUE;
#if 0
        printf( "Using band %d of source image as alpha.\n", 
                GDALGetRasterCount(hSrcDS) );
#endif
    }
        

/* -------------------------------------------------------------------- */
/*      If not, we need to create it.                                   */
/* -------------------------------------------------------------------- */
    if( hDstDS == NULL )
    {
        hDstDS = GDALWarpCreateOutput( hSrcDS, pszDstFilename, pszFormat, 
                                       pszSourceSRS, pszTargetSRS, nOrder,
                                       papszCreateOptions, eOutputType );
        bCreateOutput = TRUE;

        papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");

        CSLDestroy( papszCreateOptions );
        papszCreateOptions = NULL;
    }

    if( hDstDS == NULL )
    {
        PyErr_Format(PyExc_IOError, 
            "Error creating destination image for projecting %s",
            pszSrcFilename);
        LEAVE_NOW( CPLE_FileIO );
    }

/* -------------------------------------------------------------------- */
/*      Create a transformation object from the source to               */
/*      destination coordinate system.                                  */
/* -------------------------------------------------------------------- */
    hTransformArg = hGenImgProjArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS, 
                                         hDstDS, pszTargetSRS, 
                                         TRUE, 1000.0, nOrder );


    if( hTransformArg == NULL )
    {
        PYTHON_CPL_ERR(PyExc_ValueError );
        LEAVE_NOW( CPLE_IllegalArg );
    }

    pfnTransformer = GDALGenImgProjTransform;

    CPLFree( pszSourceSRS );
    pszSourceSRS = NULL;

    CPLFree( pszTargetSRS );
    pszTargetSRS = NULL;

/* -------------------------------------------------------------------- */
/*      Warp the transformer with a linear approximator unless the      */
/*      acceptable error is zero.                                       */
/* -------------------------------------------------------------------- */
    if( dfErrorThreshold != 0.0 )
    {
        hTransformArg = hApproxArg = 
            GDALCreateApproxTransformer( GDALGenImgProjTransform, 
                                         hGenImgProjArg, dfErrorThreshold );
        pfnTransformer = GDALApproxTransform;
    }

/* -------------------------------------------------------------------- */
/*      Setup warp options.                                             */
/* -------------------------------------------------------------------- */
    psWO = GDALCreateWarpOptions();

    psWO->papszWarpOptions = papszWarpOptions;
    psWO->eWorkingDataType = eWorkingType;
    psWO->eResampleAlg = eResampleAlg;

    psWO->hSrcDS = hSrcDS;
    psWO->hDstDS = hDstDS;

    psWO->pfnTransformer = pfnTransformer;
    psWO->pTransformerArg = hTransformArg;

    if( dfWarpMemoryLimit != 0.0 )
        psWO->dfWarpMemoryLimit = dfWarpMemoryLimit;

/* -------------------------------------------------------------------- */
/*      Setup band mapping.                                             */
/* -------------------------------------------------------------------- */
    if( bEnableSrcAlpha )
        psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
    else
        psWO->nBandCount = GDALGetRasterCount(hSrcDS);

    psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
    psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));

    for( i = 0; i < psWO->nBandCount; i++ )
    {
        psWO->panSrcBands[i] = i+1;
        psWO->panDstBands[i] = i+1;
    }

/* -------------------------------------------------------------------- */
/*      Setup alpha bands used if any.                                  */
/* -------------------------------------------------------------------- */
    if( bEnableSrcAlpha )
        psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);

    if( !bEnableDstAlpha 
        && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1 
        && GDALGetRasterColorInterpretation( 
            GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS))) 
        == GCI_AlphaBand )
    {
#if 0
        printf( "Using band %d of destination image as alpha.\n", 
                GDALGetRasterCount(hDstDS) );
#endif
                
        bEnableDstAlpha = TRUE;
    }

    if( bEnableDstAlpha )
        psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);

/* -------------------------------------------------------------------- */
/*      Setup NODATA options.                                           */
/* -------------------------------------------------------------------- */
#if 0
    if( pszSrcNodata != NULL )
    {
        char **papszTokens = CSLTokenizeString( pszSrcNodata );
        int  nTokenCount = CSLCount(papszTokens);

        psWO->padfSrcNoDataReal = (double *) 
            CPLMalloc(psWO->nBandCount*sizeof(double));
        psWO->padfSrcNoDataImag = (double *) 
            CPLMalloc(psWO->nBandCount*sizeof(double));

        for( i = 0; i < psWO->nBandCount; i++ )
        {
            if( i < nTokenCount )
            {
                CPLStringToComplex( papszTokens[i], 
                                    psWO->padfSrcNoDataReal + i,
                                    psWO->padfSrcNoDataImag + i );
            }
            else
            {
                psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i-1];
                psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i-1];
            }
        }

        CSLDestroy( papszTokens );
    }

/* -------------------------------------------------------------------- */
/*      If the output dataset was created, and we have a destination    */
/*      nodata value, go through marking the bands with the information.*/
/* -------------------------------------------------------------------- */
    if( pszDstNodata != NULL && bCreateOutput )
    {
        char **papszTokens = CSLTokenizeString( pszDstNodata );
        int  nTokenCount = CSLCount(papszTokens);

        psWO->padfDstNoDataReal = (double *) 
            CPLMalloc(psWO->nBandCount*sizeof(double));
        psWO->padfDstNoDataImag = (double *) 
            CPLMalloc(psWO->nBandCount*sizeof(double));

        for( i = 0; i < psWO->nBandCount; i++ )
        {
            if( i < nTokenCount )
            {
                CPLStringToComplex( papszTokens[i], 
                                    psWO->padfDstNoDataReal + i,
                                    psWO->padfDstNoDataImag + i );
            }
            else
            {
                psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i-1];
                psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i-1];
            }

            if( bCreateOutput )
            {
                GDALSetRasterNoDataValue( 
                    GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ), 
                    psWO->padfDstNoDataReal[i] );
            }
        }

        CSLDestroy( papszTokens );
    }
#endif

/* -------------------------------------------------------------------- */
/*      Initialize and execute the warp.                                */
/* -------------------------------------------------------------------- */

    if( oWO.Initialize( psWO ) == CE_None )
    {
        if( bMulti )
            oWO.ChunkAndWarpMulti( 0, 0, 
                                   GDALGetRasterXSize( hDstDS ),
                                   GDALGetRasterYSize( hDstDS ) );
        else
            oWO.ChunkAndWarpImage( 0, 0, 
                                   GDALGetRasterXSize( hDstDS ),
                                   GDALGetRasterYSize( hDstDS ) );
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
/* -------------------------------------------------------------------- */
    if( hApproxArg != NULL )
        GDALDestroyApproxTransformer( hApproxArg );

    if( hGenImgProjArg != NULL )
        GDALDestroyGenImgProjTransformer( hGenImgProjArg );

/* -------------------------------------------------------------------- */
/*      Cleanup.                                                        */
/* -------------------------------------------------------------------- */
    CPLErrorReset();
    err = 0;

    imglen = 0;
    imgbuf = NULL;
    masklen = 0;
    maskbuf = NULL;

    if ( GetImageData((GDALDataset *)hDstDS, 
                      &imgbuf, &imglen, &maskbuf, &masklen) == CE_None)
    {
        pyImageData = PyString_FromStringAndSize( ( char * )imgbuf,  imglen);

        pyReturnData = PyTuple_New(3);
        PyTuple_SetItem(pyReturnData, 0, pyImageData);

        if (bMakeAlpha && maskbuf != NULL)
        {
            pyMaskData = PyString_FromStringAndSize( ( char * )maskbuf,masklen);
            PyTuple_SetItem(pyReturnData, 1, Py_None);
            PyTuple_SetItem(pyReturnData, 2, pyMaskData);
        }
        else if (bMakeMask && maskbuf != NULL)
        {
            pyMaskData = PyString_FromStringAndSize( ( char * )maskbuf,masklen);
            PyTuple_SetItem(pyReturnData, 1, pyMaskData);
            PyTuple_SetItem(pyReturnData, 2, Py_None);
        }
        else
        {
            PyTuple_SetItem(pyReturnData, 1, Py_None);
            PyTuple_SetItem(pyReturnData, 2, Py_None);
        }

        if (imgbuf  != NULL) { CPLFree(imgbuf);  imgbuf  = NULL; }
        if (maskbuf != NULL) { CPLFree(maskbuf); maskbuf = NULL; }

    }

getOut:

    GDALClose( hDstDS );
    GDALClose( hSrcDS );

    if (psWO != NULL) GDALDestroyWarpOptions( psWO );

    //GDALDumpOpenDatasets( stderr );

    GDALDestroyDriverManager();
    
    if ( !err && CPLGetLastErrorNo() )
    {
        PYTHON_CPL_ERR( PyExc_StandardError );
        return NULL;
    }

    if ( err ) 
    {
        return NULL;
    }
    else
    {
        return pyReturnData;
    }
}

/************************************************************************/
/*                        GDALWarpCreateOutput()                        */
/*                                                                      */
/*      Create the output file based on various commandline options,    */
/*      and the input file.                                             */
/************************************************************************/

static GDALDatasetH 
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename, 
                      const char *pszFormat, const char *pszSourceSRS, 
                      const char *pszTargetSRS, int nOrder,
                      char **papszCreateOptions, GDALDataType eDT )


{
    GDALDriverH hDriver;
    GDALDatasetH hDstDS;
    void *hTransformArg;
    double adfDstGeoTransform[6];
    int nPixels=0, nLines=0;

    CPLErrorReset();

    if( eDT == GDT_Unknown )
        eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));

/* -------------------------------------------------------------------- */
/*      Find the output driver.                                         */
/* -------------------------------------------------------------------- */
    hDriver = GDALGetDriverByName( pszFormat );
    if( hDriver == NULL 
        || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
    {
#if 0
        int	iDr;
        
        printf( "Output driver `%s' not recognised or does not support\n", 
                pszFormat );
        printf( "direct output file creation.  The following format drivers are configured\n"
                "and support direct output:\n" );

        for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
        {
            GDALDriverH hDriver = GDALGetDriver(iDr);

            if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
            {
                printf( "  %s: %s\n",
                        GDALGetDriverShortName( hDriver  ),
                        GDALGetDriverLongName( hDriver ) );
            }
        }
        printf( "\n" );
#endif

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Create a transformation object from the source to               */
/*      destination coordinate system.                                  */
/* -------------------------------------------------------------------- */
    hTransformArg = 
        GDALCreateGenImgProjTransformer( hSrcDS, pszSourceSRS, 
                                         NULL, pszTargetSRS, 
                                         TRUE, 1000.0, nOrder );

    if( hTransformArg == NULL )
        return NULL;

    //
    // This could happen if the proj library didn't load correctly
    // Fixes RF#2947
    //
    if (CPLGetLastErrorNo() != CE_None)
    {
        GDALDestroyGenImgProjTransformer(hTransformArg);
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Get approximate output definition.                              */
/* -------------------------------------------------------------------- */
    if( GDALSuggestedWarpOutput( hSrcDS, 
                                 GDALGenImgProjTransform, hTransformArg, 
                                 adfDstGeoTransform, &nPixels, &nLines )
        != CE_None )
        return NULL;

    GDALDestroyGenImgProjTransformer( hTransformArg );

/* -------------------------------------------------------------------- */
/*      Did the user override some parameters?                          */
/* -------------------------------------------------------------------- */
    if( dfXRes != 0.0 && dfYRes != 0.0 )
    {
        CPLAssert( nForcePixels == 0 && nForceLines == 0 );
        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
        {
            dfMinX = adfDstGeoTransform[0];
            dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
            dfMaxY = adfDstGeoTransform[3];
            dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
        }

        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
        adfDstGeoTransform[0] = dfMinX;
        adfDstGeoTransform[3] = dfMaxY;
        adfDstGeoTransform[1] = dfXRes;
        adfDstGeoTransform[5] = -dfYRes;
    }

    else if( nForcePixels != 0 && nForceLines != 0 )
    {
        if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
        {
            dfMinX = adfDstGeoTransform[0];
            dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
            dfMaxY = adfDstGeoTransform[3];
            dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
        }

        dfXRes = (dfMaxX - dfMinX) / nForcePixels;
        dfYRes = (dfMaxY - dfMinY) / nForceLines;

        adfDstGeoTransform[0] = dfMinX;
        adfDstGeoTransform[3] = dfMaxY;
        adfDstGeoTransform[1] = dfXRes;
        adfDstGeoTransform[5] = -dfYRes;

        nPixels = nForcePixels;
        nLines = nForceLines;
    }

    else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
    {
        dfXRes = adfDstGeoTransform[1];
        dfYRes = fabs(adfDstGeoTransform[5]);

        nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
        nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);

        adfDstGeoTransform[0] = dfMinX;
        adfDstGeoTransform[3] = dfMaxY;
    }

/* -------------------------------------------------------------------- */
/*      Do we want to generate an alpha band in the output file?        */
/* -------------------------------------------------------------------- */
    int nDstBandCount = GDALGetRasterCount(hSrcDS);

    if( bEnableSrcAlpha )
        nDstBandCount--;

    if( bEnableDstAlpha )
        nDstBandCount++;

/* -------------------------------------------------------------------- */
/*      Create the output file.                                         */
/* -------------------------------------------------------------------- */
    hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines, 
                         nDstBandCount, eDT, papszCreateOptions );
    
    if( hDstDS == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Write out the projection definition.                            */
/* -------------------------------------------------------------------- */
    GDALSetProjection( hDstDS, pszTargetSRS );
    GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

/* -------------------------------------------------------------------- */
/*      Try to set color interpretation of output file alpha band.      */
/*      TODO: We should likely try to copy the other bands too.         */
/* -------------------------------------------------------------------- */
    if( bEnableDstAlpha )
    {
        GDALSetRasterColorInterpretation( 
            GDALGetRasterBand( hDstDS, nDstBandCount ), 
            GCI_AlphaBand );
    }

/* -------------------------------------------------------------------- */
/*      Copy the color table, if required.                              */
/* -------------------------------------------------------------------- */
    GDALColorTableH hCT;

    hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
    if( hCT != NULL )
        GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );

    return hDstDS;
}
    
