jonathan: thuban/libraries/thuban gdalwarp.cpp,1.2,1.3
cvs@intevation.de
cvs at intevation.de
Thu Jan 27 15:17:03 CET 2005
Author: jonathan
Update of /thubanrepository/thuban/libraries/thuban
In directory doto:/tmp/cvs-serv16816
Modified Files:
gdalwarp.cpp
Log Message:
Replace the old gdalwarp.cpp code with the non-simple version supplied with
gdal. This allows added features such as creating an alpha band.
Index: gdalwarp.cpp
===================================================================
RCS file: /thubanrepository/thuban/libraries/thuban/gdalwarp.cpp,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -d -r1.2 -r1.3
--- gdalwarp.cpp 21 Jan 2005 14:01:25 -0000 1.2
+++ gdalwarp.cpp 27 Jan 2005 14:17:01 -0000 1.3
@@ -1,9 +1,8 @@
/******************************************************************************
* $Id$
*
- * Project: Mapinfo Image Warper
- * Purpose: Commandline program for doing a variety of image warps, including
- * image reprojection.
+ * Project: High Performance Image Reprojector
+ * Purpose: Test program for high performance warper API.
* Author: Frank Warmerdam <warmerdam at pobox.com>
*
******************************************************************************
@@ -30,67 +29,84 @@
******************************************************************************
*
* $Log$
- * Revision 1.2 2005/01/21 14:01:25 jonathan
- * Improved rendering raster layers by changing the return format of
- * the ProjectRasterFile function.
+ * 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.1 2003/08/19 21:32:24 jan
- * These files have been moved here from thuban/extensions/thuban/
- * See there in the Attic for the older history.
+ * Revision 1.13 2004/11/14 04:57:04 fwarmerdam
+ * added -srcalpha switch, and automatic alpha detection
*
- * Revision 1.7 2003/07/10 14:55:11 jonathan
- * (ProjectRasterFile): Make sure
- * the function doesn't return NULL without first setting a Python
- * Error.
+ * Revision 1.12 2004/11/05 06:15:08 fwarmerdam
+ * Don't double free the warpoptions array.
*
- * Revision 1.6 2003/06/26 17:00:59 jonathan
- * (get_gdal_version): New.
- * Returns the version of gdal library being used as a string.
+ * Revision 1.11 2004/11/05 05:53:43 fwarmerdam
+ * Avoid various memory leaks.
*
- * Revision 1.5 2003/06/17 15:25:56 jonathan
- * (MFILENAME): Padding should be to twice sizeof(void*) because there are
- * two digits for each hex byte.
+ * Revision 1.10 2004/10/07 15:53:42 fwarmerdam
+ * added preliminary alpha band support
*
- * Revision 1.4 2003/06/12 17:33:25 jonathan
- * Removed debug printing as the image is rendered. Fixes RTbug #1937.
+ * 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.3 2003/05/30 06:30:43 jonathan
- * Removed code that allowed gdalwarp
- * to be compiled as a standalone app. Now the code can only be
- * called from Python which simplifies the parameter passing.
- * (ProjectRasterFile): Handle Python arguments. Remove code that
- * checks for a destination dataset. Add more clean up code.
+ * Revision 1.8 2004/08/11 21:10:29 warmerda
+ * Removed extra dumpopendatasets call.
*
- * Revision 1.2 2003/05/21 17:25:16 jonathan
- * Improved the error handling
- * and passed GDAL messages back up to the Python layer. Also
- * tried to fix some memory leaks that were present in the original
- * utility but didn't matter because the program aborted.
+ * Revision 1.7 2004/08/11 20:11:24 warmerda
+ * Added special VRT mode
*
- * Revision 1.1 2003/05/20 15:26:17 jonathan
- * New, but basically a direct copy
- * of the gdalwarp utility provided in GDAL. Added function calls
- * that can be accessed from python.
+ * Revision 1.6 2004/07/28 17:56:00 warmerda
+ * use return instead of exit() to avoid lame warnings on windows
*
- * Revision 1.7 2003/03/18 17:33:00 warmerda
- * Copy colortable if there is one on the source file.
+ * Revision 1.5 2004/04/02 17:33:22 warmerda
+ * added GDALGeneralCmdLineProcessor()
*
- * Revision 1.6 2002/12/17 18:23:19 warmerda
- * copy GCP projection if the main projection isn't set meaningfully
+ * Revision 1.4 2004/04/01 19:51:18 warmerda
+ * Added the -dstnodata commandline switch.
*
- * Revision 1.5 2002/12/13 04:57:15 warmerda
- * implementing remaining commandline options
+ * Revision 1.3 2004/03/17 05:49:26 warmerda
+ * Fixed assert check in GDALWarpCreateOutput().
*
- * Revision 1.4 2002/12/09 16:08:57 warmerda
- * added approximating transformer
+ * Revision 1.2 2003/09/19 17:52:21 warmerda
+ * removed planned -gcp option
*
- * Revision 1.3 2002/12/07 22:58:09 warmerda
- * pass initialization warp option
+ * Revision 1.1 2003/09/19 17:40:46 warmerda
+ * Renamed from gdalwarptest.cpp to gdalwarp.cpp, replacing old gdalwarp.
*
- * Revision 1.2 2002/12/07 17:09:13 warmerda
- * added -order flag
+ * Revision 1.12 2003/07/04 11:53:14 dron
+ * Added `-rcs' option to select bicubic B-spline resampling.
*
- * Revision 1.1 2002/12/07 16:13:38 warmerda
+ * 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
*
*/
@@ -100,17 +116,18 @@
#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_ERR(x) \
+#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 PYTHON_ERR(n, str) \
+ {str != NULL ? PyErr_SetString(n, str) : PyErr_SetString(n, "");}
-#define NO_DATA_VAL 255
-#define NO_DATA_VAL_STR "255"
+#define LEAVE_NOW(e) { err = e; goto getOut; }
CPL_CVSID("$Id$");
@@ -118,13 +135,14 @@
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename,
const char *pszFormat, const char *pszSourceSRS,
const char *pszTargetSRS, int nOrder,
- char **papszCreateOptions );
+ 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 PyMethodDef gdalwarp_methods[] = {
{ "ProjectRasterFile", ProjectRasterFile, METH_VARARGS },
@@ -174,10 +192,21 @@
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;
@@ -203,11 +232,12 @@
/************************************************************************/
static CPLErr GetImageData(GDALDataset *ds,
- unsigned char **imgbuf, unsigned int *imglen)
+ unsigned char **imgbuf,
+ unsigned char **maskbuf,
+ unsigned int *imglen)
{
CPLErr ret = CE_None;
- GDALColorEntry color;
GDALColorTable *pal = NULL;
ds->FlushCache();
@@ -223,27 +253,32 @@
//
*imglen = 3 * nRasterXSize * nRasterYSize;
*imgbuf = new unsigned char[*imglen];
- if ( *imgbuf == NULL )
- return CE_Failure;
+ if ( *imgbuf == NULL ) { ret=CE_Failure; goto no_mem; }
- memset(*imgbuf, 170, *imglen);
+ if (bEnableDstAlpha)
+ {
+ *maskbuf = new unsigned char[*imglen];
+ if ( *maskbuf == NULL ) { ret=CE_Failure; goto no_mem; }
+ }
//
- // if there are three bands assume for now that they are RGB bands
+ // if there are three or more bands assume that the first three
+ // are for RGB, unless told otherwise
//
- if (rasterCount == 3)
+ if (rasterCount >= 3)
{
- for (int i=1; i <= rasterCount; i++)
+ 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 = i-1; break;
+ default: offs = -1; break;
}
//
@@ -251,15 +286,18 @@
// so we first copy over all Red values, then all Green
// values, and then all Blue values
//
- ret = band->RasterIO(GF_Read, 0, 0,
- nRasterXSize, nRasterYSize,
- *imgbuf+offs, nRasterXSize, nRasterYSize,
- GDT_Byte, 3, 0);
-
- if (ret == CE_Failure) break;
+
+ 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) break;
+ }
}
}
- else if (rasterCount == 1)
+ else if (rasterCount >= 1)
{
//
// one band is either a palette based image, or greyscale
@@ -275,8 +313,11 @@
// FIXME: what happens if pal is NULL? ignore?
// test this once to see if conversion is supported.
- if (pal && pal->GetColorEntryAsRGB(0, &color))
+ if (pal != NULL)
{
+ GDALPaletteInterp pal_interp
+ = pal->GetPaletteInterpretation();
+
//
// copy over all the palette indices and then
// loop through the buffer replacing the values
@@ -285,7 +326,7 @@
ret = band->RasterIO(GF_Read, 0, 0,
nRasterXSize, nRasterYSize,
*imgbuf, nRasterXSize, nRasterYSize,
- GDT_Byte, 3, 0);
+ GDT_UInt16, 3, 0);
if (ret == CE_Failure) break;
@@ -293,19 +334,22 @@
data != (*imgbuf+*imglen);
data += 3)
{
- if (*data == NO_DATA_VAL)
+
+ unsigned short int val = *((unsigned short int *)data);
+
+ const GDALColorEntry *color = pal->GetColorEntry(val);
+
+ if (pal_interp == GPI_Gray)
{
- *(data + 0) = 3;
- *(data + 1) = 1;
- *(data + 2) = 4;
+ *(data + 0) = color->c1;
+ *(data + 1) = color->c1;
+ *(data + 2) = color->c1;
}
else
{
- pal->GetColorEntryAsRGB(*data, &color);
-
- *(data + 0) = color.c1;
- *(data + 1) = color.c2;
- *(data + 2) = color.c3;
+ *(data + 0) = color->c1;
+ *(data + 1) = color->c2;
+ *(data + 2) = color->c3;
}
}
}
@@ -330,20 +374,11 @@
data != (*imgbuf+*imglen);
data += 3)
{
- if (*data == NO_DATA_VAL)
- {
- *(data + 0) = 3;
- *(data + 1) = 1;
- *(data + 2) = 4;
- }
- else
- {
- //pal->GetColorEntry(*data, &color);
+ //pal->GetColorEntry(*data, &color);
- //*(data + 0) = *data; // already correct
- *(data + 1) = *data;
- *(data + 2) = *data;
- }
+ //*(data + 0) = *data; // already correct
+ *(data + 1) = *data;
+ *(data + 2) = *data;
}
break;
@@ -355,49 +390,91 @@
break;
}
}
+ else
+ {
+ memset(*imgbuf, 170, *imglen);
+ if (bEnableDstAlpha) memset(*maskbuf, 255, *imglen);
+ }
+
+ if (bEnableDstAlpha && rasterCount > 1)
+ {
+ GDALRasterBand *band = ds->GetRasterBand(rasterCount);
+
+ //
+ // Unfortunately, wxImage requires the data be in RGB format,
+ // so here the data is tripled
+ //
+ for (int i=0; i < 3; i++)
+ {
+ ret = band->RasterIO(GF_Read, 0, 0,
+ nRasterXSize, nRasterYSize,
+ *maskbuf+i, nRasterXSize, nRasterYSize,
+ GDT_Byte, 3, 0);
+ }
+ }
+
+no_mem:
if (ret != CE_None)
- if (imgbuf != NULL) delete imgbuf;
+ {
+ if (imgbuf != NULL) { delete imgbuf; imgbuf = NULL; }
+ if (maskbuf != NULL) { delete maskbuf; maskbuf = NULL; }
+ }
return ret;
}
+/************************************************************************/
+/* ProjectRasterFile */
+/************************************************************************/
+
static PyObject*
ProjectRasterFile(PyObject *, PyObject * args)
{
- GDALDatasetH hSrcDS = NULL, hDstDS = NULL;
+ GDALDatasetH hSrcDS=NULL, hDstDS=NULL;
const char *pszFormat = "MEM";
- const char *pszDstFilename = "MEM:::";
char *pszTargetSRS = NULL;
char *pszSourceSRS = NULL;
- const char *pszSrcFilename = NULL;
- int bCreateOutput = FALSE, nOrder = 0;
+ const char *pszSrcFilename = NULL, *pszDstFilename = "MEM:::";
+ 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;
+
+ GDALWarpOperation oWO;
+ GDALWarpOptions *psWO = NULL;
PyObject * pyImageData = NULL;
+ PyObject * pyMaskData = NULL;
PyObject * filename;
PyObject * srcImageArgs;
PyObject * dstImageArgs;
PyObject * extents;
PyObject * resolution;
PyObject * imageRes;
+ PyObject * mask;
+ PyObject * pyReturnData = NULL;
- if (!PyArg_ParseTuple(args, "OOOOOO", &filename, &srcImageArgs,
+
+ if (!PyArg_ParseTuple(args, "OOOOOOO", &filename, &srcImageArgs,
&dstImageArgs,
&extents,
- &resolution, &imageRes))
+ &resolution, &imageRes,
+ &mask))
{
return NULL;
}
-
dfXRes=0.0; dfYRes=0.0;
pszSrcFilename = PyString_AsString( filename );
@@ -413,22 +490,18 @@
nForcePixels = ( int )PyInt_AsLong( PyTuple_GetItem( imageRes, 0 ) );
nForceLines = ( int )PyInt_AsLong( PyTuple_GetItem( imageRes, 1 ) );
- bCreateOutput = TRUE;
+ bEnableDstAlpha = mask == Py_True;
GDALAllRegister();
- GDALRegister_MEM();
/* -------------------------------------------------------------------- */
/* Open source dataset. */
/* -------------------------------------------------------------------- */
-
- CPLPushErrorHandler( CPLQuietErrorHandler );
-
hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
if( hSrcDS == NULL )
{
- PYTHON_ERR( PyExc_IOError );
+ PYTHON_CPL_ERR( PyExc_IOError );
LEAVE_NOW( CPLGetLastErrorNo() );
}
@@ -446,20 +519,48 @@
pszSourceSRS = CPLStrdup("");
}
+ if( pszTargetSRS != NULL && strlen(pszSourceSRS) == 0 )
+ {
+ PYTHON_ERR(PyExc_ValueError,
+ "A target projection was specified, but there is no source projection." );
+ LEAVE_NOW( 1 );
+ }
+
if( pszTargetSRS == NULL )
pszTargetSRS = CPLStrdup(pszSourceSRS);
- hDstDS = GDALWarpCreateOutput( hSrcDS, pszDstFilename, pszFormat,
- pszSourceSRS, pszTargetSRS, nOrder,
- papszCreateOptions );
- papszWarpOptions = CSLSetNameValue( papszWarpOptions,
- "INIT", NO_DATA_VAL_STR );
- CSLDestroy( papszCreateOptions );
- papszCreateOptions = NULL;
+ 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 )
{
- PYTHON_ERR( PyExc_IOError );
+ 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 )
+ {
+ PYTHON_CPL_ERR( PyExc_IOError );
LEAVE_NOW( CPLE_FileIO );
}
@@ -474,12 +575,18 @@
if( hTransformArg == NULL )
{
- PYTHON_ERR( PyExc_ValueError );
+ 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. */
@@ -492,17 +599,155 @@
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;
/* -------------------------------------------------------------------- */
-/* Now actually invoke the warper to do the work. */
+/* Setup band mapping. */
/* -------------------------------------------------------------------- */
- GDALSimpleImageWarp( hSrcDS, hDstDS, 0, NULL,
- pfnTransformer, hTransformArg,
- GDALDummyProgress, NULL, papszWarpOptions );
+ if( bEnableSrcAlpha )
+ psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
+ else
+ psWO->nBandCount = GDALGetRasterCount(hSrcDS);
- CSLDestroy( papszWarpOptions );
- papszWarpOptions = NULL;
+ 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 );
@@ -517,48 +762,57 @@
imglen = 0;
imgbuf = NULL;
+ maskbuf = NULL;
- if ( GetImageData((GDALDataset *)hDstDS, &imgbuf, &imglen) == CE_None)
+#if 1
+
+ if ( GetImageData((GDALDataset *)hDstDS, &imgbuf, &maskbuf, &imglen)
+ == CE_None)
{
- //printf("imglen: %i\n", imglen);
+ pyImageData = PyString_FromStringAndSize( ( char * )imgbuf, imglen);
- pyImageData = PyString_FromStringAndSize( ( char * )imgbuf, imglen);
+ if (bEnableDstAlpha)
+ pyMaskData = PyString_FromStringAndSize( ( char * )maskbuf, imglen);
+ else
+ pyMaskData = Py_None;
- delete imgbuf;
- imgbuf = NULL;
+ pyReturnData = PyTuple_New(2);
+ PyTuple_SetItem(pyReturnData, 0, pyImageData);
+ PyTuple_SetItem(pyReturnData, 1, pyMaskData);
+
+ if (imgbuf != NULL) { delete imgbuf; imgbuf = NULL; }
+ if (maskbuf != NULL) { delete maskbuf; maskbuf = NULL; }
}
else
{
- PYTHON_ERR( PyExc_StandardError );
+ PYTHON_CPL_ERR( PyExc_StandardError );
}
+#endif
getOut:
GDALClose( hDstDS );
GDALClose( hSrcDS );
- if ( papszWarpOptions ) CSLDestroy( papszWarpOptions );
- if ( pszTargetSRS != NULL ) CPLFree(pszTargetSRS);
- if ( pszSourceSRS != NULL ) CPLFree(pszSourceSRS);
+ if (psWO != NULL) GDALDestroyWarpOptions( psWO );
+
+ //GDALDumpOpenDatasets( stderr );
- GDALDumpOpenDatasets( stderr );
GDALDestroyDriverManager();
- CPLPopErrorHandler();
if ( !err && CPLGetLastErrorNo() )
{
- PYTHON_ERR( PyExc_StandardError );
- err = CPLGetLastErrorNo();
+ PYTHON_CPL_ERR( PyExc_StandardError );
+ return NULL;
}
if ( err )
{
- //CPLError( CE_Failure, err, "%s\n", CPLGetLastErrorMsg() );
return NULL;
}
else
{
- return pyImageData;
+ return pyReturnData;
}
}
@@ -573,7 +827,8 @@
GDALWarpCreateOutput( GDALDatasetH hSrcDS, const char *pszFilename,
const char *pszFormat, const char *pszSourceSRS,
const char *pszTargetSRS, int nOrder,
- char **papszCreateOptions )
+ char **papszCreateOptions, GDALDataType eDT )
+
{
GDALDriverH hDriver;
@@ -581,7 +836,9 @@
void *hTransformArg;
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
- GDALColorTableH hCT;
+
+ if( eDT == GDT_Unknown )
+ eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
/* -------------------------------------------------------------------- */
/* Find the output driver. */
@@ -590,14 +847,13 @@
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++ )
{
@@ -605,14 +861,14 @@
if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
{
- /*
printf( " %s: %s\n",
GDALGetDriverShortName( hDriver ),
GDALGetDriverLongName( hDriver ) );
- */
}
}
- //printf( "\n" );
+ printf( "\n" );
+#endif
+
return NULL;
}
@@ -644,7 +900,7 @@
/* -------------------------------------------------------------------- */
if( dfXRes != 0.0 && dfYRes != 0.0 )
{
- CPLAssert( nPixels == 0 && nLines == 0 );
+ CPLAssert( nForcePixels == 0 && nForceLines == 0 );
if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
{
dfMinX = adfDstGeoTransform[0];
@@ -696,15 +952,22 @@
}
/* -------------------------------------------------------------------- */
-/* Create the output file. */
+/* Do we want to generate an alpha band in the output file? */
/* -------------------------------------------------------------------- */
- //printf( "Creating output file is that %dP x %dL.\n", nPixels, nLines );
+ int nDstBandCount = GDALGetRasterCount(hSrcDS);
- hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
- GDALGetRasterCount(hSrcDS),
- GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1)),
- papszCreateOptions );
+ if( bEnableSrcAlpha )
+ nDstBandCount--;
+ if( bEnableDstAlpha )
+ nDstBandCount++;
+
+/* -------------------------------------------------------------------- */
+/* Create the output file. */
+/* -------------------------------------------------------------------- */
+ hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
+ nDstBandCount, eDT, papszCreateOptions );
+
if( hDstDS == NULL )
return NULL;
@@ -715,8 +978,21 @@
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 );
More information about the Thuban-devel
mailing list
This site is hosted by Intevation GmbH (Datenschutzerklärung und Impressum | Privacy Policy and Imprint)