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