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