1 : /******************************************************************************
2 : * $Id: gdalwarp.cpp 19692 2010-05-13 17:16:55Z rouault $
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 :
36 : CPL_CVSID("$Id: gdalwarp.cpp 19692 2010-05-13 17:16:55Z rouault $");
37 :
38 : static void
39 : LoadCutline( const char *pszCutlineDSName, const char *pszCLayer,
40 : const char *pszCWHERE, const char *pszCSQL,
41 : void **phCutlineRet );
42 : static void
43 : TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
44 : char ***ppapszWarpOptions, char **papszTO );
45 : static GDALDatasetH
46 : GDALWarpCreateOutput( char **papszSrcFiles, const char *pszFilename,
47 : const char *pszFormat, char **papszTO,
48 : char ***ppapszCreateOptions, GDALDataType eDT );
49 :
50 : static double dfMinX=0.0, dfMinY=0.0, dfMaxX=0.0, dfMaxY=0.0;
51 : static double dfXRes=0.0, dfYRes=0.0;
52 : static int nForcePixels=0, nForceLines=0, bQuiet = FALSE;
53 : static int bEnableDstAlpha = FALSE, bEnableSrcAlpha = FALSE;
54 :
55 : static int bVRT = FALSE;
56 :
57 : /******************************************************************************/
58 : /*! \page gdalwarp gdalwarp
59 :
60 : image reprojection and warping utility
61 :
62 : \section wsynopsis SYNOPSIS
63 :
64 : \htmlonly
65 : Usage:
66 : \endhtmlonly
67 :
68 : \verbatim
69 : gdalwarp [--help-general] [--formats]
70 : [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]
71 : [-order n] [-tps] [-rpc] [-geoloc] [-et err_threshold]
72 : [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]
73 : [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
74 : [-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha
75 : [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
76 : [-cutline datasource] [-cl layer] [-cwhere expression]
77 : [-csql statement] [-cblend dist_in_pixels]
78 : [-of format] [-co "NAME=VALUE"]*
79 : srcfile* dstfile
80 : \endverbatim
81 :
82 : \section wdescription DESCRIPTION
83 :
84 : <p>
85 : The gdalwarp utility is an image mosaicing, reprojection and warping
86 : utility. The program can reproject to any supported projection,
87 : and can also apply GCPs stored with the image if the image is "raw"
88 : with control information.
89 :
90 : <p>
91 : <dl>
92 : <dt> <b>-s_srs</b> <em>srs def</em>:</dt><dd> source spatial reference set.
93 : The coordinate systems that can be passed are anything supported by the
94 : OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes
95 : (ie. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file
96 : containing well known text.</dd>
97 : <dt> <b>-t_srs</b> <em>srs_def</em>:</dt><dd> target spatial reference set.
98 : The coordinate systems that can be passed are anything supported by the
99 : OGRSpatialReference.SetFromUserInput() call, which includes EPSG PCS and GCSes
100 : (ie. EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file
101 : containing well known text.</dd>
102 : <dt> <b>-to</b> <em>NAME=VALUE</em>:</dt><dd> set a transformer option suitable
103 : to pass to GDALCreateGenImgProjTransformer2(). </dd>
104 : <dt> <b>-order</b> <em>n</em>:</dt><dd> order of polynomial used for warping
105 : (1 to 3). The default is to select a polynomial order based on the number of
106 : GCPs.</dd>
107 : <dt> <b>-tps</b>:</dt><dd>Force use of thin plate spline transformer based on
108 : available GCPs.</dd>
109 : <dt> <b>-rpc</b>:</dt> <dd>Force use of RPCs.</dd>
110 : <dt> <b>-geoloc</b>:</dt><dd>Force use of Geolocation Arrays.</dd>
111 : <dt> <b>-et</b> <em>err_threshold</em>:</dt><dd> error threshold for
112 : transformation approximation (in pixel units - defaults to 0.125).</dd>
113 : <dt> <b>-te</b> <em>xmin ymin xmax ymax</em>:</dt><dd> set georeferenced
114 : extents of output file to be created (in target SRS).</dd>
115 : <dt> <b>-tr</b> <em>xres yres</em>:</dt><dd> set output file resolution (in
116 : target georeferenced units)</dd>
117 : <dt> <b>-ts</b> <em>width height</em>:</dt><dd> set output file size in
118 : pixels and lines. If width or height is set to 0, the other dimension will be
119 : guessed from the computed resolution. Note that -ts cannot be used with -tr</dd>
120 : <dt> <b>-wo</b> <em>"NAME=VALUE"</em>:</dt><dd> Set a warp options. The
121 : GDALWarpOptions::papszWarpOptions docs show all options. Multiple
122 : <b>-wo</b> options may be listed.</dd>
123 : <dt> <b>-ot</b> <em>type</em>:</dt><dd> For the output bands to be of the
124 : indicated data type.</dd>
125 : <dt> <b>-wt</b> <em>type</em>:</dt><dd> Working pixel data type. The data type
126 : of pixels in the source image and destination image buffers.</dd>
127 : <dt> <b>-r</b> <em>resampling_method</em>:</dt><dd> Resampling method to use. Available methods are:
128 : <dl>
129 : <dt><b>near</b></dt>: <dd>nearest neighbour resampling (default, fastest
130 : algorithm, worst interpolation quality).</dd>
131 : <dt><b>bilinear</b></dt>: <dd>bilinear resampling.</dd>
132 : <dt><b>cubic</b></dt>: <dd>cubic resampling.</dd>
133 : <dt><b>cubicspline</b></dt>: <dd>cubic spline resampling.</dd>
134 : <dt><b>lanczos</b></dt>: <dd>Lanczos windowed sinc resampling.</dd>
135 : </dl>
136 : <dt> <b>-srcnodata</b> <em>value [value...]</em>:</dt><dd> Set nodata masking
137 : values for input bands (different values can be supplied for each band). If
138 : more than one value is supplied all values should be quoted to keep them
139 : together as a single operating system argument. Masked values will not be
140 : used in interpolation. Use a value of <tt>None</tt> to ignore intrinsic nodata settings on the source dataset.</dd>
141 : <dt> <b>-dstnodata</b> <em>value [value...]</em>:</dt><dd> Set nodata values
142 : for output bands (different values can be supplied for each band). If more
143 : than one value is supplied all values should be quoted to keep them together
144 : as a single operating system argument. New files will be initialized to this
145 : value and if possible the nodata value will be recorded in the output
146 : file.</dd>
147 : <dt> <b>-dstalpha</b>:</dt><dd> Create an output alpha band to identify
148 : nodata (unset/transparent) pixels. </dd>
149 : <dt> <b>-wm</b> <em>memory_in_mb</em>:</dt><dd> Set the amount of memory (in
150 : megabytes) that the warp API is allowed to use for caching.</dd>
151 : <dt> <b>-multi</b>:</dt><dd> Use multithreaded warping implementation.
152 : Multiple threads will be used to process chunks of image and perform
153 : input/output operation simultaneously.</dd>
154 : <dt> <b>-q</b>:</dt><dd> Be quiet.</dd>
155 : <dt> <b>-of</b> <em>format</em>:</dt><dd> Select the output format. The default is GeoTIFF (GTiff). Use the short format name. </dd>
156 : <dt> <b>-co</b> <em>"NAME=VALUE"</em>:</dt><dd> passes a creation option to
157 : the output format driver. Multiple <b>-co</b> options may be listed. See
158 : format specific documentation for legal creation options for each format.
159 : </dd>
160 :
161 : <dt> <b>-cutline</b> <em>datasource</em>:</dt><dd>Enable use of a blend cutline from the name OGR support datasource.</dd>
162 : <dt> <b>-cl</b> <em>layername</em>:</dt><dd>Select the named layer from the
163 : cutline datasource.</dd>
164 : <dt> <b>-cwhere</b> <em>expression</em>:</dt><dd>Restrict desired cutline features based on attribute query.</dd>
165 : <dt> <b>-csql</b> <em>query</em>:</dt><dd>Select cutline features using an SQL query instead of from a layer with -cl.</dd>
166 : <dt> <b>-cblend</b> <em>distance</em>:</dt><dd>Set a blend distance to use to blend over cutlines (in pixels).</dd>
167 :
168 : <dt> <em>srcfile</em>:</dt><dd> The source file name(s). </dd>
169 : <dt> <em>dstfile</em>:</dt><dd> The destination file name. </dd>
170 : </dl>
171 :
172 : Mosaicing into an existing output file is supported if the output file
173 : already exists. The spatial extent of the existing file will not
174 : be modified to accomodate new data, so you may have to remove it in that case.
175 :
176 : Polygon cutlines may be used to restrict the the area of the destination file
177 : that may be updated, including blending. Cutline features must be in the
178 : georeferenced units of the destination file.
179 :
180 : <p>
181 : \section wexample EXAMPLE
182 :
183 : For instance, an eight bit spot scene stored in GeoTIFF with
184 : control points mapping the corners to lat/long could be warped to a UTM
185 : projection with a command like this:<p>
186 :
187 : \verbatim
188 : gdalwarp -t_srs '+proj=utm +zone=11 +datum=WGS84' raw_spot.tif utm11.tif
189 : \endverbatim
190 :
191 : For instance, the second channel of an ASTER image stored in HDF with
192 : control points mapping the corners to lat/long could be warped to a UTM
193 : projection with a command like this:<p>
194 :
195 : \verbatim
196 : gdalwarp HDF4_SDS:ASTER_L1B:"pg-PR1B0000-2002031402_100_001":2 pg-PR1B0000-2002031402_100_001_2.tif
197 : \endverbatim
198 :
199 : \if man
200 : \section wauthor AUTHORS
201 : Frank Warmerdam <warmerdam@pobox.com>, Silke Reimer <silke@intevation.de>
202 : \endif
203 : */
204 :
205 : /************************************************************************/
206 : /* Usage() */
207 : /************************************************************************/
208 :
209 0 : static void Usage()
210 :
211 : {
212 : printf(
213 : "Usage: gdalwarp [--help-general] [--formats]\n"
214 : " [-s_srs srs_def] [-t_srs srs_def] [-to \"NAME=VALUE\"]\n"
215 : " [-order n] [-tps] [-rpc] [-geoloc] [-et err_threshold]\n"
216 : " [-te xmin ymin xmax ymax] [-tr xres yres] [-ts width height]\n"
217 : " [-wo \"NAME=VALUE\"] [-ot Byte/Int16/...] [-wt Byte/Int16]\n"
218 : " [-srcnodata \"value [value...]\"] [-dstnodata \"value [value...]\"] -dstalpha\n"
219 : " [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]\n"
220 : " [-cutline datasource] [-cl layer] [-cwhere expression]\n"
221 : " [-csql statement] [-cblend dist_in_pixels]\n"
222 :
223 : " [-of format] [-co \"NAME=VALUE\"]*\n"
224 : " srcfile* dstfile\n"
225 : "\n"
226 : "Available resampling methods:\n"
227 0 : " near (default), bilinear, cubic, cubicspline, lanczos.\n" );
228 0 : exit( 1 );
229 : }
230 :
231 : /************************************************************************/
232 : /* SanitizeSRS */
233 : /************************************************************************/
234 :
235 10 : char *SanitizeSRS( const char *pszUserInput )
236 :
237 : {
238 : OGRSpatialReferenceH hSRS;
239 10 : char *pszResult = NULL;
240 :
241 10 : CPLErrorReset();
242 :
243 10 : hSRS = OSRNewSpatialReference( NULL );
244 10 : if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
245 10 : OSRExportToWkt( hSRS, &pszResult );
246 : else
247 : {
248 : CPLError( CE_Failure, CPLE_AppDefined,
249 : "Translating source or target SRS failed:\n%s",
250 0 : pszUserInput );
251 0 : exit( 1 );
252 : }
253 :
254 10 : OSRDestroySpatialReference( hSRS );
255 :
256 10 : return pszResult;
257 : }
258 :
259 : /************************************************************************/
260 : /* main() */
261 : /************************************************************************/
262 :
263 261 : int main( int argc, char ** argv )
264 :
265 : {
266 : GDALDatasetH hDstDS;
267 261 : const char *pszFormat = "GTiff";
268 261 : char **papszSrcFiles = NULL;
269 261 : char *pszDstFilename = NULL;
270 261 : int bCreateOutput = FALSE, i;
271 261 : void *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL;
272 261 : char **papszWarpOptions = NULL;
273 261 : double dfErrorThreshold = 0.125;
274 261 : double dfWarpMemoryLimit = 0.0;
275 261 : GDALTransformerFunc pfnTransformer = NULL;
276 261 : char **papszCreateOptions = NULL;
277 261 : GDALDataType eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown;
278 261 : GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
279 261 : const char *pszSrcNodata = NULL;
280 261 : const char *pszDstNodata = NULL;
281 261 : int bMulti = FALSE;
282 261 : char **papszTO = NULL;
283 261 : char *pszCutlineDSName = NULL;
284 261 : char *pszCLayer = NULL, *pszCWHERE = NULL, *pszCSQL = NULL;
285 261 : void *hCutline = NULL;
286 :
287 : /* Check that we are running against at least GDAL 1.6 */
288 : /* Note to developers : if we use newer API, please change the requirement */
289 261 : if (atoi(GDALVersionInfo("VERSION_NUM")) < 1600)
290 : {
291 : fprintf(stderr, "At least, GDAL >= 1.6.0 is required for this version of %s, "
292 0 : "which was compiled against GDAL %s\n", argv[0], GDAL_RELEASE_NAME);
293 0 : exit(1);
294 : }
295 :
296 : /* Must process GDAL_SKIP before GDALAllRegister(), but we can't call */
297 : /* GDALGeneralCmdLineProcessor before it needs the drivers to be registered */
298 : /* for the --format or --formats options */
299 1376 : for( i = 1; i < argc; i++ )
300 : {
301 1115 : if( EQUAL(argv[i],"--config") && i + 2 < argc && EQUAL(argv[i + 1], "GDAL_SKIP") )
302 : {
303 0 : CPLSetConfigOption( argv[i+1], argv[i+2] );
304 :
305 0 : i += 2;
306 : }
307 : }
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Register standard GDAL drivers, and process generic GDAL */
311 : /* command options. */
312 : /* -------------------------------------------------------------------- */
313 261 : GDALAllRegister();
314 261 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
315 261 : if( argc < 1 )
316 0 : exit( -argc );
317 :
318 : /* -------------------------------------------------------------------- */
319 : /* Parse arguments. */
320 : /* -------------------------------------------------------------------- */
321 1072 : for( i = 1; i < argc; i++ )
322 : {
323 812 : if( EQUAL(argv[i], "--utility_version") )
324 : {
325 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
326 1 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
327 1 : return 0;
328 : }
329 823 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
330 : {
331 12 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
332 12 : bCreateOutput = TRUE;
333 : }
334 805 : else if( EQUAL(argv[i],"-wo") && i < argc-1 )
335 : {
336 6 : papszWarpOptions = CSLAddString( papszWarpOptions, argv[++i] );
337 : }
338 793 : else if( EQUAL(argv[i],"-multi") )
339 : {
340 1 : bMulti = TRUE;
341 : }
342 792 : else if( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet"))
343 : {
344 0 : bQuiet = TRUE;
345 : }
346 792 : else if( EQUAL(argv[i],"-dstalpha") )
347 : {
348 1 : bEnableDstAlpha = TRUE;
349 : }
350 791 : else if( EQUAL(argv[i],"-srcalpha") )
351 : {
352 0 : bEnableSrcAlpha = TRUE;
353 : }
354 794 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
355 : {
356 3 : pszFormat = argv[++i];
357 3 : bCreateOutput = TRUE;
358 3 : if( EQUAL(pszFormat,"VRT") )
359 2 : bVRT = TRUE;
360 : }
361 798 : else if( EQUAL(argv[i],"-t_srs") && i < argc-1 )
362 : {
363 10 : char *pszSRS = SanitizeSRS(argv[++i]);
364 10 : papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszSRS );
365 10 : CPLFree( pszSRS );
366 : }
367 778 : else if( EQUAL(argv[i],"-s_srs") && i < argc-1 )
368 : {
369 0 : char *pszSRS = SanitizeSRS(argv[++i]);
370 0 : papszTO = CSLSetNameValue( papszTO, "SRC_SRS", pszSRS );
371 0 : CPLFree( pszSRS );
372 : }
373 778 : else if( EQUAL(argv[i],"-order") && i < argc-1 )
374 : {
375 0 : papszTO = CSLSetNameValue( papszTO, "MAX_GCP_ORDER", argv[++i] );
376 : }
377 778 : else if( EQUAL(argv[i],"-tps") )
378 : {
379 1 : papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
380 : }
381 777 : else if( EQUAL(argv[i],"-rpc") )
382 : {
383 0 : papszTO = CSLSetNameValue( papszTO, "METHOD", "RPC" );
384 : }
385 777 : else if( EQUAL(argv[i],"-geoloc") )
386 : {
387 0 : papszTO = CSLSetNameValue( papszTO, "METHOD", "GEOLOC_ARRAY" );
388 : }
389 777 : else if( EQUAL(argv[i],"-to") && i < argc-1 )
390 : {
391 0 : papszTO = CSLAddString( papszTO, argv[++i] );
392 : }
393 779 : else if( EQUAL(argv[i],"-et") && i < argc-1 )
394 : {
395 2 : dfErrorThreshold = CPLAtofM(argv[++i]);
396 : }
397 784 : else if( EQUAL(argv[i],"-wm") && i < argc-1 )
398 : {
399 9 : if( CPLAtofM(argv[i+1]) < 10000 )
400 1 : dfWarpMemoryLimit = CPLAtofM(argv[i+1]) * 1024 * 1024;
401 : else
402 8 : dfWarpMemoryLimit = CPLAtofM(argv[i+1]);
403 9 : i++;
404 : }
405 767 : else if( EQUAL(argv[i],"-srcnodata") && i < argc-1 )
406 : {
407 1 : pszSrcNodata = argv[++i];
408 : }
409 766 : else if( EQUAL(argv[i],"-dstnodata") && i < argc-1 )
410 : {
411 1 : pszDstNodata = argv[++i];
412 : }
413 766 : else if( EQUAL(argv[i],"-tr") && i < argc-2 )
414 : {
415 2 : dfXRes = CPLAtofM(argv[++i]);
416 2 : dfYRes = fabs(CPLAtofM(argv[++i]));
417 2 : if( dfXRes == 0 || dfYRes == 0 )
418 : {
419 0 : printf( "Wrong value for -tr parameters\n");
420 0 : Usage();
421 0 : exit( 2 );
422 : }
423 2 : bCreateOutput = TRUE;
424 : }
425 763 : else if( EQUAL(argv[i],"-ot") && i < argc-1 )
426 : {
427 : int iType;
428 :
429 12 : for( iType = 1; iType < GDT_TypeCount; iType++ )
430 : {
431 11 : if( GDALGetDataTypeName((GDALDataType)iType) != NULL
432 : && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
433 : argv[i+1]) )
434 : {
435 1 : eOutputType = (GDALDataType) iType;
436 : }
437 : }
438 :
439 1 : if( eOutputType == GDT_Unknown )
440 : {
441 0 : printf( "Unknown output pixel type: %s\n", argv[i+1] );
442 0 : Usage();
443 0 : exit( 2 );
444 : }
445 1 : i++;
446 1 : bCreateOutput = TRUE;
447 : }
448 761 : else if( EQUAL(argv[i],"-wt") && i < argc-1 )
449 : {
450 : int iType;
451 :
452 0 : for( iType = 1; iType < GDT_TypeCount; iType++ )
453 : {
454 0 : if( GDALGetDataTypeName((GDALDataType)iType) != NULL
455 : && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
456 : argv[i+1]) )
457 : {
458 0 : eWorkingType = (GDALDataType) iType;
459 : }
460 : }
461 :
462 0 : if( eWorkingType == GDT_Unknown )
463 : {
464 0 : printf( "Unknown output pixel type: %s\n", argv[i+1] );
465 0 : Usage();
466 0 : exit( 2 );
467 : }
468 0 : i++;
469 : }
470 769 : else if( EQUAL(argv[i],"-ts") && i < argc-2 )
471 : {
472 8 : nForcePixels = atoi(argv[++i]);
473 8 : nForceLines = atoi(argv[++i]);
474 8 : bCreateOutput = TRUE;
475 : }
476 754 : else if( EQUAL(argv[i],"-te") && i < argc-4 )
477 : {
478 1 : dfMinX = CPLAtofM(argv[++i]);
479 1 : dfMinY = CPLAtofM(argv[++i]);
480 1 : dfMaxX = CPLAtofM(argv[++i]);
481 1 : dfMaxY = CPLAtofM(argv[++i]);
482 1 : bCreateOutput = TRUE;
483 : }
484 752 : else if( EQUAL(argv[i],"-rn") )
485 1 : eResampleAlg = GRA_NearestNeighbour;
486 :
487 751 : else if( EQUAL(argv[i],"-rb") )
488 1 : eResampleAlg = GRA_Bilinear;
489 :
490 750 : else if( EQUAL(argv[i],"-rc") )
491 1 : eResampleAlg = GRA_Cubic;
492 :
493 749 : else if( EQUAL(argv[i],"-rcs") )
494 1 : eResampleAlg = GRA_CubicSpline;
495 :
496 970 : else if( EQUAL(argv[i],"-r") && i < argc - 1 )
497 : {
498 222 : if ( EQUAL(argv[++i], "near") )
499 44 : eResampleAlg = GRA_NearestNeighbour;
500 178 : else if ( EQUAL(argv[i], "bilinear") )
501 45 : eResampleAlg = GRA_Bilinear;
502 133 : else if ( EQUAL(argv[i], "cubic") )
503 44 : eResampleAlg = GRA_Cubic;
504 89 : else if ( EQUAL(argv[i], "cubicspline") )
505 44 : eResampleAlg = GRA_CubicSpline;
506 45 : else if ( EQUAL(argv[i], "lanczos") )
507 45 : eResampleAlg = GRA_Lanczos;
508 : else
509 : {
510 0 : printf( "Unknown resampling method: \"%s\".\n", argv[i] );
511 0 : Usage();
512 : }
513 : }
514 :
515 529 : else if( EQUAL(argv[i],"-cutline") && i < argc-1 )
516 : {
517 3 : pszCutlineDSName = argv[++i];
518 : }
519 523 : else if( EQUAL(argv[i],"-cwhere") && i < argc-1 )
520 : {
521 0 : pszCWHERE = argv[++i];
522 : }
523 526 : else if( EQUAL(argv[i],"-cl") && i < argc-1 )
524 : {
525 3 : pszCLayer = argv[++i];
526 : }
527 520 : else if( EQUAL(argv[i],"-csql") && i < argc-1 )
528 : {
529 0 : pszCSQL = argv[++i];
530 : }
531 520 : else if( EQUAL(argv[i],"-cblend") && i < argc-1 )
532 : {
533 : papszWarpOptions =
534 : CSLSetNameValue( papszWarpOptions,
535 0 : "CUTLINE_BLEND_DIST", argv[++i] );
536 : }
537 520 : else if( argv[i][0] == '-' )
538 0 : Usage();
539 :
540 : else
541 520 : papszSrcFiles = CSLAddString( papszSrcFiles, argv[i] );
542 : }
543 : /* -------------------------------------------------------------------- */
544 : /* Check that incompatible options are not used */
545 : /* -------------------------------------------------------------------- */
546 :
547 260 : if ((nForcePixels != 0 || nForceLines != 0) &&
548 : (dfXRes != 0 && dfYRes != 0))
549 : {
550 0 : printf( "-tr and -ts options cannot be used at the same time\n");
551 0 : Usage();
552 0 : exit( 2 );
553 : }
554 :
555 : /* -------------------------------------------------------------------- */
556 : /* The last filename in the file list is really our destination */
557 : /* file. */
558 : /* -------------------------------------------------------------------- */
559 260 : if( CSLCount(papszSrcFiles) > 1 )
560 : {
561 260 : pszDstFilename = papszSrcFiles[CSLCount(papszSrcFiles)-1];
562 260 : papszSrcFiles[CSLCount(papszSrcFiles)-1] = NULL;
563 : }
564 :
565 260 : if( pszDstFilename == NULL )
566 0 : Usage();
567 :
568 260 : if( bVRT && CSLCount(papszSrcFiles) > 1 )
569 : {
570 : fprintf(stderr, "Warning: gdalwarp -of VRT just takes into account "
571 : "the first source dataset.\nIf all source datasets "
572 : "are in the same projection, try making a mosaic of\n"
573 : "them with gdalbuildvrt, and use the resulting "
574 0 : "VRT file as the input of\ngdalwarp -of VRT.\n");
575 : }
576 :
577 : /* -------------------------------------------------------------------- */
578 : /* Does the output dataset already exist? */
579 : /* -------------------------------------------------------------------- */
580 260 : CPLPushErrorHandler( CPLQuietErrorHandler );
581 260 : hDstDS = GDALOpen( pszDstFilename, GA_Update );
582 260 : CPLPopErrorHandler();
583 :
584 260 : if( hDstDS != NULL && bCreateOutput )
585 : {
586 : fprintf( stderr,
587 : "Output dataset %s exists,\n"
588 : "but some commandline options were provided indicating a new dataset\n"
589 : "should be created. Please delete existing dataset and run again.\n",
590 0 : pszDstFilename );
591 0 : exit( 1 );
592 : }
593 :
594 : /* Avoid overwriting an existing destination file that cannot be opened in */
595 : /* update mode with a new GTiff file */
596 260 : if ( hDstDS == NULL )
597 : {
598 260 : CPLPushErrorHandler( CPLQuietErrorHandler );
599 260 : hDstDS = GDALOpen( pszDstFilename, GA_ReadOnly );
600 260 : CPLPopErrorHandler();
601 :
602 260 : if (hDstDS)
603 : {
604 : fprintf( stderr,
605 : "Output dataset %s exists, but cannot be opened in update mode\n",
606 0 : pszDstFilename );
607 0 : GDALClose(hDstDS);
608 0 : exit( 1 );
609 : }
610 : }
611 :
612 : /* -------------------------------------------------------------------- */
613 : /* If not, we need to create it. */
614 : /* -------------------------------------------------------------------- */
615 260 : int bInitDestSetForFirst = FALSE;
616 :
617 260 : if( hDstDS == NULL )
618 : {
619 : hDstDS = GDALWarpCreateOutput( papszSrcFiles, pszDstFilename,pszFormat,
620 : papszTO, &papszCreateOptions,
621 260 : eOutputType );
622 260 : bCreateOutput = TRUE;
623 :
624 260 : if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL
625 : && pszDstNodata == NULL )
626 : {
627 : papszWarpOptions = CSLSetNameValue(papszWarpOptions,
628 258 : "INIT_DEST", "0");
629 258 : bInitDestSetForFirst = TRUE;
630 : }
631 2 : else if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL )
632 : {
633 : papszWarpOptions = CSLSetNameValue(papszWarpOptions,
634 1 : "INIT_DEST", "NO_DATA" );
635 1 : bInitDestSetForFirst = TRUE;
636 : }
637 :
638 260 : CSLDestroy( papszCreateOptions );
639 260 : papszCreateOptions = NULL;
640 : }
641 :
642 260 : if( hDstDS == NULL )
643 0 : exit( 1 );
644 :
645 : /* -------------------------------------------------------------------- */
646 : /* If we have a cutline datasource read it and attach it in the */
647 : /* warp options. */
648 : /* -------------------------------------------------------------------- */
649 260 : if( pszCutlineDSName != NULL )
650 : {
651 : LoadCutline( pszCutlineDSName, pszCLayer, pszCWHERE, pszCSQL,
652 3 : &hCutline );
653 : }
654 :
655 : /* -------------------------------------------------------------------- */
656 : /* Loop over all source files, processing each in turn. */
657 : /* -------------------------------------------------------------------- */
658 : int iSrc;
659 :
660 1036 : for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
661 : {
662 : GDALDatasetH hSrcDS;
663 :
664 : /* -------------------------------------------------------------------- */
665 : /* Open this file. */
666 : /* -------------------------------------------------------------------- */
667 260 : hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
668 :
669 260 : if( hSrcDS == NULL )
670 0 : exit( 2 );
671 :
672 : /* -------------------------------------------------------------------- */
673 : /* Check that there's at least one raster band */
674 : /* -------------------------------------------------------------------- */
675 260 : if ( GDALGetRasterCount(hSrcDS) == 0 )
676 : {
677 0 : fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
678 0 : exit( 1 );
679 : }
680 :
681 260 : if( !bQuiet )
682 260 : printf( "Processing input file %s.\n", papszSrcFiles[iSrc] );
683 :
684 : /* -------------------------------------------------------------------- */
685 : /* Warns if the file has a color table and something more */
686 : /* complicated than nearest neighbour resampling is asked */
687 : /* -------------------------------------------------------------------- */
688 :
689 260 : if ( eResampleAlg != GRA_NearestNeighbour &&
690 : GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != NULL)
691 : {
692 0 : if( !bQuiet )
693 : fprintf( stderr, "Warning: Input file %s has a color table, which will likely lead to "
694 : "bad results when using a resampling method other than "
695 : "nearest neighbour. Converting the dataset prior to 24/32 bit "
696 0 : "is advised.\n", papszSrcFiles[iSrc] );
697 : }
698 :
699 : /* -------------------------------------------------------------------- */
700 : /* Do we have a source alpha band? */
701 : /* -------------------------------------------------------------------- */
702 260 : if( GDALGetRasterColorInterpretation(
703 : GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) )
704 : == GCI_AlphaBand
705 : && !bEnableSrcAlpha )
706 : {
707 1 : bEnableSrcAlpha = TRUE;
708 1 : if( !bQuiet )
709 : printf( "Using band %d of source image as alpha.\n",
710 1 : GDALGetRasterCount(hSrcDS) );
711 : }
712 :
713 : /* -------------------------------------------------------------------- */
714 : /* Create a transformation object from the source to */
715 : /* destination coordinate system. */
716 : /* -------------------------------------------------------------------- */
717 : hTransformArg = hGenImgProjArg =
718 260 : GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszTO );
719 :
720 260 : if( hTransformArg == NULL )
721 0 : exit( 1 );
722 :
723 260 : pfnTransformer = GDALGenImgProjTransform;
724 :
725 : /* -------------------------------------------------------------------- */
726 : /* Warp the transformer with a linear approximator unless the */
727 : /* acceptable error is zero. */
728 : /* -------------------------------------------------------------------- */
729 260 : if( dfErrorThreshold != 0.0 )
730 : {
731 : hTransformArg = hApproxArg =
732 : GDALCreateApproxTransformer( GDALGenImgProjTransform,
733 258 : hGenImgProjArg, dfErrorThreshold);
734 258 : pfnTransformer = GDALApproxTransform;
735 : }
736 :
737 : /* -------------------------------------------------------------------- */
738 : /* Clear temporary INIT_DEST settings after the first image. */
739 : /* -------------------------------------------------------------------- */
740 260 : if( bInitDestSetForFirst && iSrc == 1 )
741 : papszWarpOptions = CSLSetNameValue( papszWarpOptions,
742 0 : "INIT_DEST", NULL );
743 :
744 : /* -------------------------------------------------------------------- */
745 : /* Setup warp options. */
746 : /* -------------------------------------------------------------------- */
747 260 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
748 :
749 260 : psWO->papszWarpOptions = CSLDuplicate(papszWarpOptions);
750 260 : psWO->eWorkingDataType = eWorkingType;
751 260 : psWO->eResampleAlg = eResampleAlg;
752 :
753 260 : psWO->hSrcDS = hSrcDS;
754 260 : psWO->hDstDS = hDstDS;
755 :
756 260 : psWO->pfnTransformer = pfnTransformer;
757 260 : psWO->pTransformerArg = hTransformArg;
758 :
759 260 : if( !bQuiet )
760 260 : psWO->pfnProgress = GDALTermProgress;
761 :
762 260 : if( dfWarpMemoryLimit != 0.0 )
763 9 : psWO->dfWarpMemoryLimit = dfWarpMemoryLimit;
764 :
765 : /* -------------------------------------------------------------------- */
766 : /* Setup band mapping. */
767 : /* -------------------------------------------------------------------- */
768 260 : if( bEnableSrcAlpha )
769 1 : psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
770 : else
771 259 : psWO->nBandCount = GDALGetRasterCount(hSrcDS);
772 :
773 260 : psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
774 260 : psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
775 :
776 542 : for( i = 0; i < psWO->nBandCount; i++ )
777 : {
778 282 : psWO->panSrcBands[i] = i+1;
779 282 : psWO->panDstBands[i] = i+1;
780 : }
781 :
782 : /* -------------------------------------------------------------------- */
783 : /* Setup alpha bands used if any. */
784 : /* -------------------------------------------------------------------- */
785 260 : if( bEnableSrcAlpha )
786 1 : psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
787 :
788 260 : if( !bEnableDstAlpha
789 : && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1
790 : && GDALGetRasterColorInterpretation(
791 : GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS)))
792 : == GCI_AlphaBand )
793 : {
794 1 : if( !bQuiet )
795 : printf( "Using band %d of destination image as alpha.\n",
796 1 : GDALGetRasterCount(hDstDS) );
797 :
798 1 : bEnableDstAlpha = TRUE;
799 : }
800 :
801 260 : if( bEnableDstAlpha )
802 2 : psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
803 :
804 : /* -------------------------------------------------------------------- */
805 : /* Setup NODATA options. */
806 : /* -------------------------------------------------------------------- */
807 260 : if( pszSrcNodata != NULL && !EQUALN(pszSrcNodata,"n",1) )
808 : {
809 1 : char **papszTokens = CSLTokenizeString( pszSrcNodata );
810 1 : int nTokenCount = CSLCount(papszTokens);
811 :
812 : psWO->padfSrcNoDataReal = (double *)
813 1 : CPLMalloc(psWO->nBandCount*sizeof(double));
814 : psWO->padfSrcNoDataImag = (double *)
815 1 : CPLMalloc(psWO->nBandCount*sizeof(double));
816 :
817 4 : for( i = 0; i < psWO->nBandCount; i++ )
818 : {
819 3 : if( i < nTokenCount )
820 : {
821 : CPLStringToComplex( papszTokens[i],
822 : psWO->padfSrcNoDataReal + i,
823 1 : psWO->padfSrcNoDataImag + i );
824 : }
825 : else
826 : {
827 2 : psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i-1];
828 2 : psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i-1];
829 : }
830 : }
831 :
832 1 : CSLDestroy( papszTokens );
833 :
834 : psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
835 1 : "UNIFIED_SRC_NODATA", "YES" );
836 : }
837 :
838 : /* -------------------------------------------------------------------- */
839 : /* If -srcnodata was not specified, but the data has nodata */
840 : /* values, use them. */
841 : /* -------------------------------------------------------------------- */
842 260 : if( pszSrcNodata == NULL )
843 : {
844 259 : int bHaveNodata = FALSE;
845 259 : double dfReal = 0.0;
846 :
847 536 : for( i = 0; !bHaveNodata && i < psWO->nBandCount; i++ )
848 : {
849 277 : GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
850 277 : dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
851 : }
852 :
853 259 : if( bHaveNodata )
854 : {
855 1 : if( !bQuiet )
856 : {
857 1 : if (CPLIsNan(dfReal))
858 : printf( "Using internal nodata values (eg. nan) for image %s.\n",
859 0 : papszSrcFiles[iSrc] );
860 : else
861 : printf( "Using internal nodata values (eg. %g) for image %s.\n",
862 1 : dfReal, papszSrcFiles[iSrc] );
863 : }
864 : psWO->padfSrcNoDataReal = (double *)
865 1 : CPLMalloc(psWO->nBandCount*sizeof(double));
866 : psWO->padfSrcNoDataImag = (double *)
867 1 : CPLMalloc(psWO->nBandCount*sizeof(double));
868 :
869 4 : for( i = 0; i < psWO->nBandCount; i++ )
870 : {
871 3 : GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
872 :
873 3 : dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
874 :
875 3 : if( bHaveNodata )
876 : {
877 3 : psWO->padfSrcNoDataReal[i] = dfReal;
878 3 : psWO->padfSrcNoDataImag[i] = 0.0;
879 : }
880 : else
881 : {
882 0 : psWO->padfSrcNoDataReal[i] = -123456.789;
883 0 : psWO->padfSrcNoDataImag[i] = 0.0;
884 : }
885 : }
886 : }
887 : }
888 :
889 : /* -------------------------------------------------------------------- */
890 : /* If the output dataset was created, and we have a destination */
891 : /* nodata value, go through marking the bands with the information.*/
892 : /* -------------------------------------------------------------------- */
893 260 : if( pszDstNodata != NULL )
894 : {
895 1 : char **papszTokens = CSLTokenizeString( pszDstNodata );
896 1 : int nTokenCount = CSLCount(papszTokens);
897 :
898 : psWO->padfDstNoDataReal = (double *)
899 1 : CPLMalloc(psWO->nBandCount*sizeof(double));
900 : psWO->padfDstNoDataImag = (double *)
901 1 : CPLMalloc(psWO->nBandCount*sizeof(double));
902 :
903 2 : for( i = 0; i < psWO->nBandCount; i++ )
904 : {
905 1 : if( i < nTokenCount )
906 : {
907 : CPLStringToComplex( papszTokens[i],
908 : psWO->padfDstNoDataReal + i,
909 1 : psWO->padfDstNoDataImag + i );
910 : }
911 : else
912 : {
913 0 : psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i-1];
914 0 : psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i-1];
915 : }
916 :
917 1 : GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 );
918 1 : int bClamped = FALSE, bRounded = FALSE;
919 :
920 : #define CLAMP(val,type,minval,maxval) \
921 : do { if (val < minval) { bClamped = TRUE; val = minval; } \
922 : else if (val > maxval) { bClamped = TRUE; val = maxval; } \
923 : else if (val != (type)val) { bRounded = TRUE; val = (type)(val + 0.5); } } \
924 : while(0)
925 :
926 1 : switch(GDALGetRasterDataType(hBand))
927 : {
928 : case GDT_Byte:
929 1 : CLAMP(psWO->padfDstNoDataReal[i], GByte,
930 : 0.0, 255.0);
931 1 : break;
932 : case GDT_Int16:
933 0 : CLAMP(psWO->padfDstNoDataReal[i], GInt16,
934 : -32768.0, 32767.0);
935 0 : break;
936 : case GDT_UInt16:
937 0 : CLAMP(psWO->padfDstNoDataReal[i], GUInt16,
938 : 0.0, 65535.0);
939 0 : break;
940 : case GDT_Int32:
941 0 : CLAMP(psWO->padfDstNoDataReal[i], GInt32,
942 : -2147483648.0, 2147483647.0);
943 0 : break;
944 : case GDT_UInt32:
945 0 : CLAMP(psWO->padfDstNoDataReal[i], GUInt32,
946 : 0.0, 4294967295.0);
947 : break;
948 : default:
949 : break;
950 : }
951 :
952 1 : if (bClamped)
953 : {
954 : printf( "for band %d, destination nodata value has been clamped "
955 : "to %.0f, the original value being out of range.\n",
956 0 : i + 1, psWO->padfDstNoDataReal[i]);
957 : }
958 1 : else if(bRounded)
959 : {
960 : printf("for band %d, destination nodata value has been rounded "
961 : "to %.0f, %s being an integer datatype.\n",
962 : i + 1, psWO->padfDstNoDataReal[i],
963 0 : GDALGetDataTypeName(GDALGetRasterDataType(hBand)));
964 : }
965 :
966 1 : if( bCreateOutput )
967 : {
968 : GDALSetRasterNoDataValue(
969 : GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ),
970 1 : psWO->padfDstNoDataReal[i] );
971 : }
972 : }
973 :
974 1 : CSLDestroy( papszTokens );
975 : }
976 :
977 : /* -------------------------------------------------------------------- */
978 : /* If we have a cutline, transform it into the source */
979 : /* pixel/line coordinate system and insert into warp options. */
980 : /* -------------------------------------------------------------------- */
981 260 : if( hCutline != NULL )
982 : {
983 : TransformCutlineToSource( hSrcDS, hCutline,
984 : &(psWO->papszWarpOptions),
985 3 : papszTO );
986 : }
987 :
988 : /* -------------------------------------------------------------------- */
989 : /* If we are producing VRT output, then just initialize it with */
990 : /* the warp options and write out now rather than proceeding */
991 : /* with the operations. */
992 : /* -------------------------------------------------------------------- */
993 260 : if( bVRT )
994 : {
995 2 : if( GDALInitializeWarpedVRT( hDstDS, psWO ) != CE_None )
996 0 : exit( 1 );
997 :
998 2 : GDALClose( hDstDS );
999 2 : GDALClose( hSrcDS );
1000 :
1001 : /* The warped VRT will clean itself the transformer used */
1002 : /* So we have only to destroy the hGenImgProjArg if we */
1003 : /* have wrapped it inside the hApproxArg */
1004 2 : if (pfnTransformer == GDALApproxTransform)
1005 : {
1006 1 : if( hGenImgProjArg != NULL )
1007 1 : GDALDestroyGenImgProjTransformer( hGenImgProjArg );
1008 : }
1009 :
1010 2 : GDALDestroyWarpOptions( psWO );
1011 :
1012 2 : CPLFree( pszDstFilename );
1013 2 : CSLDestroy( argv );
1014 2 : CSLDestroy( papszSrcFiles );
1015 2 : CSLDestroy( papszWarpOptions );
1016 2 : CSLDestroy( papszTO );
1017 :
1018 2 : GDALDumpOpenDatasets( stderr );
1019 :
1020 2 : GDALDestroyDriverManager();
1021 :
1022 2 : return 0;
1023 : }
1024 :
1025 : /* -------------------------------------------------------------------- */
1026 : /* Initialize and execute the warp. */
1027 : /* -------------------------------------------------------------------- */
1028 258 : GDALWarpOperation oWO;
1029 :
1030 258 : if( oWO.Initialize( psWO ) == CE_None )
1031 : {
1032 258 : if( bMulti )
1033 : oWO.ChunkAndWarpMulti( 0, 0,
1034 : GDALGetRasterXSize( hDstDS ),
1035 1 : GDALGetRasterYSize( hDstDS ) );
1036 : else
1037 : oWO.ChunkAndWarpImage( 0, 0,
1038 : GDALGetRasterXSize( hDstDS ),
1039 257 : GDALGetRasterYSize( hDstDS ) );
1040 : }
1041 :
1042 : /* -------------------------------------------------------------------- */
1043 : /* Cleanup */
1044 : /* -------------------------------------------------------------------- */
1045 258 : if( hApproxArg != NULL )
1046 257 : GDALDestroyApproxTransformer( hApproxArg );
1047 :
1048 258 : if( hGenImgProjArg != NULL )
1049 258 : GDALDestroyGenImgProjTransformer( hGenImgProjArg );
1050 :
1051 258 : GDALDestroyWarpOptions( psWO );
1052 :
1053 258 : GDALClose( hSrcDS );
1054 258 : }
1055 :
1056 : /* -------------------------------------------------------------------- */
1057 : /* Final Cleanup. */
1058 : /* -------------------------------------------------------------------- */
1059 258 : GDALClose( hDstDS );
1060 :
1061 258 : CPLFree( pszDstFilename );
1062 258 : CSLDestroy( argv );
1063 258 : CSLDestroy( papszSrcFiles );
1064 258 : CSLDestroy( papszWarpOptions );
1065 258 : CSLDestroy( papszTO );
1066 :
1067 258 : GDALDumpOpenDatasets( stderr );
1068 :
1069 258 : GDALDestroyDriverManager();
1070 :
1071 : #ifdef OGR_ENABLED
1072 258 : if( hCutline != NULL )
1073 3 : OGR_G_DestroyGeometry( (OGRGeometryH) hCutline );
1074 258 : OGRCleanupAll();
1075 : #endif
1076 :
1077 258 : return 0;
1078 : }
1079 :
1080 : /************************************************************************/
1081 : /* GDALWarpCreateOutput() */
1082 : /* */
1083 : /* Create the output file based on various commandline options, */
1084 : /* and the input file. */
1085 : /************************************************************************/
1086 :
1087 : static GDALDatasetH
1088 : GDALWarpCreateOutput( char **papszSrcFiles, const char *pszFilename,
1089 : const char *pszFormat, char **papszTO,
1090 260 : char ***ppapszCreateOptions, GDALDataType eDT )
1091 :
1092 :
1093 : {
1094 : GDALDriverH hDriver;
1095 : GDALDatasetH hDstDS;
1096 : void *hTransformArg;
1097 260 : GDALColorTableH hCT = NULL;
1098 260 : double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0;
1099 260 : double dfWrkResX=0, dfWrkResY=0;
1100 260 : int nDstBandCount = 0;
1101 :
1102 : /* -------------------------------------------------------------------- */
1103 : /* Find the output driver. */
1104 : /* -------------------------------------------------------------------- */
1105 260 : hDriver = GDALGetDriverByName( pszFormat );
1106 260 : if( hDriver == NULL
1107 : || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
1108 : {
1109 : int iDr;
1110 :
1111 : printf( "Output driver `%s' not recognised or does not support\n",
1112 0 : pszFormat );
1113 : printf( "direct output file creation. The following format drivers are configured\n"
1114 0 : "and support direct output:\n" );
1115 :
1116 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
1117 : {
1118 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
1119 :
1120 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
1121 : {
1122 : printf( " %s: %s\n",
1123 : GDALGetDriverShortName( hDriver ),
1124 0 : GDALGetDriverLongName( hDriver ) );
1125 : }
1126 : }
1127 0 : printf( "\n" );
1128 0 : exit( 1 );
1129 : }
1130 :
1131 : /* -------------------------------------------------------------------- */
1132 : /* For virtual output files, we have to set a special subclass */
1133 : /* of dataset to create. */
1134 : /* -------------------------------------------------------------------- */
1135 260 : if( bVRT )
1136 : *ppapszCreateOptions =
1137 : CSLSetNameValue( *ppapszCreateOptions, "SUBCLASS",
1138 2 : "VRTWarpedDataset" );
1139 :
1140 : /* -------------------------------------------------------------------- */
1141 : /* Loop over all input files to collect extents. */
1142 : /* -------------------------------------------------------------------- */
1143 : int iSrc;
1144 260 : char *pszThisTargetSRS = (char*)CSLFetchNameValue( papszTO, "DST_SRS" );
1145 260 : if( pszThisTargetSRS != NULL )
1146 10 : pszThisTargetSRS = CPLStrdup( pszThisTargetSRS );
1147 :
1148 520 : for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
1149 : {
1150 : GDALDatasetH hSrcDS;
1151 260 : const char *pszThisSourceSRS = CSLFetchNameValue(papszTO,"SRC_SRS");
1152 :
1153 260 : hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
1154 260 : if( hSrcDS == NULL )
1155 0 : exit( 1 );
1156 :
1157 : /* -------------------------------------------------------------------- */
1158 : /* Check that there's at least one raster band */
1159 : /* -------------------------------------------------------------------- */
1160 260 : if ( GDALGetRasterCount(hSrcDS) == 0 )
1161 : {
1162 0 : fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
1163 0 : exit( 1 );
1164 : }
1165 :
1166 260 : if( eDT == GDT_Unknown )
1167 259 : eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
1168 :
1169 : /* -------------------------------------------------------------------- */
1170 : /* If we are processing the first file, and it has a color */
1171 : /* table, then we will copy it to the destination file. */
1172 : /* -------------------------------------------------------------------- */
1173 260 : if( iSrc == 0 )
1174 : {
1175 260 : nDstBandCount = GDALGetRasterCount(hSrcDS);
1176 260 : hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
1177 260 : if( hCT != NULL )
1178 : {
1179 0 : hCT = GDALCloneColorTable( hCT );
1180 0 : if( !bQuiet )
1181 : printf( "Copying color table from %s to new file.\n",
1182 0 : papszSrcFiles[iSrc] );
1183 : }
1184 : }
1185 :
1186 : /* -------------------------------------------------------------------- */
1187 : /* Get the sourcesrs from the dataset, if not set already. */
1188 : /* -------------------------------------------------------------------- */
1189 260 : if( pszThisSourceSRS == NULL )
1190 : {
1191 260 : const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" );
1192 :
1193 260 : if( GDALGetProjectionRef( hSrcDS ) != NULL
1194 : && strlen(GDALGetProjectionRef( hSrcDS )) > 0
1195 : && (pszMethod == NULL || EQUAL(pszMethod,"GEOTRANSFORM")) )
1196 245 : pszThisSourceSRS = GDALGetProjectionRef( hSrcDS );
1197 :
1198 15 : else if( GDALGetGCPProjection( hSrcDS ) != NULL
1199 : && strlen(GDALGetGCPProjection(hSrcDS)) > 0
1200 : && GDALGetGCPCount( hSrcDS ) > 1
1201 : && (pszMethod == NULL || EQUALN(pszMethod,"GCP_",4)) )
1202 15 : pszThisSourceSRS = GDALGetGCPProjection( hSrcDS );
1203 0 : else if( pszMethod != NULL && EQUAL(pszMethod,"RPC") )
1204 0 : pszThisSourceSRS = SRS_WKT_WGS84;
1205 : else
1206 0 : pszThisSourceSRS = "";
1207 : }
1208 :
1209 260 : if( pszThisTargetSRS == NULL )
1210 250 : pszThisTargetSRS = CPLStrdup( pszThisSourceSRS );
1211 :
1212 : /* -------------------------------------------------------------------- */
1213 : /* Create a transformation object from the source to */
1214 : /* destination coordinate system. */
1215 : /* -------------------------------------------------------------------- */
1216 : hTransformArg =
1217 260 : GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
1218 :
1219 260 : if( hTransformArg == NULL )
1220 : {
1221 0 : CPLFree( pszThisTargetSRS );
1222 0 : GDALClose( hSrcDS );
1223 0 : return NULL;
1224 : }
1225 :
1226 : /* -------------------------------------------------------------------- */
1227 : /* Get approximate output definition. */
1228 : /* -------------------------------------------------------------------- */
1229 : double adfThisGeoTransform[6];
1230 : double adfExtent[4];
1231 : int nThisPixels, nThisLines;
1232 :
1233 260 : if( GDALSuggestedWarpOutput2( hSrcDS,
1234 : GDALGenImgProjTransform, hTransformArg,
1235 : adfThisGeoTransform,
1236 : &nThisPixels, &nThisLines,
1237 : adfExtent, 0 ) != CE_None )
1238 : {
1239 0 : CPLFree( pszThisTargetSRS );
1240 0 : GDALClose( hSrcDS );
1241 0 : return NULL;
1242 : }
1243 :
1244 260 : if (CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", NULL ) == NULL)
1245 : {
1246 260 : double MinX = adfExtent[0];
1247 260 : double MaxX = adfExtent[2];
1248 260 : double MaxY = adfExtent[3];
1249 260 : double MinY = adfExtent[1];
1250 260 : int bSuccess = TRUE;
1251 :
1252 : /* Check that the the edges of the target image are in the validity area */
1253 : /* of the target projection */
1254 : #define N_STEPS 20
1255 : int i,j;
1256 5560 : for(i=0;i<=N_STEPS && bSuccess;i++)
1257 : {
1258 116440 : for(j=0;j<=N_STEPS && bSuccess;j++)
1259 : {
1260 111140 : double dfRatioI = i * 1.0 / N_STEPS;
1261 111140 : double dfRatioJ = j * 1.0 / N_STEPS;
1262 111140 : double expected_x = (1 - dfRatioI) * MinX + dfRatioI * MaxX;
1263 111140 : double expected_y = (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
1264 111140 : double x = expected_x;
1265 111140 : double y = expected_y;
1266 111140 : double z = 0;
1267 : /* Target SRS coordinates to source image pixel coordinates */
1268 111140 : if (!GDALGenImgProjTransform(hTransformArg, TRUE, 1, &x, &y, &z, &bSuccess) || !bSuccess)
1269 0 : bSuccess = FALSE;
1270 : /* Source image pixel coordinates to target SRS coordinates */
1271 111140 : if (!GDALGenImgProjTransform(hTransformArg, FALSE, 1, &x, &y, &z, &bSuccess) || !bSuccess)
1272 0 : bSuccess = FALSE;
1273 111140 : if (fabs(x - expected_x) > (MaxX - MinX) / nThisPixels ||
1274 : fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
1275 8 : bSuccess = FALSE;
1276 : }
1277 : }
1278 :
1279 : /* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces ogrct.cpp */
1280 : /* to check the consistency of each requested projection result with the */
1281 : /* invert projection */
1282 260 : if (!bSuccess)
1283 : {
1284 8 : CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" );
1285 8 : CPLDebug("WARP", "Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE");
1286 8 : GDALDestroyGenImgProjTransformer(hTransformArg);
1287 : hTransformArg =
1288 8 : GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
1289 :
1290 8 : if( GDALSuggestedWarpOutput2( hSrcDS,
1291 : GDALGenImgProjTransform, hTransformArg,
1292 : adfThisGeoTransform,
1293 : &nThisPixels, &nThisLines,
1294 : adfExtent, 0 ) != CE_None )
1295 : {
1296 0 : CPLFree( pszThisTargetSRS );
1297 0 : GDALClose( hSrcDS );
1298 0 : return NULL;
1299 : }
1300 : }
1301 : }
1302 :
1303 : /* -------------------------------------------------------------------- */
1304 : /* Expand the working bounds to include this region, ensure the */
1305 : /* working resolution is no more than this resolution. */
1306 : /* -------------------------------------------------------------------- */
1307 520 : if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 )
1308 : {
1309 260 : dfWrkMinX = adfExtent[0];
1310 260 : dfWrkMaxX = adfExtent[2];
1311 260 : dfWrkMaxY = adfExtent[3];
1312 260 : dfWrkMinY = adfExtent[1];
1313 260 : dfWrkResX = adfThisGeoTransform[1];
1314 260 : dfWrkResY = ABS(adfThisGeoTransform[5]);
1315 : }
1316 : else
1317 : {
1318 0 : dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]);
1319 0 : dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]);
1320 0 : dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]);
1321 0 : dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]);
1322 0 : dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]);
1323 0 : dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5]));
1324 : }
1325 :
1326 260 : GDALDestroyGenImgProjTransformer( hTransformArg );
1327 :
1328 260 : GDALClose( hSrcDS );
1329 : }
1330 :
1331 : /* -------------------------------------------------------------------- */
1332 : /* Did we have any usable sources? */
1333 : /* -------------------------------------------------------------------- */
1334 260 : if( nDstBandCount == 0 )
1335 : {
1336 : CPLError( CE_Failure, CPLE_AppDefined,
1337 0 : "No usable source images." );
1338 0 : CPLFree( pszThisTargetSRS );
1339 0 : return NULL;
1340 : }
1341 :
1342 : /* -------------------------------------------------------------------- */
1343 : /* Turn the suggested region into a geotransform and suggested */
1344 : /* number of pixels and lines. */
1345 : /* -------------------------------------------------------------------- */
1346 : double adfDstGeoTransform[6];
1347 : int nPixels, nLines;
1348 :
1349 260 : adfDstGeoTransform[0] = dfWrkMinX;
1350 260 : adfDstGeoTransform[1] = dfWrkResX;
1351 260 : adfDstGeoTransform[2] = 0.0;
1352 260 : adfDstGeoTransform[3] = dfWrkMaxY;
1353 260 : adfDstGeoTransform[4] = 0.0;
1354 260 : adfDstGeoTransform[5] = -1 * dfWrkResY;
1355 :
1356 260 : nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
1357 260 : nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
1358 :
1359 : /* -------------------------------------------------------------------- */
1360 : /* Did the user override some parameters? */
1361 : /* -------------------------------------------------------------------- */
1362 262 : if( dfXRes != 0.0 && dfYRes != 0.0 )
1363 : {
1364 2 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1365 : {
1366 2 : dfMinX = adfDstGeoTransform[0];
1367 2 : dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
1368 2 : dfMaxY = adfDstGeoTransform[3];
1369 2 : dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
1370 : }
1371 :
1372 2 : nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1373 2 : nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1374 2 : adfDstGeoTransform[0] = dfMinX;
1375 2 : adfDstGeoTransform[3] = dfMaxY;
1376 2 : adfDstGeoTransform[1] = dfXRes;
1377 2 : adfDstGeoTransform[5] = -dfYRes;
1378 : }
1379 :
1380 266 : else if( nForcePixels != 0 && nForceLines != 0 )
1381 : {
1382 8 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1383 : {
1384 8 : dfMinX = dfWrkMinX;
1385 8 : dfMaxX = dfWrkMaxX;
1386 8 : dfMaxY = dfWrkMaxY;
1387 8 : dfMinY = dfWrkMinY;
1388 : }
1389 :
1390 8 : dfXRes = (dfMaxX - dfMinX) / nForcePixels;
1391 8 : dfYRes = (dfMaxY - dfMinY) / nForceLines;
1392 :
1393 8 : adfDstGeoTransform[0] = dfMinX;
1394 8 : adfDstGeoTransform[3] = dfMaxY;
1395 8 : adfDstGeoTransform[1] = dfXRes;
1396 8 : adfDstGeoTransform[5] = -dfYRes;
1397 :
1398 8 : nPixels = nForcePixels;
1399 8 : nLines = nForceLines;
1400 : }
1401 :
1402 250 : else if( nForcePixels != 0 )
1403 : {
1404 0 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1405 : {
1406 0 : dfMinX = dfWrkMinX;
1407 0 : dfMaxX = dfWrkMaxX;
1408 0 : dfMaxY = dfWrkMaxY;
1409 0 : dfMinY = dfWrkMinY;
1410 : }
1411 :
1412 0 : dfXRes = (dfMaxX - dfMinX) / nForcePixels;
1413 0 : dfYRes = dfXRes;
1414 :
1415 0 : adfDstGeoTransform[0] = dfMinX;
1416 0 : adfDstGeoTransform[3] = dfMaxY;
1417 0 : adfDstGeoTransform[1] = dfXRes;
1418 0 : adfDstGeoTransform[5] = -dfYRes;
1419 :
1420 0 : nPixels = nForcePixels;
1421 0 : nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1422 : }
1423 :
1424 250 : else if( nForceLines != 0 )
1425 : {
1426 0 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1427 : {
1428 0 : dfMinX = dfWrkMinX;
1429 0 : dfMaxX = dfWrkMaxX;
1430 0 : dfMaxY = dfWrkMaxY;
1431 0 : dfMinY = dfWrkMinY;
1432 : }
1433 :
1434 0 : dfYRes = (dfMaxY - dfMinY) / nForceLines;
1435 0 : dfXRes = dfYRes;
1436 :
1437 0 : adfDstGeoTransform[0] = dfMinX;
1438 0 : adfDstGeoTransform[3] = dfMaxY;
1439 0 : adfDstGeoTransform[1] = dfXRes;
1440 0 : adfDstGeoTransform[5] = -dfYRes;
1441 :
1442 0 : nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1443 0 : nLines = nForceLines;
1444 : }
1445 :
1446 250 : else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
1447 : {
1448 1 : dfXRes = adfDstGeoTransform[1];
1449 1 : dfYRes = fabs(adfDstGeoTransform[5]);
1450 :
1451 1 : nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1452 1 : nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1453 :
1454 1 : dfXRes = (dfMaxX - dfMinX) / nPixels;
1455 1 : dfYRes = (dfMaxY - dfMinY) / nLines;
1456 :
1457 1 : adfDstGeoTransform[0] = dfMinX;
1458 1 : adfDstGeoTransform[3] = dfMaxY;
1459 1 : adfDstGeoTransform[1] = dfXRes;
1460 1 : adfDstGeoTransform[5] = -dfYRes;
1461 : }
1462 :
1463 : /* -------------------------------------------------------------------- */
1464 : /* Do we want to generate an alpha band in the output file? */
1465 : /* -------------------------------------------------------------------- */
1466 260 : if( bEnableSrcAlpha )
1467 0 : nDstBandCount--;
1468 :
1469 260 : if( bEnableDstAlpha )
1470 1 : nDstBandCount++;
1471 :
1472 : /* -------------------------------------------------------------------- */
1473 : /* Create the output file. */
1474 : /* -------------------------------------------------------------------- */
1475 260 : if( !bQuiet )
1476 260 : printf( "Creating output file that is %dP x %dL.\n", nPixels, nLines );
1477 :
1478 : hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
1479 260 : nDstBandCount, eDT, *ppapszCreateOptions );
1480 :
1481 260 : if( hDstDS == NULL )
1482 : {
1483 0 : CPLFree( pszThisTargetSRS );
1484 0 : return NULL;
1485 : }
1486 :
1487 : /* -------------------------------------------------------------------- */
1488 : /* Write out the projection definition. */
1489 : /* -------------------------------------------------------------------- */
1490 260 : GDALSetProjection( hDstDS, pszThisTargetSRS );
1491 260 : GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
1492 :
1493 : /* -------------------------------------------------------------------- */
1494 : /* Try to set color interpretation of output file alpha band. */
1495 : /* TODO: We should likely try to copy the other bands too. */
1496 : /* -------------------------------------------------------------------- */
1497 260 : if( bEnableDstAlpha )
1498 : {
1499 : GDALSetRasterColorInterpretation(
1500 : GDALGetRasterBand( hDstDS, nDstBandCount ),
1501 1 : GCI_AlphaBand );
1502 : }
1503 :
1504 : /* -------------------------------------------------------------------- */
1505 : /* Copy the color table, if required. */
1506 : /* -------------------------------------------------------------------- */
1507 260 : if( hCT != NULL )
1508 : {
1509 0 : GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
1510 0 : GDALDestroyColorTable( hCT );
1511 : }
1512 :
1513 260 : CPLFree( pszThisTargetSRS );
1514 260 : return hDstDS;
1515 : }
1516 :
1517 : /************************************************************************/
1518 : /* GeoTransform_Transformer() */
1519 : /* */
1520 : /* Convert points from georef coordinates to pixel/line based */
1521 : /* on a geotransform. */
1522 : /************************************************************************/
1523 :
1524 : class CutlineTransformer : public OGRCoordinateTransformation
1525 6 : {
1526 : public:
1527 :
1528 : void *hSrcImageTransformer;
1529 :
1530 0 : virtual OGRSpatialReference *GetSourceCS() { return NULL; }
1531 9 : virtual OGRSpatialReference *GetTargetCS() { return NULL; }
1532 :
1533 : virtual int Transform( int nCount,
1534 3 : double *x, double *y, double *z = NULL ) {
1535 : int nResult;
1536 :
1537 3 : int *pabSuccess = (int *) CPLCalloc(sizeof(int),nCount);
1538 3 : nResult = TransformEx( nCount, x, y, z, pabSuccess );
1539 3 : CPLFree( pabSuccess );
1540 :
1541 3 : return nResult;
1542 : }
1543 :
1544 : virtual int TransformEx( int nCount,
1545 : double *x, double *y, double *z = NULL,
1546 3 : int *pabSuccess = NULL ) {
1547 : return GDALGenImgProjTransform( hSrcImageTransformer, TRUE,
1548 3 : nCount, x, y, z, pabSuccess );
1549 : }
1550 : };
1551 :
1552 :
1553 : /************************************************************************/
1554 : /* LoadCutline() */
1555 : /* */
1556 : /* Load blend cutline from OGR datasource and attach in warp */
1557 : /* options, after potentially transforming to destination */
1558 : /* pixel/line coordinates. */
1559 : /************************************************************************/
1560 :
1561 : static void
1562 : LoadCutline( const char *pszCutlineDSName, const char *pszCLayer,
1563 : const char *pszCWHERE, const char *pszCSQL,
1564 3 : void **phCutlineRet )
1565 :
1566 : {
1567 : #ifndef OGR_ENABLED
1568 : CPLError( CE_Failure, CPLE_AppDefined,
1569 : "Request to load a cutline failed, this build does not support OGR features.\n" );
1570 : exit( 1 );
1571 : #else // def OGR_ENABLED
1572 3 : OGRRegisterAll();
1573 :
1574 : /* -------------------------------------------------------------------- */
1575 : /* Open source vector dataset. */
1576 : /* -------------------------------------------------------------------- */
1577 : OGRDataSourceH hSrcDS;
1578 :
1579 3 : hSrcDS = OGROpen( pszCutlineDSName, FALSE, NULL );
1580 3 : if( hSrcDS == NULL )
1581 0 : exit( 1 );
1582 :
1583 : /* -------------------------------------------------------------------- */
1584 : /* Get the source layer */
1585 : /* -------------------------------------------------------------------- */
1586 3 : OGRLayerH hLayer = NULL;
1587 :
1588 3 : if( pszCSQL != NULL )
1589 0 : hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszCSQL, NULL, NULL );
1590 3 : else if( pszCLayer != NULL )
1591 3 : hLayer = OGR_DS_GetLayerByName( hSrcDS, pszCLayer );
1592 : else
1593 0 : hLayer = OGR_DS_GetLayer( hSrcDS, 0 );
1594 :
1595 3 : if( hLayer == NULL )
1596 : {
1597 0 : fprintf( stderr, "Failed to identify source layer from datasource.\n" );
1598 0 : exit( 1 );
1599 : }
1600 :
1601 : /* -------------------------------------------------------------------- */
1602 : /* Apply WHERE clause if there is one. */
1603 : /* -------------------------------------------------------------------- */
1604 3 : if( pszCWHERE != NULL )
1605 0 : OGR_L_SetAttributeFilter( hLayer, pszCWHERE );
1606 :
1607 : /* -------------------------------------------------------------------- */
1608 : /* Collect the geometries from this layer, and build list of */
1609 : /* burn values. */
1610 : /* -------------------------------------------------------------------- */
1611 : OGRFeatureH hFeat;
1612 3 : OGRGeometryH hMultiPolygon = OGR_G_CreateGeometry( wkbMultiPolygon );
1613 :
1614 3 : OGR_L_ResetReading( hLayer );
1615 :
1616 9 : while( (hFeat = OGR_L_GetNextFeature( hLayer )) != NULL )
1617 : {
1618 3 : OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat);
1619 :
1620 3 : if( hGeom == NULL )
1621 : {
1622 0 : fprintf( stderr, "ERROR: Cutline feature without a geometry.\n" );
1623 0 : exit( 1 );
1624 : }
1625 :
1626 3 : OGRwkbGeometryType eType = wkbFlatten(OGR_G_GetGeometryType( hGeom ));
1627 :
1628 3 : if( eType == wkbPolygon )
1629 3 : OGR_G_AddGeometry( hMultiPolygon, hGeom );
1630 0 : else if( eType == wkbMultiPolygon )
1631 : {
1632 : int iGeom;
1633 :
1634 0 : for( iGeom = 0; iGeom < OGR_G_GetGeometryCount( hGeom ); iGeom++ )
1635 : {
1636 : OGR_G_AddGeometry( hMultiPolygon,
1637 0 : OGR_G_GetGeometryRef(hGeom,iGeom) );
1638 : }
1639 : }
1640 : else
1641 : {
1642 0 : fprintf( stderr, "ERROR: Cutline not of polygon type.\n" );
1643 0 : exit( 1 );
1644 : }
1645 :
1646 3 : OGR_F_Destroy( hFeat );
1647 : }
1648 :
1649 3 : if( OGR_G_GetGeometryCount( hMultiPolygon ) == 0 )
1650 : {
1651 0 : fprintf( stderr, "ERROR: Did not get any cutline features.\n" );
1652 0 : exit( 1 );
1653 : }
1654 :
1655 : /* -------------------------------------------------------------------- */
1656 : /* Ensure the coordinate system gets set on the geometry. */
1657 : /* -------------------------------------------------------------------- */
1658 : OGR_G_AssignSpatialReference(
1659 3 : hMultiPolygon, OGR_L_GetSpatialRef(hLayer) );
1660 :
1661 3 : *phCutlineRet = (void *) hMultiPolygon;
1662 :
1663 : /* -------------------------------------------------------------------- */
1664 : /* Cleanup */
1665 : /* -------------------------------------------------------------------- */
1666 3 : if( pszCSQL != NULL )
1667 0 : OGR_DS_ReleaseResultSet( hSrcDS, hLayer );
1668 :
1669 3 : OGR_DS_Destroy( hSrcDS );
1670 : #endif
1671 3 : }
1672 :
1673 : /************************************************************************/
1674 : /* TransformCutlineToSource() */
1675 : /************************************************************************/
1676 :
1677 : static void
1678 : TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
1679 3 : char ***ppapszWarpOptions, char **papszTO_In )
1680 :
1681 : {
1682 : #ifdef OGR_ENABLED
1683 3 : OGRGeometryH hMultiPolygon = OGR_G_Clone( (OGRGeometryH) hCutline );
1684 3 : char **papszTO = CSLDuplicate( papszTO_In );
1685 :
1686 : /* -------------------------------------------------------------------- */
1687 : /* Checkout that SRS are the same. */
1688 : /* -------------------------------------------------------------------- */
1689 3 : OGRSpatialReferenceH hRasterSRS = NULL;
1690 3 : const char *pszProjection = NULL;
1691 :
1692 3 : if( GDALGetProjectionRef( hSrcDS ) != NULL
1693 : && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
1694 3 : pszProjection = GDALGetProjectionRef( hSrcDS );
1695 0 : else if( GDALGetGCPProjection( hSrcDS ) != NULL )
1696 0 : pszProjection = GDALGetGCPProjection( hSrcDS );
1697 :
1698 3 : if( pszProjection != NULL )
1699 : {
1700 3 : hRasterSRS = OSRNewSpatialReference(NULL);
1701 3 : if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None )
1702 : {
1703 0 : OSRDestroySpatialReference(hRasterSRS);
1704 0 : hRasterSRS = NULL;
1705 : }
1706 : }
1707 :
1708 3 : OGRSpatialReferenceH hSrcSRS = OGR_G_GetSpatialReference( hMultiPolygon );
1709 3 : if( hRasterSRS != NULL && hSrcSRS != NULL )
1710 : {
1711 : /* ok, we will reproject */
1712 : }
1713 0 : else if( hRasterSRS != NULL && hSrcSRS == NULL )
1714 : {
1715 : fprintf(stderr,
1716 : "Warning : the source raster dataset has a SRS, but the input vector layer\n"
1717 0 : "not. Cutline results may be incorrect.\n");
1718 : }
1719 0 : else if( hRasterSRS == NULL && hSrcSRS != NULL )
1720 : {
1721 : fprintf(stderr,
1722 : "Warning : the input vector layer has a SRS, but the source raster dataset does not.\n"
1723 0 : "Cutline results may be incorrect.\n");
1724 : }
1725 :
1726 3 : if( hRasterSRS != NULL )
1727 3 : OSRDestroySpatialReference(hRasterSRS);
1728 :
1729 : /* -------------------------------------------------------------------- */
1730 : /* Extract the cutline SRS WKT. */
1731 : /* -------------------------------------------------------------------- */
1732 3 : if( hSrcSRS != NULL )
1733 : {
1734 3 : char *pszCutlineSRS_WKT = NULL;
1735 :
1736 3 : OSRExportToWkt( hSrcSRS, &pszCutlineSRS_WKT );
1737 3 : papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszCutlineSRS_WKT );
1738 3 : CPLFree( pszCutlineSRS_WKT );
1739 : }
1740 : else
1741 : {
1742 0 : int iDstSRS = CSLFindString( papszTO, "DST_SRS" );
1743 0 : if( iDstSRS >= 0 )
1744 0 : papszTO = CSLRemoveStrings( papszTO, iDstSRS, 1, NULL );
1745 : }
1746 :
1747 : /* -------------------------------------------------------------------- */
1748 : /* Transform the geometry to pixel/line coordinates. */
1749 : /* -------------------------------------------------------------------- */
1750 3 : CutlineTransformer oTransformer;
1751 :
1752 : oTransformer.hSrcImageTransformer =
1753 3 : GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
1754 :
1755 3 : CSLDestroy( papszTO );
1756 :
1757 3 : if( oTransformer.hSrcImageTransformer == NULL )
1758 0 : exit( 1 );
1759 :
1760 : OGR_G_Transform( hMultiPolygon,
1761 3 : (OGRCoordinateTransformationH) &oTransformer );
1762 :
1763 3 : GDALDestroyGenImgProjTransformer( oTransformer.hSrcImageTransformer );
1764 :
1765 : /* -------------------------------------------------------------------- */
1766 : /* Convert aggregate geometry into WKT. */
1767 : /* -------------------------------------------------------------------- */
1768 3 : char *pszWKT = NULL;
1769 :
1770 3 : OGR_G_ExportToWkt( hMultiPolygon, &pszWKT );
1771 3 : OGR_G_DestroyGeometry( hMultiPolygon );
1772 :
1773 : *ppapszWarpOptions = CSLSetNameValue( *ppapszWarpOptions,
1774 3 : "CUTLINE", pszWKT );
1775 3 : CPLFree( pszWKT );
1776 : #endif
1777 3 : }
|
