service.c @ 920

Last change on this file since 920 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

Line
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 @ 920

Download in other formats:

Download now !

Support ZOO-Project

Search

Context Navigation

Main Navigation

ZOO Sponsors

Knowledge partners

Related links