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