service.c @ 940

Last change on this file since 940 was 369, checked in by djay, 12 years ago
Add GCP option to gdal_translate service. Add warp service based on gdal_warp.
File size: 81.0 KB

Rev	Line
[369]	1	/******************************************************************************
	2	* $Id: gdalwarp.cpp 24214 2012-04-08 20:17:17Z etourigny $
	3	*
	4	* Project: High Performance Image Reprojector
	5	* Purpose: Test program for high performance warper API.
	6	* Author: Frank Warmerdam <warmerdam@pobox.com>
	7	*
	8	******************************************************************************
	9	* Copyright (c) 2002, i3 - information integration and imaging
	10	* Fort Collin, CO
	11	*
	12	* Permission is hereby granted, free of charge, to any person obtaining a
	13	* copy of this software and associated documentation files (the "Software"),
	14	* to deal in the Software without restriction, including without limitation
	15	* the rights to use, copy, modify, merge, publish, distribute, sublicense,
	16	* and/or sell copies of the Software, and to permit persons to whom the
	17	* Software is furnished to do so, subject to the following conditions:
	18	*
	19	* The above copyright notice and this permission notice shall be included
	20	* in all copies or substantial portions of the Software.
	21	*
	22	* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	23	* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	24	* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	25	* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	26	* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
	27	* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	28	* DEALINGS IN THE SOFTWARE.
	29	****************************************************************************/
	30
	31	#include "gdalwarper.h"
	32	#include "cpl_string.h"
	33	#include "ogr_spatialref.h"
	34	#include "ogr_api.h"
	35	#include <vector>
	36	#include "service.h"
	37
	38	CPL_CVSID("$Id: gdalwarp.cpp 24214 2012-04-08 20:17:17Z etourigny $");
	39
	40	extern "C"{
	41	static void
	42	LoadCutline( const char pszCutlineDSName, const char pszCLayer,
	43	const char pszCWHERE, const char pszCSQL,
	44	void **phCutlineRet );
	45	static void
	46	TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
	47	char *ppapszWarpOptions, char papszTO );
	48
	49	static GDALDatasetH
	50	GDALWarpCreateOutput( maps& conf,char papszSrcFiles, const char pszFilename,
	51	const char pszFormat, char *papszTO,
	52	char ***ppapszCreateOptions, GDALDataType eDT,
	53	void ** phTransformArg,
	54	GDALDatasetH* phSrcDS );
	55
	56	static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
	57	static double dfXRes=0.0, dfYRes=0.0;
	58	static int bTargetAlignedPixels = FALSE;
	59	static int nForcePixels=0, nForceLines=0, bQuiet = FALSE;
	60	static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
	61
	62	static int bVRT = FALSE;
	63
	64	/******************************************************************************/
	65	/*! \page gdalwarp gdalwarp
	66
	67	image reprojection and warping utility
	68
	69	\section gdalwarp_synopsis SYNOPSIS
	70
	71	\htmlonly
	72	Usage:
	73	\endhtmlonly
	74
	75	\verbatim
	76	gdalwarp [--help-general] [--formats]
	77	[-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]
	78	[-order n \| -tps \| -rpc \| -geoloc] [-et err_threshold]
	79	[-refine_gcps tolerance [minimum_gcps]]
	80	[-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
	81	[-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
	82	[-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha
	83	[-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
	84	[-cutline datasource] [-cl layer] [-cwhere expression]
	85	[-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
	86	[-of format] [-co "NAME=VALUE"]* [-overwrite]
	87	srcfile* dstfile
	88	\endverbatim
	89
	90	\section gdalwarp_description DESCRIPTION
	91
	92	<p>
	93	The gdalwarp utility is an image mosaicing, reprojection and warping
	94	utility. The program can reproject to any supported projection,
	95	and can also apply GCPs stored with the image if the image is "raw"
	96	with control information.
	97
	98	<p>
	99	<dl>
	100	<dt> <b>-s_srs</b> <em>srs def</em>:</dt><dd> source spatial reference set.
	101	The coordinate systems that can be passed are anything supported by the
	102	OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes
	103	(ie. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file
	104	containing well known text.</dd>
	105	<dt> <b>-t_srs</b> <em>srs_def</em>:</dt><dd> target spatial reference set.
	106	The coordinate systems that can be passed are anything supported by the
	107	OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes
	108	(ie. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file
	109	containing well known text.</dd>
	110	<dt> <b>-to</b> <em>NAME=VALUE</em>:</dt><dd> set a transformer option suitable
	111	to pass to GDALCreateGenImgProjTransformer2(). </dd>
	112	<dt> <b>-order</b> <em>n</em>:</dt><dd> order of polynomial used for warping
	113	(1 to 3). The default is to select a polynomial order based on the number of
	114	GCPs.</dd>
	115	<dt> <b>-tps</b>:</dt><dd>Force use of thin plate spline transformer based on
	116	available GCPs.</dd>
	117	<dt> <b>-rpc</b>:</dt> <dd>Force use of RPCs.</dd>
	118	<dt> <b>-geoloc</b>:</dt><dd>Force use of Geolocation Arrays.</dd>
	119	<dt> <b>-et</b> <em>err_threshold</em>:</dt><dd> error threshold for
	120	transformation approximation (in pixel units - defaults to 0.125).</dd>
	121	<dt> <b>-refine_gcps</b> <em>tolerance minimum_gcps</em>:</dt><dd> (GDAL >= 1.9.0) refines the GCPs by automatically eliminating outliers.
	122	Outliers will be eliminated until minimum_gcps are left or when no outliers can be detected.
	123	The tolerance is passed to adjust when a GCP will be eliminated.
	124	Not that GCP refinement only works with polynomial interpolation.
	125	The tolerance is in pixel units if no projection is available, otherwise it is in SRS units.
	126	If minimum_gcps is not provided, the minimum GCPs according to the polynomial model is used.</dd>
	127	<dt> <b>-te</b> <em>xmin ymin xmax ymax</em>:</dt><dd> set georeferenced
	128	extents of output file to be created (in target SRS).</dd>
	129	<dt> <b>-tr</b> <em>xres yres</em>:</dt><dd> set output file resolution (in
	130	target georeferenced units)</dd>
	131	<dt> <b>-tap</b>:</dt><dd> (GDAL >= 1.8.0) (target aligned pixels) align
	132	the coordinates of the extent of the output file to the values of the -tr,
	133	such that the aligned extent includes the minimum extent.</dd>
	134	<dt> <b>-ts</b> <em>width height</em>:</dt><dd> set output file size in
	135	pixels and lines. If width or height is set to 0, the other dimension will be
	136	guessed from the computed resolution. Note that -ts cannot be used with -tr</dd>
	137	<dt> <b>-wo</b> <em>"NAME=VALUE"</em>:</dt><dd> Set a warp options. The
	138	GDALWarpOptions::papszWarpOptions docs show all options. Multiple
	139	<b>-wo</b> options may be listed.</dd>
	140	<dt> <b>-ot</b> <em>type</em>:</dt><dd> For the output bands to be of the
	141	indicated data type.</dd>
	142	<dt> <b>-wt</b> <em>type</em>:</dt><dd> Working pixel data type. The data type
	143	of pixels in the source image and destination image buffers.</dd>
	144	<dt> <b>-srcnodata</b> <em>value [value...]</em>:</dt><dd> Set nodata masking
	145	values for input bands (different values can be supplied for each band). If
	146	more than one value is supplied all values should be quoted to keep them
	147	together as a single operating system argument. Masked values will not be
	148	used in interpolation. Use a value of <tt>None</tt> to ignore intrinsic nodata settings on the source dataset.</dd>
	149	<dt> <b>-dstnodata</b> <em>value [value...]</em>:</dt><dd> Set nodata values
	150	for output bands (different values can be supplied for each band). If more
	151	than one value is supplied all values should be quoted to keep them together
	152	as a single operating system argument. New files will be initialized to this
	153	value and if possible the nodata value will be recorded in the output
	154	file.</dd>
	155	<dt> <b>-dstalpha</b>:</dt><dd> Create an output alpha band to identify
	156	nodata (unset/transparent) pixels. </dd>
	157	<dt> <b>-wm</b> <em>memory_in_mb</em>:</dt><dd> Set the amount of memory (in
	158	megabytes) that the warp API is allowed to use for caching.</dd>
	159	<dt> <b>-multi</b>:</dt><dd> Use multithreaded warping implementation.
	160	Multiple threads will be used to process chunks of image and perform
	161	input/output operation simultaneously.</dd>
	162	<dt> <b>-q</b>:</dt><dd> Be quiet.</dd>
	163	<dt> <b>-of</b> <em>format</em>:</dt><dd> Select the output format. The default is GeoTIFF (GTiff). Use the short format name. </dd>
	164	<dt> <b>-co</b> <em>"NAME=VALUE"</em>:</dt><dd> passes a creation option to
	165	the output format driver. Multiple <b>-co</b> options may be listed. See
	166	format specific documentation for legal creation options for each format.
	167	</dd>
	168
	169	<dt> <b>-cutline</b> <em>datasource</em>:</dt><dd>Enable use of a blend cutline from the name OGR support datasource.</dd>
	170	<dt> <b>-cl</b> <em>layername</em>:</dt><dd>Select the named layer from the
	171	cutline datasource.</dd>
	172	<dt> <b>-cwhere</b> <em>expression</em>:</dt><dd>Restrict desired cutline features based on attribute query.</dd>
	173	<dt> <b>-csql</b> <em>query</em>:</dt><dd>Select cutline features using an SQL query instead of from a layer with -cl.</dd>
	174	<dt> <b>-cblend</b> <em>distance</em>:</dt><dd>Set a blend distance to use to blend over cutlines (in pixels).</dd>
	175	<dt> <b>-crop_to_cutline</b>:</dt><dd>(GDAL >= 1.8.0) Crop the extent of the target dataset to the extent of the cutline.</dd>
	176	<dt> <b>-overwrite</b>:</dt><dd>(GDAL >= 1.8.0) Overwrite the target dataset if it already exists.</dd>
	177
	178	<dt> <em>srcfile</em>:</dt><dd> The source file name(s). </dd>
	179	<dt> <em>dstfile</em>:</dt><dd> The destination file name. </dd>
	180	</dl>
	181
	182	Mosaicing into an existing output file is supported if the output file
	183	already exists. The spatial extent of the existing file will not
	184	be modified to accomodate new data, so you may have to remove it in that case, or
	185	use the -overwrite option.
	186
	187	Polygon cutlines may be used as a mask to restrict the area of the destination file
	188	that may be updated, including blending. If the OGR layer containing the cutline
	189	features has no explicit SRS, the cutline features must be in the georeferenced
	190	units of the destination file. When outputing to a not yet existing target dataset,
	191	its extent will be the one of the original raster unless -te or -crop_to_cutline are
	192	specified.
	193
	194	<p>
	195	\section gdalwarp_example EXAMPLE
	196
	197	For instance, an eight bit spot scene stored in GeoTIFF with
	198	control points mapping the corners to lat/long could be warped to a UTM
	199	projection with a command like this:<p>
	200
	201	\verbatim
	202	gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif
	203	\endverbatim
	204
	205	For instance, the second channel of an ASTER image stored in HDF with
	206	control points mapping the corners to lat/long could be warped to a UTM
	207	projection with a command like this:<p>
	208
	209	\verbatim
	210	gdalwarp HDF4_SDS:ASTER_L1B:"pg-PR1B0000-2002031402_100_001":2 pg-PR1B0000-2002031402_100_001_2.tif
	211	\endverbatim
	212
	213	\if man
	214	\section gdalwarp_author AUTHORS
	215	Frank Warmerdam <warmerdam@pobox.com>, Silke Reimer <silke@intevation.de>
	216	\endif
	217	*/
	218
	219	/************************************************************************/
	220	/* GDALExit() */
	221	/* This function exits and cleans up GDAL and OGR resources */
	222	/* Perhaps it should be added to C api and used in all apps? */
	223	/************************************************************************/
	224
	225	static int GDALExit( int nCode )
	226	{
	227	const char *pszDebug = CPLGetConfigOption("CPL_DEBUG",NULL);
	228	if( pszDebug && (EQUAL(pszDebug,"ON") \|\| EQUAL(pszDebug,"") ) )
	229	{
	230	GDALDumpOpenDatasets( stderr );
	231	CPLDumpSharedList( NULL );
	232	}
	233
	234	GDALDestroyDriverManager();
	235
	236	#ifdef OGR_ENABLED
	237	OGRCleanupAll();
	238	#endif
	239
	240	return SERVICE_FAILED;
	241	}
	242
	243	/************************************************************************/
	244	/* Usage() */
	245	/************************************************************************/
	246
	247	static void Usage()
	248
	249	{
	250	printf(
	251	"Usage: gdalwarp [--help-general] [--formats]\n"
	252	" [-s_srs srs_def] [-t_srs srs_def] [-to \"NAME=VALUE\"]\n"
	253	" [-order n \| -tps \| -rpc \| -geoloc] [-et err_threshold]\n"
	254	" [-refine_gcps tolerance [minimum_gcps]]\n"
	255	" [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]\n"
	256	" [-wo \"NAME=VALUE\"] [-ot Byte/Int16/...] [-wt Byte/Int16]\n"
	257	" [-srcnodata \"value [value...]\"] [-dstnodata \"value [value...]\"] -dstalpha\n"
	258	" [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n"
	259	" [-cutline datasource] [-cl layer] [-cwhere expression]\n"
	260	" [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]\n"
	261	" [-of format] [-co \"NAME=VALUE\"]* [-overwrite]\n"
	262	" srcfile* dstfile\n"
	263	"\n"
	264	"Available resampling methods:\n"
	265	" near (default), bilinear, cubic, cubicspline, lanczos.\n" );
	266	GDALExit( 1 );
	267	}
	268
	269	/************************************************************************/
	270	/* SanitizeSRS */
	271	/************************************************************************/
	272
	273	char SanitizeSRS( const char pszUserInput )
	274
	275	{
	276	OGRSpatialReferenceH hSRS;
	277	char *pszResult = NULL;
	278
	279	CPLErrorReset();
	280
	281	hSRS = OSRNewSpatialReference( NULL );
	282	if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
	283	OSRExportToWkt( hSRS, &pszResult );
	284	else
	285	{
	286	CPLError( CE_Failure, CPLE_AppDefined,
	287	"Translating source or target SRS failed:\n%s",
	288	pszUserInput );
	289	GDALExit( 1 );
	290	}
	291
	292	OSRDestroySpatialReference( hSRS );
	293
	294	return pszResult;
	295	}
	296
	297	/************************************************************************/
	298	/* main() */
	299	/************************************************************************/
	300
	301	#ifdef WIN32
	302	__declspec(dllexport)
	303	#endif
	304	int Gdal_Warp( maps& conf,maps& inputs,maps*& outputs )
	305
	306	{
	307	GDALDatasetH hDstDS;
	308	const char *pszFormat = "GTiff";
	309	int bFormatExplicitelySet = FALSE;
	310	char **papszSrcFiles = NULL;
	311	char *pszDstFilename = NULL;
	312	int bCreateOutput = FALSE, i;
	313	void hTransformArg, hGenImgProjArg=NULL, *hApproxArg=NULL;
	314	char **papszWarpOptions = NULL;
	315	double dfErrorThreshold = 0.125;
	316	double dfWarpMemoryLimit = 0.0;
	317	GDALTransformerFunc pfnTransformer = NULL;
	318	char **papszCreateOptions = NULL;
	319	GDALDataType eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown;
	320	GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
	321	const char *pszSrcNodata = NULL;
	322	const char *pszDstNodata = NULL;
	323	int bMulti = FALSE;
	324	char **papszTO = NULL;
	325	char *pszCutlineDSName = NULL;
	326	char pszCLayer = NULL, pszCWHERE = NULL, *pszCSQL = NULL;
	327	void *hCutline = NULL;
	328	int bHasGotErr = FALSE;
	329	int bCropToCutline = FALSE;
	330	int bOverwrite = FALSE;
	331
	332
	333	/* Check that we are running against at least GDAL 1.6 */
	334	/* Note to developers : if we use newer API, please change the requirement */
	335	if (atoi(GDALVersionInfo("VERSION_NUM")) < 1600)
	336	{
	337	fprintf(stderr, "At least, GDAL >= 1.6.0 is required for this version of warp_service.zo, "
	338	"which was compiled against GDAL %s\n", GDAL_RELEASE_NAME);
	339	GDALExit(1);
	340	}
	341
	342	/* -------------------------------------------------------------------- */
	343	/* Register standard GDAL drivers, and process generic GDAL */
	344	/* command options. */
	345	/* -------------------------------------------------------------------- */
	346	GDALAllRegister();
	347
	348	/* -------------------------------------------------------------------- */
	349	/* Parse arguments. */
	350	/* -------------------------------------------------------------------- */
	351	char *dataPath;
	352	map* tmpMap=getMapFromMaps(conf,"main","dataPath");
	353	if(tmpMap!=NULL)
	354	dataPath=strdup(tmpMap->value);
	355	tmpMap=NULL;
	356
	357	char *tempPath;
	358	tmpMap=getMapFromMaps(conf,"main","tmpPath");
	359	if(tmpMap!=NULL){
	360	tempPath=strdup(tmpMap->value);
	361	}
	362
	363	tmpMap=NULL;
	364	tmpMap=getMapFromMaps(inputs,"InputDSN","value");
	365	if(tmpMap!=NULL){
	366	map* length=getMapFromMaps(inputs,"InputDSN","length");
	367	if(length!=NULL){
	368	int i=0;
	369	int len=atoi(length->value);
	370	maps* currentMaps=getMaps(inputs,"InputDSN");
	371	for(i=0;i<len;i++){
	372	tmpMap=getMapArray(currentMaps->content,"value",i);
	373	papszSrcFiles[i]=(char)CPLMalloc(sizeof(char)(strlen(dataPath)+strlen(tmpMap->value)+4));
	374	sprintf((char*)papszSrcFiles[i],"%s/%s",dataPath,tmpMap->value);
	375	}
	376	}else{
	377	char tmp=(char)CPLMalloc(sizeof(char)*(strlen(dataPath)+strlen(tmpMap->value)+4));
	378	sprintf(tmp,"%s/%s",dataPath,tmpMap->value);
	379	papszSrcFiles=CSLAddString( papszSrcFiles, tmp );
	380	free(tmp);
	381	}
	382	}
	383
	384	tmpMap=NULL;
	385	tmpMap=getMapFromMaps(inputs,"OutputDSN","value");
	386	if(tmpMap!=NULL){
	387	pszDstFilename=(char)CPLMalloc(sizeof(char)(strlen(tempPath)+strlen(tmpMap->value)+4));
	388	char *ext=new char[4];
	389	ext="tif";
	390	if(strncasecmp(pszFormat,"AAIGRID",7)==0)
	391	ext="csv";
	392	else
	393	if(strncasecmp(pszFormat,"PNG",3)==0)
	394	ext="png";
	395	else
	396	if(strncasecmp(pszFormat,"GIF",3)==0)
	397	ext="gif";
	398	else
	399	if(strncasecmp(pszFormat,"JPEG",4)==0)
	400	ext="jpg";
	401	sprintf(pszDstFilename,"%s/%s.%s",tempPath,tmpMap->value,ext);
	402	}
	403
	404
	405	tmpMap=NULL;
	406	tmpMap=getMapFromMaps(inputs,"r","value");
	407	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	408	if ( EQUAL(tmpMap->value, "near") )
	409	eResampleAlg = GRA_NearestNeighbour;
	410	else if ( EQUAL(tmpMap->value, "bilinear") )
	411	eResampleAlg = GRA_Bilinear;
	412	else if ( EQUAL(tmpMap->value, "cubic") )
	413	eResampleAlg = GRA_Cubic;
	414	else if ( EQUAL(tmpMap->value, "cubicspline") )
	415	eResampleAlg = GRA_CubicSpline;
	416	else if ( EQUAL(tmpMap->value, "lanczos") )
	417	eResampleAlg = GRA_Lanczos;
	418	else
	419	{
	420	char tmp[1024];
	421	sprintf(tmp,"Unknown resampling method: \"%s\".", tmpMap->value );
	422	setMapInMaps(conf,"lenv","message",tmp);
	423	return SERVICE_FAILED;
	424	}
	425	}
	426
	427	tmpMap=NULL;
	428	tmpMap=getMapFromMaps(inputs,"t_srs","value");
	429	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	430	char *pszSRS = SanitizeSRS(tmpMap->value);
	431	papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszSRS );
	432	CPLFree( pszSRS );
	433	}
	434
	435	tmpMap=getMapFromMaps(inputs,"s_srs","value");
	436	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	437	char *pszSRS = SanitizeSRS(tmpMap->value);
	438	papszTO = CSLSetNameValue( papszTO, "SRC_SRS", pszSRS );
	439	CPLFree( pszSRS );
	440	}
	441
	442	tmpMap=getMapFromMaps(inputs,"order","value");
	443	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	444	const char* pszMethod = CSLFetchNameValue(papszTO, "METHOD");
	445	if (pszMethod)
	446	fprintf(stderr, "Warning: only one METHOD can be used. Method %s is already defined\n",
	447	pszMethod);
	448	papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", tmpMap->value );
	449	}
	450
	451	tmpMap=getMapFromMaps(inputs,"co","value");
	452	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	453	papszCreateOptions = CSLAddString( papszCreateOptions, tmpMap->value );
	454	bCreateOutput = TRUE;
	455	}
	456
	457	tmpMap=getMapFromMaps(inputs,"wo","value");
	458	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	459	papszWarpOptions = CSLAddString( papszWarpOptions, tmpMap->value );
	460	}
	461
	462	tmpMap=getMapFromMaps(inputs,"multi","value");
	463	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	464	bMulti = TRUE;
	465	}
	466
	467	tmpMap=getMapFromMaps(inputs,"dstalpha","value");
	468	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	469	bEnableDstAlpha = TRUE;
	470	}
	471
	472	tmpMap=getMapFromMaps(inputs,"srcalpha","value");
	473	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	474	bEnableSrcAlpha = TRUE;
	475	}
	476
	477	tmpMap=getMapFromMaps(inputs,"tap","value");
	478	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	479	bTargetAlignedPixels = TRUE;
	480	}
	481
	482
	483	tmpMap=getMapFromMaps(inputs,"crop_to_cutline","value");
	484	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	485	bCropToCutline = TRUE;
	486	bCreateOutput = TRUE;
	487	}
	488
	489
	490	tmpMap=getMapFromMaps(inputs,"tps","value");
	491	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	492	papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
	493	}
	494
	495	tmpMap=getMapFromMaps(inputs,"rpc","value");
	496	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	497	papszTO = CSLSetNameValue( papszTO, "METHOD", "RPC" );
	498	}
	499
	500	tmpMap=getMapFromMaps(inputs,"geoloc","value");
	501	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"TRUE",4)==0 ){
	502	papszTO = CSLSetNameValue( papszTO, "METHOD", "GEOLOC_ARRAY" );
	503	}
	504
	505	tmpMap=getMapFromMaps(inputs,"wt","value");
	506	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	507	int iType;
	508
	509	for( iType = 1; iType < GDT_TypeCount; iType++ )
	510	{
	511	if( GDALGetDataTypeName((GDALDataType)iType) != NULL
	512	&& EQUAL(GDALGetDataTypeName((GDALDataType)iType),
	513	tmpMap->value) )
	514	{
	515	eWorkingType = (GDALDataType) iType;
	516	}
	517	}
	518
	519	if( eWorkingType == GDT_Unknown )
	520	{
	521	char tmp[1024];
	522	sprintf(tmp,"Unknown output pixel type: %s\n", tmpMap->value );
	523	setMapInMaps(conf,"lenv","message",tmp);
	524	return SERVICE_FAILED;
	525	}
	526	}
	527
	528	tmpMap=getMapFromMaps(inputs,"srcnodata","value");
	529	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	530	pszSrcNodata = strdup(tmpMap->value);
	531	}
	532
	533	tmpMap=getMapFromMaps(inputs,"dstnodata","value");
	534	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	535	pszDstNodata = strdup(tmpMap->value);
	536	}
	537
	538
	539	tmpMap=getMapFromMaps(inputs,"tr","value");
	540	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	541	char *tmpV=strdup(tmpMap->value);
	542	char *res=strtok(tmpV,",");
	543	int i=0;
	544	while(res!=NULL){
	545	if(i==0)
	546	dfXRes = CPLAtofM(res);
	547	else
	548	dfYRes = fabs(CPLAtofM(res));
	549	res=strtok(NULL,",");
	550	i++;
	551	}
	552	if( dfXRes == 0 \|\| dfYRes == 0 )
	553	{
	554	setMapInMaps(conf,"lenv","message","Wrong value for TR parameters");
	555	}
	556	bCreateOutput = TRUE;
	557	}
	558
	559	tmpMap=getMapFromMaps(inputs,"cblend","value");
	560	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	561	papszWarpOptions =
	562	CSLSetNameValue( papszWarpOptions,
	563	"CUTLINE_BLEND_DIST", tmpMap->value );
	564	}
	565
	566	tmpMap=getMapFromMaps(inputs,"ot","value");
	567	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	568	int iType;
	569
	570	for( iType = 1; iType < GDT_TypeCount; iType++ )
	571	{
	572	if( GDALGetDataTypeName((GDALDataType)iType) != NULL
	573	&& EQUAL(GDALGetDataTypeName((GDALDataType)iType),
	574	tmpMap->value) )
	575	{
	576	eOutputType = (GDALDataType) iType;
	577	}
	578	}
	579
	580	if( eOutputType == GDT_Unknown )
	581	{
	582	char tmp[128];
	583	sprintf( tmp,"Unknown output pixel type: %s\n", tmpMap->value );
	584	setMapInMaps(conf,"lenv","message",tmp);
	585	return SERVICE_FAILED;
	586	}
	587	bCreateOutput = TRUE;
	588	}
	589
	590
	591	tmpMap=getMapFromMaps(inputs,"ts","value");
	592	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	593	char *tmpV=strdup(tmpMap->value);
	594	char *res=strtok(tmpV,",");
	595	int i=0;
	596	while(res!=NULL){
	597	if(i==0)
	598	nForcePixels = atoi(res);
	599	else
	600	nForceLines = atoi(res);
	601	res=strtok(NULL,",");
	602	i++;
	603	}
	604	free(tmpV);
	605	bCreateOutput = TRUE;
	606	}
	607
	608	tmpMap=getMapFromMaps(inputs,"te","value");
	609	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	610	char *tmpV=strdup(tmpMap->value);
	611	char *res=strtok(tmpV,",");
	612	int i=0;
	613	while(res!=NULL){
	614	switch(i){
	615	case 0:
	616	dfMinX = CPLAtofM(res);
	617	break;
	618	case 1:
	619	dfMinY = CPLAtofM(res);
	620	break;
	621	case 2:
	622	dfMaxX = CPLAtofM(res);
	623	break;
	624	case 3:
	625	dfMaxY = CPLAtofM(res);
	626	break;
	627	}
	628	if(i==0)
	629	nForcePixels = atoi(res);
	630	else
	631	nForceLines = atoi(res);
	632	res=strtok(NULL,",");
	633	i++;
	634	}
	635	free(tmpV);
	636	bCreateOutput = TRUE;
	637	}
	638
	639	tmpMap=getMapFromMaps(inputs,"te","value");
	640	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	641	pszCutlineDSName = strdup(tmpMap->value);
	642	}
	643
	644	tmpMap=getMapFromMaps(inputs,"cwhere","value");
	645	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	646	pszCWHERE = strdup(tmpMap->value);
	647	}
	648
	649	tmpMap=getMapFromMaps(inputs,"cl","value");
	650	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	651	pszCLayer = strdup(tmpMap->value);
	652	}
	653
	654	tmpMap=getMapFromMaps(inputs,"csql","value");
	655	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	656	pszCSQL = strdup(tmpMap->value);
	657	}
	658
	659	tmpMap=getMapFromMaps(inputs,"of","value");
	660	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	661	pszFormat = strdup(tmpMap->value);
	662	bFormatExplicitelySet = TRUE;
	663	bCreateOutput = TRUE;
	664	if( EQUAL(pszFormat,"VRT") )
	665	bVRT = TRUE;
	666	}
	667
	668	tmpMap=getMapFromMaps(inputs,"to","value");
	669	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	670	papszTO = CSLAddString( papszTO, tmpMap->value );
	671	}
	672
	673	tmpMap=getMapFromMaps(inputs,"et","value");
	674	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	675	dfErrorThreshold = CPLAtofM(tmpMap->value);
	676	}
	677
	678
	679	tmpMap=getMapFromMaps(inputs,"refine_gpcs","value");
	680	if(tmpMap!=NULL && strncasecmp(tmpMap->value,"NULL",4)!=0 ){
	681	char *tmpV=strdup(tmpMap->value);
	682	char *res=strtok(tmpV,",");
	683	int i=0;
	684	while(res!=NULL){
	685	if(i==0){
	686	papszTO = CSLSetNameValue( papszTO, "REFINE_TOLERANCE", res );
	687	if(atof(res) < 0){
	688	setMapInMaps(conf,"lenv","message","The tolerance for -refine_gcps may not be negative\n");
	689	return SERVICE_FAILED;
	690	}
	691	}else if(i==1){
	692	if(atoi(res) >= 0 && isdigit(res[0]))
	693	papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", res );
	694	else
	695	papszTO = CSLSetNameValue( papszTO, "REFINE_MINIMUM_GCPS", "-1" );
	696	}
	697	res=strtok(NULL,",");
	698	i++;
	699	}
	700	}
	701
	702	bQuiet = TRUE;
	703
	704	/* -------------------------------------------------------------------- */
	705	/* Check that incompatible options are not used */
	706	/* -------------------------------------------------------------------- */
	707
	708	if ((nForcePixels != 0 \|\| nForceLines != 0) &&
	709	(dfXRes != 0 && dfYRes != 0))
	710	{
	711	setMapInMaps(conf,"lenv","message","tr and ts options cannot be used at the same time\n");
	712	return SERVICE_FAILED;
	713	}
	714
	715	if (bTargetAlignedPixels && dfXRes == 0 && dfYRes == 0)
	716	{
	717	setMapInMaps(conf,"lenv","message","tap option cannot be used without using tr");
	718	return SERVICE_FAILED;
	719	}
	720
	721
	722	/* -------------------------------------------------------------------- */
	723	/* Does the output dataset already exist? */
	724	/* -------------------------------------------------------------------- */
	725
	726	/* FIXME ? source filename=target filename and -overwrite is definitely */
	727	/* an error. But I can't imagine of a valid case (without -overwrite), */
	728	/* where it would make sense. In doubt, let's keep that dubious possibility... */
	729	if ( CSLCount(papszSrcFiles) == 1 &&
	730	strcmp(papszSrcFiles[0], pszDstFilename) == 0 && bOverwrite)
	731	{
	732	setMapInMaps(conf,"lenv","message","Source and destination datasets must be different.");
	733	return GDALExit(1);
	734	}
	735
	736	CPLPushErrorHandler( CPLQuietErrorHandler );
	737	hDstDS = GDALOpen( pszDstFilename, GA_Update );
	738	CPLPopErrorHandler();
	739
	740	if( hDstDS != NULL && bOverwrite )
	741	{
	742	GDALClose(hDstDS);
	743	hDstDS = NULL;
	744	}
	745
	746	if( hDstDS != NULL && bCreateOutput )
	747	{
	748	char tmp[024];
	749	sprintf( tmp,
	750	"Output dataset %s exists,\n"
	751	"but some commandline options were provided indicating a new dataset\n"
	752	"should be created. Please delete existing dataset and run again.\n",
	753	pszDstFilename );
	754	setMapInMaps(conf,"lenv","message",tmp);
	755	return GDALExit(1);
	756	}
	757
	758	/* Avoid overwriting an existing destination file that cannot be opened in */
	759	/* update mode with a new GTiff file */
	760	if ( hDstDS == NULL && !bOverwrite )
	761	{
	762	CPLPushErrorHandler( CPLQuietErrorHandler );
	763	hDstDS = GDALOpen( pszDstFilename, GA_ReadOnly );
	764	CPLPopErrorHandler();
	765
	766	if (hDstDS)
	767	{
	768	char tmp[1024];
	769	sprintf( tmp,
	770	"Output dataset %s exists, but cannot be opened in update mode\n",
	771	pszDstFilename );
	772	GDALClose(hDstDS);
	773	setMapInMaps(conf,"lenv","message",tmp);
	774	return GDALExit( 1 );
	775	}
	776	}
	777
	778	/* -------------------------------------------------------------------- */
	779	/* If we have a cutline datasource read it and attach it in the */
	780	/* warp options. */
	781	/* -------------------------------------------------------------------- */
	782	if( pszCutlineDSName != NULL )
	783	{
	784	LoadCutline( pszCutlineDSName, pszCLayer, pszCWHERE, pszCSQL,
	785	&hCutline );
	786	}
	787
	788	#ifdef OGR_ENABLED
	789	if ( bCropToCutline && hCutline != NULL )
	790	{
	791
	792
	793	OGRGeometryH hCutlineGeom = OGR_G_Clone( (OGRGeometryH) hCutline );
	794	OGRSpatialReferenceH hCutlineSRS = OGR_G_GetSpatialReference( hCutlineGeom );
	795	const char *pszThisTargetSRS = CSLFetchNameValue( papszTO, "DST_SRS" );
	796	OGRCoordinateTransformationH hCT = NULL;
	797	if (hCutlineSRS == NULL)
	798	{
	799	/* We suppose it is in target coordinates */
	800	}
	801	else if (pszThisTargetSRS != NULL)
	802	{
	803	OGRSpatialReferenceH hTargetSRS = OSRNewSpatialReference(NULL);
	804	if( OSRImportFromWkt( hTargetSRS, (char **)&pszThisTargetSRS ) != CE_None )
	805	{
	806	setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.\n");
	807	return GDALExit(1);
	808	}
	809
	810	hCT = OCTNewCoordinateTransformation(hCutlineSRS, hTargetSRS);
	811
	812	OSRDestroySpatialReference(hTargetSRS);
	813	}
	814	else if (pszThisTargetSRS == NULL)
	815	{
	816	if (papszSrcFiles[0] != NULL)
	817	{
	818	GDALDatasetH hSrcDS = GDALOpen(papszSrcFiles[0], GA_ReadOnly);
	819	if (hSrcDS == NULL)
	820	{
	821	setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.");
	822	return GDALExit(1);
	823	}
	824
	825	OGRSpatialReferenceH hRasterSRS = NULL;
	826	const char *pszProjection = NULL;
	827
	828	if( GDALGetProjectionRef( hSrcDS ) != NULL
	829	&& strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
	830	pszProjection = GDALGetProjectionRef( hSrcDS );
	831	else if( GDALGetGCPProjection( hSrcDS ) != NULL )
	832	pszProjection = GDALGetGCPProjection( hSrcDS );
	833
	834	if( pszProjection == NULL )
	835	{
	836	setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.");
	837	return GDALExit(1);
	838	}
	839
	840	hRasterSRS = OSRNewSpatialReference(NULL);
	841	if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None )
	842	{
	843	setMapInMaps(conf,"lenv","message","Cannot compute bounding box of cutline.");
	844	return GDALExit(1);
	845	}
	846
	847	hCT = OCTNewCoordinateTransformation(hCutlineSRS, hRasterSRS);
	848
	849	OSRDestroySpatialReference(hRasterSRS);
	850
	851	GDALClose(hSrcDS);
	852	}
	853	else
	854	{
	855	setMapInMaps(conf,"levn","message", "Cannot compute bounding box of cutline.\n");
	856	return GDALExit(1);
	857	}
	858	}
	859
	860	if (hCT)
	861	{
	862	OGR_G_Transform( hCutlineGeom, hCT );
	863
	864	OCTDestroyCoordinateTransformation(hCT);
	865	}
	866
	867	OGREnvelope sEnvelope;
	868	OGR_G_GetEnvelope(hCutlineGeom, &sEnvelope);
	869
	870	dfMinX = sEnvelope.MinX;
	871	dfMinY = sEnvelope.MinY;
	872	dfMaxX = sEnvelope.MaxX;
	873	dfMaxY = sEnvelope.MaxY;
	874
	875	OGR_G_DestroyGeometry(hCutlineGeom);
	876	}
	877	#endif
	878
	879	/* -------------------------------------------------------------------- */
	880	/* If not, we need to create it. */
	881	/* -------------------------------------------------------------------- */
	882	int bInitDestSetForFirst = FALSE;
	883
	884	void* hUniqueTransformArg = NULL;
	885	GDALDatasetH hUniqueSrcDS = NULL;
	886
	887	if( hDstDS == NULL )
	888	{
	889
	890	hDstDS = GDALWarpCreateOutput( conf,papszSrcFiles, pszDstFilename,pszFormat,
	891	papszTO, &papszCreateOptions,
	892	eOutputType, &hUniqueTransformArg,
	893	&hUniqueSrcDS);
	894	bCreateOutput = TRUE;
	895
	896
	897	if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL
	898	&& pszDstNodata == NULL )
	899	{
	900	papszWarpOptions = CSLSetNameValue(papszWarpOptions,
	901	"INIT_DEST", "0");
	902	bInitDestSetForFirst = TRUE;
	903	}
	904	else if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL )
	905	{
	906	papszWarpOptions = CSLSetNameValue(papszWarpOptions,
	907	"INIT_DEST", "NO_DATA" );
	908	bInitDestSetForFirst = TRUE;
	909	}
	910
	911	CSLDestroy( papszCreateOptions );
	912	papszCreateOptions = NULL;
	913	}
	914
	915
	916	if( hDstDS == NULL )
	917	GDALExit( 1 );
	918
	919	/* -------------------------------------------------------------------- */
	920	/* Loop over all source files, processing each in turn. */
	921	/* -------------------------------------------------------------------- */
	922	int iSrc;
	923
	924	for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
	925	{
	926	GDALDatasetH hSrcDS;
	927
	928	/* -------------------------------------------------------------------- */
	929	/* Open this file. */
	930	/* -------------------------------------------------------------------- */
	931	if (hUniqueSrcDS)
	932	hSrcDS = hUniqueSrcDS;
	933	else
	934	hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
	935
	936	if( hSrcDS == NULL ){
	937	setMapInMaps(conf,"lenv","message",CPLGetLastErrorMsg());
	938	return GDALExit( 2 );
	939	}
	940
	941	/* -------------------------------------------------------------------- */
	942	/* Check that there's at least one raster band */
	943	/* -------------------------------------------------------------------- */
	944	if ( GDALGetRasterCount(hSrcDS) == 0 )
	945	{
	946	char tmp[1024];
	947	sprintf(tmp, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
	948	setMapInMaps(conf,"lenv","message",tmp);
	949	return GDALExit( 1 );
	950	}
	951
	952	if( !bQuiet )
	953	printf( "Processing input file %s.\n", papszSrcFiles[iSrc] );
	954
	955	/* -------------------------------------------------------------------- */
	956	/* Warns if the file has a color table and something more */
	957	/* complicated than nearest neighbour resampling is asked */
	958	/* -------------------------------------------------------------------- */
	959
	960	if ( eResampleAlg != GRA_NearestNeighbour &&
	961	GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
	962	{
	963	fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to "
	964	"bad results when using a resampling method other than "
	965	"nearest neighbour. Converting the dataset prior to 24/32 bit "
	966	"is advised.\n", papszSrcFiles[iSrc] );
	967	}
	968
	969	/* -------------------------------------------------------------------- */
	970	/* Do we have a source alpha band? */
	971	/* -------------------------------------------------------------------- */
	972	if( GDALGetRasterColorInterpretation(
	973	GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) )
	974	== GCI_AlphaBand
	975	&& !bEnableSrcAlpha )
	976	{
	977	bEnableSrcAlpha = TRUE;
	978	if( !bQuiet )
	979	printf( "Using band %d of source image as alpha.\n",
	980	GDALGetRasterCount(hSrcDS) );
	981	}
	982
	983	/* -------------------------------------------------------------------- */
	984	/* Create a transformation object from the source to */
	985	/* destination coordinate system. */
	986	/* -------------------------------------------------------------------- */
	987	if (hUniqueTransformArg)
	988	hTransformArg = hGenImgProjArg = hUniqueTransformArg;
	989	else
	990	hTransformArg = hGenImgProjArg =
	991	GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszTO );
	992
	993	if( hTransformArg == NULL ){
	994	setMapInMaps(conf,"lenv","message",CPLGetLastErrorMsg());
	995	return GDALExit( 1 );
	996	}
	997
	998	pfnTransformer = GDALGenImgProjTransform;
	999
	1000	/* -------------------------------------------------------------------- */
	1001	/* Warp the transformer with a linear approximator unless the */
	1002	/* acceptable error is zero. */
	1003	/* -------------------------------------------------------------------- */
	1004	if( dfErrorThreshold != 0.0 )
	1005	{
	1006	hTransformArg = hApproxArg =
	1007	GDALCreateApproxTransformer( GDALGenImgProjTransform,
	1008	hGenImgProjArg, dfErrorThreshold);
	1009	pfnTransformer = GDALApproxTransform;
	1010	}
	1011
	1012	/* -------------------------------------------------------------------- */
	1013	/* Clear temporary INIT_DEST settings after the first image. */
	1014	/* -------------------------------------------------------------------- */
	1015	if( bInitDestSetForFirst && iSrc == 1 )
	1016	papszWarpOptions = CSLSetNameValue( papszWarpOptions,
	1017	"INIT_DEST", NULL );
	1018
	1019	/* -------------------------------------------------------------------- */
	1020	/* Setup warp options. */
	1021	/* -------------------------------------------------------------------- */
	1022	GDALWarpOptions *psWO = GDALCreateWarpOptions();
	1023
	1024	psWO->papszWarpOptions = CSLDuplicate(papszWarpOptions);
	1025	psWO->eWorkingDataType = eWorkingType;
	1026	psWO->eResampleAlg = eResampleAlg;
	1027
	1028	psWO->hSrcDS = hSrcDS;
	1029	psWO->hDstDS = hDstDS;
	1030
	1031	psWO->pfnTransformer = pfnTransformer;
	1032	psWO->pTransformerArg = hTransformArg;
	1033
	1034	if( !bQuiet )
	1035	psWO->pfnProgress = GDALTermProgress;
	1036
	1037	if( dfWarpMemoryLimit != 0.0 )
	1038	psWO->dfWarpMemoryLimit = dfWarpMemoryLimit;
	1039
	1040	/* -------------------------------------------------------------------- */
	1041	/* Setup band mapping. */
	1042	/* -------------------------------------------------------------------- */
	1043	if( bEnableSrcAlpha )
	1044	psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
	1045	else
	1046	psWO->nBandCount = GDALGetRasterCount(hSrcDS);
	1047
	1048	psWO->panSrcBands = (int ) CPLMalloc(psWO->nBandCountsizeof(int));
	1049	psWO->panDstBands = (int ) CPLMalloc(psWO->nBandCountsizeof(int));
	1050
	1051	for( i = 0; i < psWO->nBandCount; i++ )
	1052	{
	1053	psWO->panSrcBands[i] = i+1;
	1054	psWO->panDstBands[i] = i+1;
	1055	}
	1056
	1057	/* -------------------------------------------------------------------- */
	1058	/* Setup alpha bands used if any. */
	1059	/* -------------------------------------------------------------------- */
	1060	if( bEnableSrcAlpha )
	1061	psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
	1062
	1063	if( !bEnableDstAlpha
	1064	&& GDALGetRasterCount(hDstDS) == psWO->nBandCount+1
	1065	&& GDALGetRasterColorInterpretation(
	1066	GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS)))
	1067	== GCI_AlphaBand )
	1068	{
	1069	if( !bQuiet )
	1070	printf( "Using band %d of destination image as alpha.\n",
	1071	GDALGetRasterCount(hDstDS) );
	1072
	1073	bEnableDstAlpha = TRUE;
	1074	}
	1075
	1076	if( bEnableDstAlpha )
	1077	psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
	1078
	1079	/* -------------------------------------------------------------------- */
	1080	/* Setup NODATA options. */
	1081	/* -------------------------------------------------------------------- */
	1082	if( pszSrcNodata != NULL && !EQUALN(pszSrcNodata,"n",1) )
	1083	{
	1084	char **papszTokens = CSLTokenizeString( pszSrcNodata );
	1085	int nTokenCount = CSLCount(papszTokens);
	1086
	1087	psWO->padfSrcNoDataReal = (double *)
	1088	CPLMalloc(psWO->nBandCount*sizeof(double));
	1089	psWO->padfSrcNoDataImag = (double *)
	1090	CPLMalloc(psWO->nBandCount*sizeof(double));
	1091
	1092	for( i = 0; i < psWO->nBandCount; i++ )
	1093	{
	1094	if( i < nTokenCount )
	1095	{
	1096	CPLStringToComplex( papszTokens[i],
	1097	psWO->padfSrcNoDataReal + i,
	1098	psWO->padfSrcNoDataImag + i );
	1099	}
	1100	else
	1101	{
	1102	psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i-1];
	1103	psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i-1];
	1104	}
	1105	}
	1106
	1107	CSLDestroy( papszTokens );
	1108
	1109	psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
	1110	"UNIFIED_SRC_NODATA", "YES" );
	1111	}
	1112
	1113	/* -------------------------------------------------------------------- */
	1114	/* If -srcnodata was not specified, but the data has nodata */
	1115	/* values, use them. */
	1116	/* -------------------------------------------------------------------- */
	1117	if( pszSrcNodata == NULL )
	1118	{
	1119	int bHaveNodata = FALSE;
	1120	double dfReal = 0.0;
	1121
	1122	for( i = 0; !bHaveNodata && i < psWO->nBandCount; i++ )
	1123	{
	1124	GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
	1125	dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
	1126	}
	1127
	1128	if( bHaveNodata )
	1129	{
	1130	if( !bQuiet )
	1131	{
	1132	if (CPLIsNan(dfReal))
	1133	printf( "Using internal nodata values (eg. nan) for image %s.\n",
	1134	papszSrcFiles[iSrc] );
	1135	else
	1136	printf( "Using internal nodata values (eg. %g) for image %s.\n",
	1137	dfReal, papszSrcFiles[iSrc] );
	1138	}
	1139	psWO->padfSrcNoDataReal = (double *)
	1140	CPLMalloc(psWO->nBandCount*sizeof(double));
	1141	psWO->padfSrcNoDataImag = (double *)
	1142	CPLMalloc(psWO->nBandCount*sizeof(double));
	1143
	1144	for( i = 0; i < psWO->nBandCount; i++ )
	1145	{
	1146	GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
	1147
	1148	dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
	1149
	1150	if( bHaveNodata )
	1151	{
	1152	psWO->padfSrcNoDataReal[i] = dfReal;
	1153	psWO->padfSrcNoDataImag[i] = 0.0;
	1154	}
	1155	else
	1156	{
	1157	psWO->padfSrcNoDataReal[i] = -123456.789;
	1158	psWO->padfSrcNoDataImag[i] = 0.0;
	1159	}
	1160	}
	1161	}
	1162	}
	1163
	1164	/* -------------------------------------------------------------------- */
	1165	/* If the output dataset was created, and we have a destination */
	1166	/* nodata value, go through marking the bands with the information.*/
	1167	/* -------------------------------------------------------------------- */
	1168	if( pszDstNodata != NULL )
	1169	{
	1170	char **papszTokens = CSLTokenizeString( pszDstNodata );
	1171	int nTokenCount = CSLCount(papszTokens);
	1172
	1173	psWO->padfDstNoDataReal = (double *)
	1174	CPLMalloc(psWO->nBandCount*sizeof(double));
	1175	psWO->padfDstNoDataImag = (double *)
	1176	CPLMalloc(psWO->nBandCount*sizeof(double));
	1177
	1178	for( i = 0; i < psWO->nBandCount; i++ )
	1179	{
	1180	if( i < nTokenCount )
	1181	{
	1182	CPLStringToComplex( papszTokens[i],
	1183	psWO->padfDstNoDataReal + i,
	1184	psWO->padfDstNoDataImag + i );
	1185	}
	1186	else
	1187	{
	1188	psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i-1];
	1189	psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i-1];
	1190	}
	1191
	1192	GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 );
	1193	int bClamped = FALSE, bRounded = FALSE;
	1194
	1195	#define CLAMP(val,type,minval,maxval) \
	1196	do { if (val < minval) { bClamped = TRUE; val = minval; } \
	1197	else if (val > maxval) { bClamped = TRUE; val = maxval; } \
	1198	else if (val != (type)val) { bRounded = TRUE; val = (type)(val + 0.5); } } \
	1199	while(0)
	1200
	1201	switch(GDALGetRasterDataType(hBand))
	1202	{
	1203	case GDT_Byte:
	1204	CLAMP(psWO->padfDstNoDataReal[i], GByte,
	1205	0.0, 255.0);
	1206	break;
	1207	case GDT_Int16:
	1208	CLAMP(psWO->padfDstNoDataReal[i], GInt16,
	1209	-32768.0, 32767.0);
	1210	break;
	1211	case GDT_UInt16:
	1212	CLAMP(psWO->padfDstNoDataReal[i], GUInt16,
	1213	0.0, 65535.0);
	1214	break;
	1215	case GDT_Int32:
	1216	CLAMP(psWO->padfDstNoDataReal[i], GInt32,
	1217	-2147483648.0, 2147483647.0);
	1218	break;
	1219	case GDT_UInt32:
	1220	CLAMP(psWO->padfDstNoDataReal[i], GUInt32,
	1221	0.0, 4294967295.0);
	1222	break;
	1223	default:
	1224	break;
	1225	}
	1226
	1227	if (bClamped)
	1228	{
	1229	printf( "for band %d, destination nodata value has been clamped "
	1230	"to %.0f, the original value being out of range.\n",
	1231	i + 1, psWO->padfDstNoDataReal[i]);
	1232	}
	1233	else if(bRounded)
	1234	{
	1235	printf("for band %d, destination nodata value has been rounded "
	1236	"to %.0f, %s being an integer datatype.\n",
	1237	i + 1, psWO->padfDstNoDataReal[i],
	1238	GDALGetDataTypeName(GDALGetRasterDataType(hBand)));
	1239	}
	1240
	1241	if( bCreateOutput )
	1242	{
	1243	GDALSetRasterNoDataValue(
	1244	GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ),
	1245	psWO->padfDstNoDataReal[i] );
	1246	}
	1247	}
	1248
	1249	CSLDestroy( papszTokens );
	1250	}
	1251
	1252	/* -------------------------------------------------------------------- */
	1253	/* If we have a cutline, transform it into the source */
	1254	/* pixel/line coordinate system and insert into warp options. */
	1255	/* -------------------------------------------------------------------- */
	1256	if( hCutline != NULL )
	1257	{
	1258	TransformCutlineToSource( hSrcDS, hCutline,
	1259	&(psWO->papszWarpOptions),
	1260	papszTO );
	1261	}
	1262
	1263	/* -------------------------------------------------------------------- */
	1264	/* If we are producing VRT output, then just initialize it with */
	1265	/* the warp options and write out now rather than proceeding */
	1266	/* with the operations. */
	1267	/* -------------------------------------------------------------------- */
	1268	if( bVRT )
	1269	{
	1270	if( GDALInitializeWarpedVRT( hDstDS, psWO ) != CE_None )
	1271	GDALExit( 1 );
	1272
	1273	GDALClose( hDstDS );
	1274	GDALClose( hSrcDS );
	1275
	1276	/* The warped VRT will clean itself the transformer used */
	1277	/* So we have only to destroy the hGenImgProjArg if we */
	1278	/* have wrapped it inside the hApproxArg */
	1279	if (pfnTransformer == GDALApproxTransform)
	1280	{
	1281	if( hGenImgProjArg != NULL )
	1282	GDALDestroyGenImgProjTransformer( hGenImgProjArg );
	1283	}
	1284
	1285	GDALDestroyWarpOptions( psWO );
	1286
	1287	setMapInMaps(outputs,"Result","value",pszDstFilename);
	1288
	1289	CPLFree( pszDstFilename );
	1290	CSLDestroy( papszSrcFiles );
	1291	CSLDestroy( papszWarpOptions );
	1292	CSLDestroy( papszTO );
	1293
	1294	GDALDumpOpenDatasets( stderr );
	1295
	1296	GDALDestroyDriverManager();
	1297
	1298	return SERVICE_SUCCEEDED;
	1299	}
	1300
	1301	/* -------------------------------------------------------------------- */
	1302	/* Initialize and execute the warp. */
	1303	/* -------------------------------------------------------------------- */
	1304	GDALWarpOperation oWO;
	1305
	1306	if( oWO.Initialize( psWO ) == CE_None )
	1307	{
	1308	CPLErr eErr;
	1309	if( bMulti )
	1310	eErr = oWO.ChunkAndWarpMulti( 0, 0,
	1311	GDALGetRasterXSize( hDstDS ),
	1312	GDALGetRasterYSize( hDstDS ) );
	1313	else
	1314	eErr = oWO.ChunkAndWarpImage( 0, 0,
	1315	GDALGetRasterXSize( hDstDS ),
	1316	GDALGetRasterYSize( hDstDS ) );
	1317	if (eErr != CE_None){
	1318	bHasGotErr = TRUE;
	1319	}
	1320	}
	1321
	1322	/* -------------------------------------------------------------------- */
	1323	/* Cleanup */
	1324	/* -------------------------------------------------------------------- */
	1325	if( hApproxArg != NULL )
	1326	GDALDestroyApproxTransformer( hApproxArg );
	1327
	1328	if( hGenImgProjArg != NULL )
	1329	GDALDestroyGenImgProjTransformer( hGenImgProjArg );
	1330
	1331	GDALDestroyWarpOptions( psWO );
	1332
	1333	GDALClose( hSrcDS );
	1334	}
	1335
	1336	/* -------------------------------------------------------------------- */
	1337	/* Final Cleanup. */
	1338	/* -------------------------------------------------------------------- */
	1339	CPLErrorReset();
	1340	GDALFlushCache( hDstDS );
	1341	if( CPLGetLastErrorType() != CE_None )
	1342	bHasGotErr = TRUE;
	1343	GDALClose( hDstDS );
	1344
	1345	setMapInMaps(outputs,"Result","value",pszDstFilename);
	1346
	1347	CPLFree( pszDstFilename );
	1348	CSLDestroy( papszSrcFiles );
	1349	CSLDestroy( papszWarpOptions );
	1350	CSLDestroy( papszTO );
	1351
	1352	GDALDumpOpenDatasets( stderr );
	1353
	1354	GDALDestroyDriverManager();
	1355
	1356	#ifdef OGR_ENABLED
	1357	if( hCutline != NULL )
	1358	OGR_G_DestroyGeometry( (OGRGeometryH) hCutline );
	1359	OGRCleanupAll();
	1360	#endif
	1361
	1362
	1363	/* -------------------------------------------------------------------- */
	1364	/* Clean allocated arguments. */
	1365	/* -------------------------------------------------------------------- */
	1366	free(dataPath);
	1367	free(tempPath);
	1368
	1369	tmpMap=getMapFromMaps(inputs,"te","value");
	1370	if(tmpMap!=NULL){
	1371	free(pszCutlineDSName);
	1372	}
	1373
	1374	tmpMap=getMapFromMaps(inputs,"cwhere","value");
	1375	if(tmpMap!=NULL){
	1376	free(pszCWHERE);
	1377	}
	1378
	1379	tmpMap=getMapFromMaps(inputs,"cl","value");
	1380	if(tmpMap!=NULL){
	1381	free(pszCLayer);
	1382	}
	1383
	1384	tmpMap=getMapFromMaps(inputs,"csql","value");
	1385	if(tmpMap!=NULL){
	1386	free(pszCSQL);
	1387	}
	1388
	1389	return (bHasGotErr) ? SERVICE_FAILED : SERVICE_SUCCEEDED;
	1390	}
	1391
	1392	/************************************************************************/
	1393	/* GDALWarpCreateOutput() */
	1394	/* */
	1395	/* Create the output file based on various commandline options, */
	1396	/* and the input file. */
	1397	/* If there's just one source file, then phTransformArg and /
	1398	/* phSrcDS will be set, in order them to be reused by main /
	1399	/* function. This saves dataset re-opening, and above all transform*/
	1400	/* recomputation, which can be expensive in the -tps case */
	1401	/************************************************************************/
	1402
	1403	static GDALDatasetH
	1404	GDALWarpCreateOutput( maps& conf,char papszSrcFiles, const char pszFilename,
	1405	const char pszFormat, char *papszTO,
	1406	char ***ppapszCreateOptions, GDALDataType eDT,
	1407	void ** phTransformArg,
	1408	GDALDatasetH* phSrcDS)
	1409
	1410
	1411	{
	1412	GDALDriverH hDriver;
	1413	GDALDatasetH hDstDS;
	1414	void *hTransformArg;
	1415	GDALColorTableH hCT = NULL;
	1416	double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0;
	1417	double dfWrkResX=0, dfWrkResY=0;
	1418	int nDstBandCount = 0;
	1419	std::vector<GDALColorInterp> apeColorInterpretations;
	1420
	1421	*phTransformArg = NULL;
	1422	*phSrcDS = NULL;
	1423
	1424	/* -------------------------------------------------------------------- */
	1425	/* Find the output driver. */
	1426	/* -------------------------------------------------------------------- */
	1427	hDriver = GDALGetDriverByName( pszFormat );
	1428	if( hDriver == NULL
	1429	\|\| GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
	1430	{
	1431	int iDr;
	1432	char tmp[2048];
	1433	sprintf( tmp, "Output driver `%s' not recognised or does not support\n"
	1434	"direct output file creation. The following format drivers are configured\n"
	1435	"and support direct output:\n",
	1436	pszFormat );
	1437
	1438	for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
	1439	{
	1440	GDALDriverH hDriver = GDALGetDriver(iDr);
	1441
	1442	if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
	1443	{
	1444	char *_tmp=strdup(tmp);
	1445	sprintf( tmp,"%s %s: %s\n",
	1446	_tmp,
	1447	GDALGetDriverShortName( hDriver ),
	1448	GDALGetDriverLongName( hDriver ) );
	1449	free(_tmp);
	1450	}
	1451	}
	1452	setMapInMaps(conf,"lenv","message",tmp);
	1453	GDALExit( 1 );
	1454	}
	1455
	1456	/* -------------------------------------------------------------------- */
	1457	/* For virtual output files, we have to set a special subclass */
	1458	/* of dataset to create. */
	1459	/* -------------------------------------------------------------------- */
	1460	if( bVRT )
	1461	*ppapszCreateOptions =
	1462	CSLSetNameValue( *ppapszCreateOptions, "SUBCLASS",
	1463	"VRTWarpedDataset" );
	1464
	1465	/* -------------------------------------------------------------------- */
	1466	/* Loop over all input files to collect extents. */
	1467	/* -------------------------------------------------------------------- */
	1468	int iSrc;
	1469	char pszThisTargetSRS = (char)CSLFetchNameValue( papszTO, "DST_SRS" );
	1470	if( pszThisTargetSRS != NULL )
	1471	pszThisTargetSRS = CPLStrdup( pszThisTargetSRS );
	1472
	1473	for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
	1474	{
	1475	GDALDatasetH hSrcDS;
	1476	const char *pszThisSourceSRS = CSLFetchNameValue(papszTO,"SRC_SRS");
	1477
	1478	hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
	1479	if( hSrcDS == NULL )
	1480	GDALExit( 1 );
	1481
	1482	/* -------------------------------------------------------------------- */
	1483	/* Check that there's at least one raster band */
	1484	/* -------------------------------------------------------------------- */
	1485	if ( GDALGetRasterCount(hSrcDS) == 0 )
	1486	{
	1487	fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
	1488	GDALExit( 1 );
	1489	}
	1490
	1491	if( eDT == GDT_Unknown )
	1492	eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
	1493
	1494	/* -------------------------------------------------------------------- */
	1495	/* If we are processing the first file, and it has a color */
	1496	/* table, then we will copy it to the destination file. */
	1497	/* -------------------------------------------------------------------- */
	1498	if( iSrc == 0 )
	1499	{
	1500	nDstBandCount = GDALGetRasterCount(hSrcDS);
	1501	hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
	1502	if( hCT != NULL )
	1503	{
	1504	hCT = GDALCloneColorTable( hCT );
	1505	if( !bQuiet )
	1506	printf( "Copying color table from %s to new file.\n",
	1507	papszSrcFiles[iSrc] );
	1508	}
	1509
	1510	for(int iBand = 0; iBand < nDstBandCount; iBand++)
	1511	{
	1512	apeColorInterpretations.push_back(
	1513	GDALGetRasterColorInterpretation(GDALGetRasterBand(hSrcDS,iBand+1)) );
	1514	}
	1515	}
	1516
	1517	/* -------------------------------------------------------------------- */
	1518	/* Get the sourcesrs from the dataset, if not set already. */
	1519	/* -------------------------------------------------------------------- */
	1520	if( pszThisSourceSRS == NULL )
	1521	{
	1522	const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" );
	1523
	1524	if( GDALGetProjectionRef( hSrcDS ) != NULL
	1525	&& strlen(GDALGetProjectionRef( hSrcDS )) > 0
	1526	&& (pszMethod == NULL \|\| EQUAL(pszMethod,"GEOTRANSFORM")) )
	1527	pszThisSourceSRS = GDALGetProjectionRef( hSrcDS );
	1528
	1529	else if( GDALGetGCPProjection( hSrcDS ) != NULL
	1530	&& strlen(GDALGetGCPProjection(hSrcDS)) > 0
	1531	&& GDALGetGCPCount( hSrcDS ) > 1
	1532	&& (pszMethod == NULL \|\| EQUALN(pszMethod,"GCP_",4)) )
	1533	pszThisSourceSRS = GDALGetGCPProjection( hSrcDS );
	1534	else if( pszMethod != NULL && EQUAL(pszMethod,"RPC") )
	1535	pszThisSourceSRS = SRS_WKT_WGS84;
	1536	else
	1537	pszThisSourceSRS = "";
	1538	}
	1539
	1540	if( pszThisTargetSRS == NULL )
	1541	pszThisTargetSRS = CPLStrdup( pszThisSourceSRS );
	1542
	1543	/* -------------------------------------------------------------------- */
	1544	/* Create a transformation object from the source to */
	1545	/* destination coordinate system. */
	1546	/* -------------------------------------------------------------------- */
	1547	hTransformArg =
	1548	GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
	1549
	1550	if( hTransformArg == NULL )
	1551	{
	1552	CPLFree( pszThisTargetSRS );
	1553	GDALClose( hSrcDS );
	1554	return NULL;
	1555	}
	1556
	1557	GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;
	1558
	1559	/* -------------------------------------------------------------------- */
	1560	/* Get approximate output definition. */
	1561	/* -------------------------------------------------------------------- */
	1562	double adfThisGeoTransform[6];
	1563	double adfExtent[4];
	1564	int nThisPixels, nThisLines;
	1565
	1566	if( GDALSuggestedWarpOutput2( hSrcDS,
	1567	psInfo->pfnTransform, hTransformArg,
	1568	adfThisGeoTransform,
	1569	&nThisPixels, &nThisLines,
	1570	adfExtent, 0 ) != CE_None )
	1571	{
	1572	CPLFree( pszThisTargetSRS );
	1573	GDALClose( hSrcDS );
	1574	return NULL;
	1575	}
	1576
	1577	if (CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", NULL ) == NULL)
	1578	{
	1579	double MinX = adfExtent[0];
	1580	double MaxX = adfExtent[2];
	1581	double MaxY = adfExtent[3];
	1582	double MinY = adfExtent[1];
	1583	int bSuccess = TRUE;
	1584
	1585	/* Check that the the edges of the target image are in the validity area */
	1586	/* of the target projection */
	1587	#define N_STEPS 20
	1588	int i,j;
	1589	for(i=0;i<=N_STEPS && bSuccess;i++)
	1590	{
	1591	for(j=0;j<=N_STEPS && bSuccess;j++)
	1592	{
	1593	double dfRatioI = i * 1.0 / N_STEPS;
	1594	double dfRatioJ = j * 1.0 / N_STEPS;
	1595	double expected_x = (1 - dfRatioI) * MinX + dfRatioI * MaxX;
	1596	double expected_y = (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
	1597	double x = expected_x;
	1598	double y = expected_y;
	1599	double z = 0;
	1600	/* Target SRS coordinates to source image pixel coordinates */
	1601	if (!psInfo->pfnTransform(hTransformArg, TRUE, 1, &x, &y, &z, &bSuccess) \|\| !bSuccess)
	1602	bSuccess = FALSE;
	1603	/* Source image pixel coordinates to target SRS coordinates */
	1604	if (!psInfo->pfnTransform(hTransformArg, FALSE, 1, &x, &y, &z, &bSuccess) \|\| !bSuccess)
	1605	bSuccess = FALSE;
	1606	if (fabs(x - expected_x) > (MaxX - MinX) / nThisPixels \|\|
	1607	fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
	1608	bSuccess = FALSE;
	1609	}
	1610	}
	1611
	1612	/* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces ogrct.cpp */
	1613	/* to check the consistency of each requested projection result with the */
	1614	/* invert projection */
	1615	if (!bSuccess)
	1616	{
	1617	CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" );
	1618	CPLDebug("WARP", "Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE");
	1619
	1620	if( GDALSuggestedWarpOutput2( hSrcDS,
	1621	psInfo->pfnTransform, hTransformArg,
	1622	adfThisGeoTransform,
	1623	&nThisPixels, &nThisLines,
	1624	adfExtent, 0 ) != CE_None )
	1625	{
	1626	CPLFree( pszThisTargetSRS );
	1627	GDALClose( hSrcDS );
	1628	return NULL;
	1629	}
	1630	}
	1631	}
	1632
	1633	/* -------------------------------------------------------------------- */
	1634	/* Expand the working bounds to include this region, ensure the */
	1635	/* working resolution is no more than this resolution. */
	1636	/* -------------------------------------------------------------------- */
	1637	if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 )
	1638	{
	1639	dfWrkMinX = adfExtent[0];
	1640	dfWrkMaxX = adfExtent[2];
	1641	dfWrkMaxY = adfExtent[3];
	1642	dfWrkMinY = adfExtent[1];
	1643	dfWrkResX = adfThisGeoTransform[1];
	1644	dfWrkResY = ABS(adfThisGeoTransform[5]);
	1645	}
	1646	else
	1647	{
	1648	dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]);
	1649	dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]);
	1650	dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]);
	1651	dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]);
	1652	dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]);
	1653	dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5]));
	1654	}
	1655
	1656	if (iSrc == 0 && papszSrcFiles[1] == NULL)
	1657	{
	1658	*phTransformArg = hTransformArg;
	1659	*phSrcDS = hSrcDS;
	1660	}
	1661	else
	1662	{
	1663	GDALDestroyGenImgProjTransformer( hTransformArg );
	1664	GDALClose( hSrcDS );
	1665	}
	1666	}
	1667
	1668	/* -------------------------------------------------------------------- */
	1669	/* Did we have any usable sources? */
	1670	/* -------------------------------------------------------------------- */
	1671	if( nDstBandCount == 0 )
	1672	{
	1673	CPLError( CE_Failure, CPLE_AppDefined,
	1674	"No usable source images." );
	1675	CPLFree( pszThisTargetSRS );
	1676	return NULL;
	1677	}
	1678
	1679	/* -------------------------------------------------------------------- */
	1680	/* Turn the suggested region into a geotransform and suggested */
	1681	/* number of pixels and lines. */
	1682	/* -------------------------------------------------------------------- */
	1683	double adfDstGeoTransform[6];
	1684	int nPixels, nLines;
	1685
	1686	adfDstGeoTransform[0] = dfWrkMinX;
	1687	adfDstGeoTransform[1] = dfWrkResX;
	1688	adfDstGeoTransform[2] = 0.0;
	1689	adfDstGeoTransform[3] = dfWrkMaxY;
	1690	adfDstGeoTransform[4] = 0.0;
	1691	adfDstGeoTransform[5] = -1 * dfWrkResY;
	1692
	1693	nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
	1694	nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
	1695
	1696	/* -------------------------------------------------------------------- */
	1697	/* Did the user override some parameters? */
	1698	/* -------------------------------------------------------------------- */
	1699	if( dfXRes != 0.0 && dfYRes != 0.0 )
	1700	{
	1701	if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
	1702	{
	1703	dfMinX = adfDstGeoTransform[0];
	1704	dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
	1705	dfMaxY = adfDstGeoTransform[3];
	1706	dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
	1707	}
	1708
	1709	if ( bTargetAlignedPixels )
	1710	{
	1711	dfMinX = floor(dfMinX / dfXRes) * dfXRes;
	1712	dfMaxX = ceil(dfMaxX / dfXRes) * dfXRes;
	1713	dfMinY = floor(dfMinY / dfYRes) * dfYRes;
	1714	dfMaxY = ceil(dfMaxY / dfYRes) * dfYRes;
	1715	}
	1716
	1717	nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
	1718	nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
	1719	adfDstGeoTransform[0] = dfMinX;
	1720	adfDstGeoTransform[3] = dfMaxY;
	1721	adfDstGeoTransform[1] = dfXRes;
	1722	adfDstGeoTransform[5] = -dfYRes;
	1723	}
	1724
	1725	else if( nForcePixels != 0 && nForceLines != 0 )
	1726	{
	1727	if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
	1728	{
	1729	dfMinX = dfWrkMinX;
	1730	dfMaxX = dfWrkMaxX;
	1731	dfMaxY = dfWrkMaxY;
	1732	dfMinY = dfWrkMinY;
	1733	}
	1734
	1735	dfXRes = (dfMaxX - dfMinX) / nForcePixels;
	1736	dfYRes = (dfMaxY - dfMinY) / nForceLines;
	1737
	1738	adfDstGeoTransform[0] = dfMinX;
	1739	adfDstGeoTransform[3] = dfMaxY;
	1740	adfDstGeoTransform[1] = dfXRes;
	1741	adfDstGeoTransform[5] = -dfYRes;
	1742
	1743	nPixels = nForcePixels;
	1744	nLines = nForceLines;
	1745	}
	1746
	1747	else if( nForcePixels != 0 )
	1748	{
	1749	if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
	1750	{
	1751	dfMinX = dfWrkMinX;
	1752	dfMaxX = dfWrkMaxX;
	1753	dfMaxY = dfWrkMaxY;
	1754	dfMinY = dfWrkMinY;
	1755	}
	1756
	1757	dfXRes = (dfMaxX - dfMinX) / nForcePixels;
	1758	dfYRes = dfXRes;
	1759
	1760	adfDstGeoTransform[0] = dfMinX;
	1761	adfDstGeoTransform[3] = dfMaxY;
	1762	adfDstGeoTransform[1] = dfXRes;
	1763	adfDstGeoTransform[5] = -dfYRes;
	1764
	1765	nPixels = nForcePixels;
	1766	nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
	1767	}
	1768
	1769	else if( nForceLines != 0 )
	1770	{
	1771	if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
	1772	{
	1773	dfMinX = dfWrkMinX;
	1774	dfMaxX = dfWrkMaxX;
	1775	dfMaxY = dfWrkMaxY;
	1776	dfMinY = dfWrkMinY;
	1777	}
	1778
	1779	dfYRes = (dfMaxY - dfMinY) / nForceLines;
	1780	dfXRes = dfYRes;
	1781
	1782	adfDstGeoTransform[0] = dfMinX;
	1783	adfDstGeoTransform[3] = dfMaxY;
	1784	adfDstGeoTransform[1] = dfXRes;
	1785	adfDstGeoTransform[5] = -dfYRes;
	1786
	1787	nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
	1788	nLines = nForceLines;
	1789	}
	1790
	1791	else if( dfMinX != 0.0 \|\| dfMinY != 0.0 \|\| dfMaxX != 0.0 \|\| dfMaxY != 0.0 )
	1792	{
	1793	dfXRes = adfDstGeoTransform[1];
	1794	dfYRes = fabs(adfDstGeoTransform[5]);
	1795
	1796	nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
	1797	nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
	1798
	1799	dfXRes = (dfMaxX - dfMinX) / nPixels;
	1800	dfYRes = (dfMaxY - dfMinY) / nLines;
	1801
	1802	adfDstGeoTransform[0] = dfMinX;
	1803	adfDstGeoTransform[3] = dfMaxY;
	1804	adfDstGeoTransform[1] = dfXRes;
	1805	adfDstGeoTransform[5] = -dfYRes;
	1806	}
	1807
	1808	/* -------------------------------------------------------------------- */
	1809	/* Do we want to generate an alpha band in the output file? */
	1810	/* -------------------------------------------------------------------- */
	1811	if( bEnableSrcAlpha )
	1812	nDstBandCount--;
	1813
	1814	if( bEnableDstAlpha )
	1815	nDstBandCount++;
	1816
	1817	/* -------------------------------------------------------------------- */
	1818	/* Create the output file. */
	1819	/* -------------------------------------------------------------------- */
	1820	if( !bQuiet )
	1821	printf( "Creating output file that is %dP x %dL.\n", nPixels, nLines );
	1822
	1823	hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
	1824	nDstBandCount, eDT, *ppapszCreateOptions );
	1825
	1826	if( hDstDS == NULL )
	1827	{
	1828	CPLFree( pszThisTargetSRS );
	1829	return NULL;
	1830	}
	1831
	1832	/* -------------------------------------------------------------------- */
	1833	/* Write out the projection definition. */
	1834	/* -------------------------------------------------------------------- */
	1835	GDALSetProjection( hDstDS, pszThisTargetSRS );
	1836	GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
	1837
	1838	if (*phTransformArg != NULL)
	1839	GDALSetGenImgProjTransformerDstGeoTransform( *phTransformArg, adfDstGeoTransform);
	1840
	1841	/* -------------------------------------------------------------------- */
	1842	/* Try to set color interpretation of source bands to target */
	1843	/* dataset. */
	1844	/* FIXME? We should likely do that for other drivers than VRT */
	1845	/* but it might create spurious .aux.xml files (at least with HFA, */
	1846	/* and netCDF) */
	1847	/* -------------------------------------------------------------------- */
	1848	if( bVRT )
	1849	{
	1850	int nBandsToCopy = (int)apeColorInterpretations.size();
	1851	if ( bEnableSrcAlpha )
	1852	nBandsToCopy --;
	1853	for(int iBand = 0; iBand < nBandsToCopy; iBand++)
	1854	{
	1855	GDALSetRasterColorInterpretation(
	1856	GDALGetRasterBand( hDstDS, iBand + 1 ),
	1857	apeColorInterpretations[iBand] );
	1858	}
	1859	}
	1860
	1861	/* -------------------------------------------------------------------- */
	1862	/* Try to set color interpretation of output file alpha band. */
	1863	/* -------------------------------------------------------------------- */
	1864	if( bEnableDstAlpha )
	1865	{
	1866	GDALSetRasterColorInterpretation(
	1867	GDALGetRasterBand( hDstDS, nDstBandCount ),
	1868	GCI_AlphaBand );
	1869	}
	1870
	1871	/* -------------------------------------------------------------------- */
	1872	/* Copy the color table, if required. */
	1873	/* -------------------------------------------------------------------- */
	1874	if( hCT != NULL )
	1875	{
	1876	GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
	1877	GDALDestroyColorTable( hCT );
	1878	}
	1879
	1880	CPLFree( pszThisTargetSRS );
	1881	return hDstDS;
	1882	}
	1883
	1884	/************************************************************************/
	1885	/* GeoTransform_Transformer() */
	1886	/* */
	1887	/* Convert points from georef coordinates to pixel/line based */
	1888	/* on a geotransform. */
	1889	/************************************************************************/
	1890
	1891	class CutlineTransformer : public OGRCoordinateTransformation
	1892	{
	1893	public:
	1894
	1895	void *hSrcImageTransformer;
	1896
	1897	virtual OGRSpatialReference *GetSourceCS() { return NULL; }
	1898	virtual OGRSpatialReference *GetTargetCS() { return NULL; }
	1899
	1900	virtual int Transform( int nCount,
	1901	double x, double y, double *z = NULL ) {
	1902	int nResult;
	1903
	1904	int pabSuccess = (int ) CPLCalloc(sizeof(int),nCount);
	1905	nResult = TransformEx( nCount, x, y, z, pabSuccess );
	1906	CPLFree( pabSuccess );
	1907
	1908	return nResult;
	1909	}
	1910
	1911	virtual int TransformEx( int nCount,
	1912	double x, double y, double *z = NULL,
	1913	int *pabSuccess = NULL ) {
	1914	return GDALGenImgProjTransform( hSrcImageTransformer, TRUE,
	1915	nCount, x, y, z, pabSuccess );
	1916	}
	1917	};
	1918
	1919
	1920	/************************************************************************/
	1921	/* LoadCutline() */
	1922	/* */
	1923	/* Load blend cutline from OGR datasource. */
	1924	/************************************************************************/
	1925
	1926	static void
	1927	LoadCutline( const char pszCutlineDSName, const char pszCLayer,
	1928	const char pszCWHERE, const char pszCSQL,
	1929	void **phCutlineRet )
	1930
	1931	{
	1932	#ifndef OGR_ENABLED
	1933	CPLError( CE_Failure, CPLE_AppDefined,
	1934	"Request to load a cutline failed, this build does not support OGR features.\n" );
	1935	GDALExit( 1 );
	1936	#else // def OGR_ENABLED
	1937	OGRRegisterAll();
	1938
	1939	/* -------------------------------------------------------------------- */
	1940	/* Open source vector dataset. */
	1941	/* -------------------------------------------------------------------- */
	1942	OGRDataSourceH hSrcDS;
	1943
	1944	hSrcDS = OGROpen( pszCutlineDSName, FALSE, NULL );
	1945	if( hSrcDS == NULL )
	1946	GDALExit( 1 );
	1947
	1948	/* -------------------------------------------------------------------- */
	1949	/* Get the source layer */
	1950	/* -------------------------------------------------------------------- */
	1951	OGRLayerH hLayer = NULL;
	1952
	1953	if( pszCSQL != NULL )
	1954	hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszCSQL, NULL, NULL );
	1955	else if( pszCLayer != NULL )
	1956	hLayer = OGR_DS_GetLayerByName( hSrcDS, pszCLayer );
	1957	else
	1958	hLayer = OGR_DS_GetLayer( hSrcDS, 0 );
	1959
	1960	if( hLayer == NULL )
	1961	{
	1962	fprintf( stderr, "Failed to identify source layer from datasource.\n" );
	1963	GDALExit( 1 );
	1964	}
	1965
	1966	/* -------------------------------------------------------------------- */
	1967	/* Apply WHERE clause if there is one. */
	1968	/* -------------------------------------------------------------------- */
	1969	if( pszCWHERE != NULL )
	1970	OGR_L_SetAttributeFilter( hLayer, pszCWHERE );
	1971
	1972	/* -------------------------------------------------------------------- */
	1973	/* Collect the geometries from this layer, and build list of */
	1974	/* burn values. */
	1975	/* -------------------------------------------------------------------- */
	1976	OGRFeatureH hFeat;
	1977	OGRGeometryH hMultiPolygon = OGR_G_CreateGeometry( wkbMultiPolygon );
	1978
	1979	OGR_L_ResetReading( hLayer );
	1980
	1981	while( (hFeat = OGR_L_GetNextFeature( hLayer )) != NULL )
	1982	{
	1983	OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat);
	1984
	1985	if( hGeom == NULL )
	1986	{
	1987	fprintf( stderr, "ERROR: Cutline feature without a geometry.\n" );
	1988	GDALExit( 1 );
	1989	}
	1990
	1991	OGRwkbGeometryType eType = wkbFlatten(OGR_G_GetGeometryType( hGeom ));
	1992
	1993	if( eType == wkbPolygon )
	1994	OGR_G_AddGeometry( hMultiPolygon, hGeom );
	1995	else if( eType == wkbMultiPolygon )
	1996	{
	1997	int iGeom;
	1998
	1999	for( iGeom = 0; iGeom < OGR_G_GetGeometryCount( hGeom ); iGeom++ )
	2000	{
	2001	OGR_G_AddGeometry( hMultiPolygon,
	2002	OGR_G_GetGeometryRef(hGeom,iGeom) );
	2003	}
	2004	}
	2005	else
	2006	{
	2007	fprintf( stderr, "ERROR: Cutline not of polygon type.\n" );
	2008	GDALExit( 1 );
	2009	}
	2010
	2011	OGR_F_Destroy( hFeat );
	2012	}
	2013
	2014	if( OGR_G_GetGeometryCount( hMultiPolygon ) == 0 )
	2015	{
	2016	fprintf( stderr, "ERROR: Did not get any cutline features.\n" );
	2017	GDALExit( 1 );
	2018	}
	2019
	2020	/* -------------------------------------------------------------------- */
	2021	/* Ensure the coordinate system gets set on the geometry. */
	2022	/* -------------------------------------------------------------------- */
	2023	OGR_G_AssignSpatialReference(
	2024	hMultiPolygon, OGR_L_GetSpatialRef(hLayer) );
	2025
	2026	phCutlineRet = (void ) hMultiPolygon;
	2027
	2028	/* -------------------------------------------------------------------- */
	2029	/* Cleanup */
	2030	/* -------------------------------------------------------------------- */
	2031	if( pszCSQL != NULL )
	2032	OGR_DS_ReleaseResultSet( hSrcDS, hLayer );
	2033
	2034	OGR_DS_Destroy( hSrcDS );
	2035	#endif
	2036	}
	2037
	2038	/************************************************************************/
	2039	/* TransformCutlineToSource() */
	2040	/* */
	2041	/* Transform cutline from its SRS to source pixel/line coordinates.*/
	2042	/************************************************************************/
	2043	static void
	2044	TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
	2045	char *ppapszWarpOptions, char papszTO_In )
	2046
	2047	{
	2048	#ifdef OGR_ENABLED
	2049	OGRGeometryH hMultiPolygon = OGR_G_Clone( (OGRGeometryH) hCutline );
	2050	char **papszTO = CSLDuplicate( papszTO_In );
	2051
	2052	/* -------------------------------------------------------------------- */
	2053	/* Checkout that SRS are the same. */
	2054	/* -------------------------------------------------------------------- */
	2055	OGRSpatialReferenceH hRasterSRS = NULL;
	2056	const char *pszProjection = NULL;
	2057
	2058	if( GDALGetProjectionRef( hSrcDS ) != NULL
	2059	&& strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
	2060	pszProjection = GDALGetProjectionRef( hSrcDS );
	2061	else if( GDALGetGCPProjection( hSrcDS ) != NULL )
	2062	pszProjection = GDALGetGCPProjection( hSrcDS );
	2063
	2064	if( pszProjection != NULL )
	2065	{
	2066	hRasterSRS = OSRNewSpatialReference(NULL);
	2067	if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None )
	2068	{
	2069	OSRDestroySpatialReference(hRasterSRS);
	2070	hRasterSRS = NULL;
	2071	}
	2072	}
	2073
	2074	OGRSpatialReferenceH hCutlineSRS = OGR_G_GetSpatialReference( hMultiPolygon );
	2075	if( hRasterSRS != NULL && hCutlineSRS != NULL )
	2076	{
	2077	/* ok, we will reproject */
	2078	}
	2079	else if( hRasterSRS != NULL && hCutlineSRS == NULL )
	2080	{
	2081	fprintf(stderr,
	2082	"Warning : the source raster dataset has a SRS, but the cutline features\n"
	2083	"not. We assume that the cutline coordinates are expressed in the destination SRS.\n"
	2084	"If not, cutline results may be incorrect.\n");
	2085	}
	2086	else if( hRasterSRS == NULL && hCutlineSRS != NULL )
	2087	{
	2088	fprintf(stderr,
	2089	"Warning : the input vector layer has a SRS, but the source raster dataset does not.\n"
	2090	"Cutline results may be incorrect.\n");
	2091	}
	2092
	2093	if( hRasterSRS != NULL )
	2094	OSRDestroySpatialReference(hRasterSRS);
	2095
	2096	/* -------------------------------------------------------------------- */
	2097	/* Extract the cutline SRS WKT. */
	2098	/* -------------------------------------------------------------------- */
	2099	if( hCutlineSRS != NULL )
	2100	{
	2101	char *pszCutlineSRS_WKT = NULL;
	2102
	2103	OSRExportToWkt( hCutlineSRS, &pszCutlineSRS_WKT );
	2104	papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszCutlineSRS_WKT );
	2105	CPLFree( pszCutlineSRS_WKT );
	2106	}
	2107
	2108	/* -------------------------------------------------------------------- */
	2109	/* It may be unwise to let the mask geometry be re-wrapped by */
	2110	/* the CENTER_LONG machinery as this can easily screw up world */
	2111	/* spanning masks and invert the mask topology. */
	2112	/* -------------------------------------------------------------------- */
	2113	papszTO = CSLSetNameValue( papszTO, "INSERT_CENTER_LONG", "FALSE" );
	2114
	2115	/* -------------------------------------------------------------------- */
	2116	/* Transform the geometry to pixel/line coordinates. */
	2117	/* -------------------------------------------------------------------- */
	2118	CutlineTransformer oTransformer;
	2119
	2120	/* The cutline transformer will invert the hSrcImageTransformer */
	2121	/* so it will convert from the cutline SRS to the source pixel/line */
	2122	/* coordinates */
	2123	oTransformer.hSrcImageTransformer =
	2124	GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
	2125
	2126	CSLDestroy( papszTO );
	2127
	2128	if( oTransformer.hSrcImageTransformer == NULL )
	2129	GDALExit( 1 );
	2130
	2131	OGR_G_Transform( hMultiPolygon,
	2132	(OGRCoordinateTransformationH) &oTransformer );
	2133
	2134	GDALDestroyGenImgProjTransformer( oTransformer.hSrcImageTransformer );
	2135
	2136	/* -------------------------------------------------------------------- */
	2137	/* Convert aggregate geometry into WKT. */
	2138	/* -------------------------------------------------------------------- */
	2139	char *pszWKT = NULL;
	2140
	2141	OGR_G_ExportToWkt( hMultiPolygon, &pszWKT );
	2142	OGR_G_DestroyGeometry( hMultiPolygon );
	2143
	2144	ppapszWarpOptions = CSLSetNameValue( ppapszWarpOptions,
	2145	"CUTLINE", pszWKT );
	2146	CPLFree( pszWKT );
	2147	#endif
	2148	}
	2149
	2150	}

Note: See TracBrowser for help on using the repository browser.

Open WPS platform

source: trunk/zoo-project/zoo-services/gdal/warp/service.c @ 940

Download in other formats:

Download now !

Support ZOO-Project

Search

Context Navigation

Main Navigation

ZOO Sponsors

Knowledge partners

Related links