1 : /******************************************************************************
2 : * $Id: gdalwarp.cpp 24213 2012-04-08 20:11:22Z etourigny $
3 : *
4 : * Project: High Performance Image Reprojector
5 : * Purpose: Test program for high performance warper API.
6 : * Author: Frank Warmerdam <warmerdam@pobox.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2002, i3 - information integration and imaging
10 : * Fort Collin, CO
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "gdalwarper.h"
32 : #include "cpl_string.h"
33 : #include "ogr_spatialref.h"
34 : #include "ogr_api.h"
35 : #include "commonutils.h"
36 : #include <vector>
37 :
38 : CPL_CVSID("$Id: gdalwarp.cpp 24213 2012-04-08 20:11:22Z etourigny $");
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 >= 2.0.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 >= 2.0.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 2 : static int GDALExit( int nCode )
244 : {
245 2 : const char *pszDebug = CPLGetConfigOption("CPL_DEBUG",NULL);
246 2 : if( pszDebug && (EQUAL(pszDebug,"ON") || EQUAL(pszDebug,"") ) )
247 : {
248 0 : GDALDumpOpenDatasets( stderr );
249 0 : CPLDumpSharedList( NULL );
250 : }
251 :
252 2 : GDALDestroyDriverManager();
253 :
254 : #ifdef OGR_ENABLED
255 2 : OGRCleanupAll();
256 : #endif
257 :
258 2 : exit( nCode );
259 : }
260 :
261 : /************************************************************************/
262 : /* Usage() */
263 : /************************************************************************/
264 :
265 2 : 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 2 : " near (default), bilinear, cubic, cubicspline, lanczos.\n" );
285 2 : GDALExit( 1 );
286 0 : }
287 :
288 : /************************************************************************/
289 : /* SanitizeSRS */
290 : /************************************************************************/
291 :
292 26 : char *SanitizeSRS( const char *pszUserInput )
293 :
294 : {
295 : OGRSpatialReferenceH hSRS;
296 26 : char *pszResult = NULL;
297 :
298 26 : CPLErrorReset();
299 :
300 26 : hSRS = OSRNewSpatialReference( NULL );
301 26 : if( OSRSetFromUserInput( hSRS, pszUserInput ) == OGRERR_NONE )
302 26 : 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 26 : OSRDestroySpatialReference( hSRS );
312 :
313 26 : return pszResult;
314 : }
315 :
316 : /************************************************************************/
317 : /* main() */
318 : /************************************************************************/
319 :
320 542 : int main( int argc, char ** argv )
321 :
322 : {
323 : GDALDatasetH hDstDS;
324 542 : const char *pszFormat = "GTiff";
325 542 : int bFormatExplicitelySet = FALSE;
326 542 : char **papszSrcFiles = NULL;
327 542 : char *pszDstFilename = NULL;
328 542 : int bCreateOutput = FALSE, i;
329 542 : void *hTransformArg, *hGenImgProjArg=NULL, *hApproxArg=NULL;
330 542 : char **papszWarpOptions = NULL;
331 542 : double dfErrorThreshold = 0.125;
332 542 : double dfWarpMemoryLimit = 0.0;
333 542 : GDALTransformerFunc pfnTransformer = NULL;
334 542 : char **papszCreateOptions = NULL;
335 542 : GDALDataType eOutputType = GDT_Unknown, eWorkingType = GDT_Unknown;
336 542 : GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
337 542 : const char *pszSrcNodata = NULL;
338 542 : const char *pszDstNodata = NULL;
339 542 : int bMulti = FALSE;
340 542 : char **papszTO = NULL;
341 542 : char *pszCutlineDSName = NULL;
342 542 : char *pszCLayer = NULL, *pszCWHERE = NULL, *pszCSQL = NULL;
343 542 : void *hCutline = NULL;
344 542 : int bHasGotErr = FALSE;
345 542 : int bCropToCutline = FALSE;
346 542 : int bOverwrite = FALSE;
347 542 : int bCopyMetadata = TRUE;
348 542 : int bCopyBandInfo = TRUE;
349 542 : 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 542 : 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 2856 : for( i = 1; i < argc; i++ )
364 : {
365 2314 : 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 542 : GDALAllRegister();
378 542 : argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );
379 542 : if( argc < 1 )
380 0 : GDALExit( -argc );
381 :
382 : /* -------------------------------------------------------------------- */
383 : /* Parse arguments. */
384 : /* -------------------------------------------------------------------- */
385 2234 : for( i = 1; i < argc; i++ )
386 : {
387 1694 : if( EQUAL(argv[i],"-tps") || EQUAL(argv[i],"-rpc") || EQUAL(argv[i],"-geoloc") )
388 : {
389 4 : const char* pszMethod = CSLFetchNameValue(papszTO, "METHOD");
390 4 : if (pszMethod)
391 : fprintf(stderr, "Warning: only one METHOD can be used. Method %s is already defined.\n",
392 0 : pszMethod);
393 4 : const char* pszMAX_GCP_ORDER = CSLFetchNameValue(papszTO, "MAX_GCP_ORDER");
394 4 : 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 1694 : if( EQUAL(argv[i], "--utility_version") )
400 : {
401 : printf("%s was compiled against GDAL %s and is running against GDAL %s\n",
402 2 : argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
403 2 : return 0;
404 : }
405 1716 : else if( EQUAL(argv[i],"-co") && i < argc-1 )
406 : {
407 24 : papszCreateOptions = CSLAddString( papszCreateOptions, argv[++i] );
408 24 : bCreateOutput = TRUE;
409 : }
410 1680 : else if( EQUAL(argv[i],"-wo") && i < argc-1 )
411 : {
412 12 : papszWarpOptions = CSLAddString( papszWarpOptions, argv[++i] );
413 : }
414 1656 : else if( EQUAL(argv[i],"-multi") )
415 : {
416 2 : bMulti = TRUE;
417 : }
418 1656 : else if( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet"))
419 : {
420 2 : bQuiet = TRUE;
421 : }
422 1652 : else if( EQUAL(argv[i],"-dstalpha") )
423 : {
424 4 : bEnableDstAlpha = TRUE;
425 : }
426 1648 : else if( EQUAL(argv[i],"-srcalpha") )
427 : {
428 0 : bEnableSrcAlpha = TRUE;
429 : }
430 1658 : else if( EQUAL(argv[i],"-of") && i < argc-1 )
431 : {
432 10 : pszFormat = argv[++i];
433 10 : bFormatExplicitelySet = TRUE;
434 10 : bCreateOutput = TRUE;
435 10 : if( EQUAL(pszFormat,"VRT") )
436 8 : bVRT = TRUE;
437 : }
438 1664 : else if( EQUAL(argv[i],"-t_srs") && i < argc-1 )
439 : {
440 26 : char *pszSRS = SanitizeSRS(argv[++i]);
441 26 : papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszSRS );
442 26 : CPLFree( pszSRS );
443 : }
444 1612 : 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 1612 : 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 1612 : 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 1612 : else if( EQUAL(argv[i],"-tps") )
476 : {
477 4 : papszTO = CSLSetNameValue( papszTO, "METHOD", "GCP_TPS" );
478 : }
479 1608 : else if( EQUAL(argv[i],"-rpc") )
480 : {
481 0 : papszTO = CSLSetNameValue( papszTO, "METHOD", "RPC" );
482 : }
483 1608 : else if( EQUAL(argv[i],"-geoloc") )
484 : {
485 0 : papszTO = CSLSetNameValue( papszTO, "METHOD", "GEOLOC_ARRAY" );
486 : }
487 1608 : else if( EQUAL(argv[i],"-to") && i < argc-1 )
488 : {
489 0 : papszTO = CSLAddString( papszTO, argv[++i] );
490 : }
491 1612 : else if( EQUAL(argv[i],"-et") && i < argc-1 )
492 : {
493 4 : dfErrorThreshold = CPLAtofM(argv[++i]);
494 : }
495 1622 : else if( EQUAL(argv[i],"-wm") && i < argc-1 )
496 : {
497 18 : if( CPLAtofM(argv[i+1]) < 10000 )
498 2 : dfWarpMemoryLimit = CPLAtofM(argv[i+1]) * 1024 * 1024;
499 : else
500 16 : dfWarpMemoryLimit = CPLAtofM(argv[i+1]);
501 18 : i++;
502 : }
503 1588 : else if( EQUAL(argv[i],"-srcnodata") && i < argc-1 )
504 : {
505 2 : pszSrcNodata = argv[++i];
506 : }
507 1586 : else if( EQUAL(argv[i],"-dstnodata") && i < argc-1 )
508 : {
509 2 : pszDstNodata = argv[++i];
510 : }
511 1588 : else if( EQUAL(argv[i],"-tr") && i < argc-2 )
512 : {
513 6 : dfXRes = CPLAtofM(argv[++i]);
514 6 : dfYRes = fabs(CPLAtofM(argv[++i]));
515 6 : if( dfXRes == 0 || dfYRes == 0 )
516 : {
517 0 : printf( "Wrong value for -tr parameters\n");
518 0 : Usage();
519 : }
520 6 : bCreateOutput = TRUE;
521 : }
522 1576 : else if( EQUAL(argv[i],"-tap") )
523 : {
524 4 : bTargetAlignedPixels = TRUE;
525 : }
526 1574 : else if( EQUAL(argv[i],"-ot") && i < argc-1 )
527 : {
528 : int iType;
529 :
530 24 : for( iType = 1; iType < GDT_TypeCount; iType++ )
531 : {
532 44 : if( GDALGetDataTypeName((GDALDataType)iType) != NULL
533 22 : && EQUAL(GDALGetDataTypeName((GDALDataType)iType),
534 : argv[i+1]) )
535 : {
536 2 : eOutputType = (GDALDataType) iType;
537 : }
538 : }
539 :
540 2 : if( eOutputType == GDT_Unknown )
541 : {
542 0 : printf( "Unknown output pixel type: %s\n", argv[i+1] );
543 0 : Usage();
544 : }
545 2 : i++;
546 2 : bCreateOutput = TRUE;
547 : }
548 1570 : 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 1586 : else if( EQUAL(argv[i],"-ts") && i < argc-2 )
570 : {
571 16 : nForcePixels = atoi(argv[++i]);
572 16 : nForceLines = atoi(argv[++i]);
573 16 : bCreateOutput = TRUE;
574 : }
575 1556 : else if( EQUAL(argv[i],"-te") && i < argc-4 )
576 : {
577 2 : dfMinX = CPLAtofM(argv[++i]);
578 2 : dfMinY = CPLAtofM(argv[++i]);
579 2 : dfMaxX = CPLAtofM(argv[++i]);
580 2 : dfMaxY = CPLAtofM(argv[++i]);
581 2 : bCreateOutput = TRUE;
582 : }
583 1552 : else if( EQUAL(argv[i],"-rn") )
584 2 : eResampleAlg = GRA_NearestNeighbour;
585 :
586 1550 : else if( EQUAL(argv[i],"-rb") )
587 4 : eResampleAlg = GRA_Bilinear;
588 :
589 1546 : else if( EQUAL(argv[i],"-rc") )
590 2 : eResampleAlg = GRA_Cubic;
591 :
592 1544 : else if( EQUAL(argv[i],"-rcs") )
593 2 : eResampleAlg = GRA_CubicSpline;
594 :
595 1986 : else if( EQUAL(argv[i],"-r") && i < argc - 1 )
596 : {
597 444 : if ( EQUAL(argv[++i], "near") )
598 88 : eResampleAlg = GRA_NearestNeighbour;
599 356 : else if ( EQUAL(argv[i], "bilinear") )
600 90 : eResampleAlg = GRA_Bilinear;
601 266 : else if ( EQUAL(argv[i], "cubic") )
602 88 : eResampleAlg = GRA_Cubic;
603 178 : else if ( EQUAL(argv[i], "cubicspline") )
604 88 : eResampleAlg = GRA_CubicSpline;
605 90 : else if ( EQUAL(argv[i], "lanczos") )
606 90 : eResampleAlg = GRA_Lanczos;
607 : else
608 : {
609 0 : printf( "Unknown resampling method: \"%s\".\n", argv[i] );
610 0 : Usage();
611 : }
612 : }
613 :
614 1104 : else if( EQUAL(argv[i],"-cutline") && i < argc-1 )
615 : {
616 6 : pszCutlineDSName = argv[++i];
617 : }
618 1092 : else if( EQUAL(argv[i],"-cwhere") && i < argc-1 )
619 : {
620 0 : pszCWHERE = argv[++i];
621 : }
622 1098 : else if( EQUAL(argv[i],"-cl") && i < argc-1 )
623 : {
624 6 : pszCLayer = argv[++i];
625 : }
626 1086 : else if( EQUAL(argv[i],"-csql") && i < argc-1 )
627 : {
628 0 : pszCSQL = argv[++i];
629 : }
630 1086 : else if( EQUAL(argv[i],"-cblend") && i < argc-1 )
631 : {
632 : papszWarpOptions =
633 : CSLSetNameValue( papszWarpOptions,
634 0 : "CUTLINE_BLEND_DIST", argv[++i] );
635 : }
636 1086 : else if( EQUAL(argv[i],"-crop_to_cutline") )
637 : {
638 0 : bCropToCutline = TRUE;
639 0 : bCreateOutput = TRUE;
640 : }
641 1086 : else if( EQUAL(argv[i],"-overwrite") )
642 4 : bOverwrite = TRUE;
643 1082 : else if( EQUAL(argv[i],"-nomd") )
644 : {
645 0 : bCopyMetadata = FALSE;
646 0 : bCopyBandInfo = FALSE;
647 : }
648 1082 : else if( EQUAL(argv[i],"-cvmd") && i < argc-1 )
649 0 : pszMDConflictValue = argv[++i];
650 :
651 1082 : else if( argv[i][0] == '-' )
652 0 : Usage();
653 :
654 : else
655 1082 : papszSrcFiles = CSLAddString( papszSrcFiles, argv[i] );
656 : }
657 : /* -------------------------------------------------------------------- */
658 : /* Check that incompatible options are not used */
659 : /* -------------------------------------------------------------------- */
660 :
661 540 : 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 540 : if (bTargetAlignedPixels && dfXRes == 0 && dfYRes == 0)
669 : {
670 2 : printf( "-tap option cannot be used without using -tr\n");
671 2 : Usage();
672 : }
673 :
674 : /* -------------------------------------------------------------------- */
675 : /* The last filename in the file list is really our destination */
676 : /* file. */
677 : /* -------------------------------------------------------------------- */
678 538 : if( CSLCount(papszSrcFiles) > 1 )
679 : {
680 538 : pszDstFilename = papszSrcFiles[CSLCount(papszSrcFiles)-1];
681 538 : papszSrcFiles[CSLCount(papszSrcFiles)-1] = NULL;
682 : }
683 :
684 538 : if( pszDstFilename == NULL )
685 0 : Usage();
686 :
687 538 : 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 1074 : if ( CSLCount(papszSrcFiles) == 1 &&
704 536 : 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 538 : CPLPushErrorHandler( CPLQuietErrorHandler );
711 538 : hDstDS = GDALOpen( pszDstFilename, GA_Update );
712 538 : CPLPopErrorHandler();
713 :
714 538 : if( hDstDS != NULL && bOverwrite )
715 : {
716 2 : GDALClose(hDstDS);
717 2 : hDstDS = NULL;
718 : }
719 :
720 538 : 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 538 : if ( hDstDS == NULL && !bOverwrite )
733 : {
734 532 : CPLPushErrorHandler( CPLQuietErrorHandler );
735 532 : hDstDS = GDALOpen( pszDstFilename, GA_ReadOnly );
736 532 : CPLPopErrorHandler();
737 :
738 532 : 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 538 : if( pszCutlineDSName != NULL )
753 : {
754 : LoadCutline( pszCutlineDSName, pszCLayer, pszCWHERE, pszCSQL,
755 6 : &hCutline );
756 : }
757 :
758 : #ifdef OGR_ENABLED
759 538 : 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 538 : int bInitDestSetForFirst = FALSE;
851 :
852 538 : void* hUniqueTransformArg = NULL;
853 538 : GDALDatasetH hUniqueSrcDS = NULL;
854 :
855 538 : if( hDstDS == NULL )
856 : {
857 536 : if (!bQuiet && !bFormatExplicitelySet)
858 526 : CheckExtensionConsistency(pszDstFilename, pszFormat);
859 :
860 : hDstDS = GDALWarpCreateOutput( papszSrcFiles, pszDstFilename,pszFormat,
861 : papszTO, &papszCreateOptions,
862 : eOutputType, &hUniqueTransformArg,
863 536 : &hUniqueSrcDS);
864 536 : bCreateOutput = TRUE;
865 :
866 536 : if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL
867 : && pszDstNodata == NULL )
868 : {
869 : papszWarpOptions = CSLSetNameValue(papszWarpOptions,
870 532 : "INIT_DEST", "0");
871 532 : bInitDestSetForFirst = TRUE;
872 : }
873 4 : else if( CSLFetchNameValue( papszWarpOptions, "INIT_DEST" ) == NULL )
874 : {
875 : papszWarpOptions = CSLSetNameValue(papszWarpOptions,
876 2 : "INIT_DEST", "NO_DATA" );
877 2 : bInitDestSetForFirst = TRUE;
878 : }
879 :
880 536 : CSLDestroy( papszCreateOptions );
881 536 : papszCreateOptions = NULL;
882 : }
883 :
884 538 : 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 2140 : for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
893 : {
894 : GDALDatasetH hSrcDS;
895 :
896 : /* -------------------------------------------------------------------- */
897 : /* Open this file. */
898 : /* -------------------------------------------------------------------- */
899 540 : if (hUniqueSrcDS)
900 534 : hSrcDS = hUniqueSrcDS;
901 : else
902 6 : hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
903 :
904 540 : if( hSrcDS == NULL )
905 0 : GDALExit( 2 );
906 :
907 : /* -------------------------------------------------------------------- */
908 : /* Check that there's at least one raster band */
909 : /* -------------------------------------------------------------------- */
910 540 : 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 540 : if( !bQuiet )
917 538 : 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 540 : if ( bCopyMetadata )
926 : {
927 540 : char **papszMetadata = NULL;
928 540 : const char *pszSrcInfo = NULL;
929 540 : const char *pszDstInfo = NULL;
930 540 : GDALRasterBandH hSrcBand = NULL;
931 540 : GDALRasterBandH hDstBand = NULL;
932 :
933 : /* copy metadata from first dataset */
934 540 : if ( iSrc == 0 )
935 : {
936 538 : CPLDebug("WARP", "Copying metadata from first source to destination dataset");
937 : /* copy dataset-level metadata */
938 538 : papszMetadata = GDALGetMetadata( hSrcDS, NULL );
939 538 : if ( CSLCount(papszMetadata) > 0 ) {
940 532 : 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 538 : if ( GDALGetRasterCount( hSrcDS ) == GDALGetRasterCount( hDstDS ) )
946 : {
947 1124 : for ( int iBand = 0; iBand < GDALGetRasterCount( hSrcDS ); iBand++ )
948 : {
949 590 : hSrcBand = GDALGetRasterBand( hSrcDS, iBand + 1 );
950 590 : hDstBand = GDALGetRasterBand( hDstDS, iBand + 1 );
951 : /* copy metadata */
952 590 : papszMetadata = GDALGetMetadata( hSrcBand, NULL);
953 590 : if ( CSLCount(papszMetadata) > 0 )
954 16 : GDALSetMetadata( hDstBand, papszMetadata, NULL );
955 : /* copy other info (Description, Unit Type) - what else? */
956 590 : if ( bCopyBandInfo ) {
957 590 : pszSrcInfo = GDALGetDescription( hSrcBand );
958 590 : if( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 )
959 0 : GDALSetDescription( hDstBand, pszSrcInfo );
960 590 : pszSrcInfo = GDALGetRasterUnitType( hSrcBand );
961 590 : 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 2 : "Removing conflicting metadata from destination dataset (source #%d)", iSrc );
972 : /* remove conflicting dataset-level metadata */
973 2 : RemoveConflictingMetadata( hDstDS, GDALGetMetadata( hSrcDS, NULL ), pszMDConflictValue );
974 :
975 : /* remove conflicting copy band-level metadata and other info */
976 2 : if ( GDALGetRasterCount( hSrcDS ) == GDALGetRasterCount( hDstDS ) )
977 : {
978 4 : for ( int iBand = 0; iBand < GDALGetRasterCount( hSrcDS ); iBand++ )
979 : {
980 2 : hSrcBand = GDALGetRasterBand( hSrcDS, iBand + 1 );
981 2 : hDstBand = GDALGetRasterBand( hDstDS, iBand + 1 );
982 : /* remove conflicting metadata */
983 2 : RemoveConflictingMetadata( hDstBand, GDALGetMetadata( hSrcBand, NULL ), pszMDConflictValue );
984 : /* remove conflicting info */
985 2 : if ( bCopyBandInfo ) {
986 2 : pszSrcInfo = GDALGetDescription( hSrcBand );
987 2 : pszDstInfo = GDALGetDescription( hDstBand );
988 2 : if( ! ( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 &&
989 : pszDstInfo != NULL && strlen(pszDstInfo) > 0 &&
990 : EQUAL( pszSrcInfo, pszDstInfo ) ) )
991 2 : GDALSetDescription( hDstBand, "" );
992 2 : pszSrcInfo = GDALGetRasterUnitType( hSrcBand );
993 2 : pszDstInfo = GDALGetRasterUnitType( hDstBand );
994 2 : if( ! ( pszSrcInfo != NULL && strlen(pszSrcInfo) > 0 &&
995 : pszDstInfo != NULL && strlen(pszDstInfo) > 0 &&
996 : EQUAL( pszSrcInfo, pszDstInfo ) ) )
997 2 : 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 540 : 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 540 : if( GDALGetRasterColorInterpretation(
1023 : GDALGetRasterBand(hSrcDS,GDALGetRasterCount(hSrcDS)) )
1024 : == GCI_AlphaBand
1025 : && !bEnableSrcAlpha )
1026 : {
1027 2 : bEnableSrcAlpha = TRUE;
1028 2 : if( !bQuiet )
1029 : printf( "Using band %d of source image as alpha.\n",
1030 2 : GDALGetRasterCount(hSrcDS) );
1031 : }
1032 :
1033 : /* -------------------------------------------------------------------- */
1034 : /* Create a transformation object from the source to */
1035 : /* destination coordinate system. */
1036 : /* -------------------------------------------------------------------- */
1037 540 : if (hUniqueTransformArg)
1038 534 : hTransformArg = hGenImgProjArg = hUniqueTransformArg;
1039 : else
1040 : hTransformArg = hGenImgProjArg =
1041 6 : GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, papszTO );
1042 :
1043 540 : if( hTransformArg == NULL )
1044 0 : GDALExit( 1 );
1045 :
1046 540 : pfnTransformer = GDALGenImgProjTransform;
1047 :
1048 : /* -------------------------------------------------------------------- */
1049 : /* Warp the transformer with a linear approximator unless the */
1050 : /* acceptable error is zero. */
1051 : /* -------------------------------------------------------------------- */
1052 540 : if( dfErrorThreshold != 0.0 )
1053 : {
1054 : hTransformArg = hApproxArg =
1055 : GDALCreateApproxTransformer( GDALGenImgProjTransform,
1056 536 : hGenImgProjArg, dfErrorThreshold);
1057 536 : pfnTransformer = GDALApproxTransform;
1058 : }
1059 :
1060 : /* -------------------------------------------------------------------- */
1061 : /* Clear temporary INIT_DEST settings after the first image. */
1062 : /* -------------------------------------------------------------------- */
1063 540 : if( bInitDestSetForFirst && iSrc == 1 )
1064 : papszWarpOptions = CSLSetNameValue( papszWarpOptions,
1065 2 : "INIT_DEST", NULL );
1066 :
1067 : /* -------------------------------------------------------------------- */
1068 : /* Setup warp options. */
1069 : /* -------------------------------------------------------------------- */
1070 540 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
1071 :
1072 540 : psWO->papszWarpOptions = CSLDuplicate(papszWarpOptions);
1073 540 : psWO->eWorkingDataType = eWorkingType;
1074 540 : psWO->eResampleAlg = eResampleAlg;
1075 :
1076 540 : psWO->hSrcDS = hSrcDS;
1077 540 : psWO->hDstDS = hDstDS;
1078 :
1079 540 : psWO->pfnTransformer = pfnTransformer;
1080 540 : psWO->pTransformerArg = hTransformArg;
1081 :
1082 540 : if( !bQuiet )
1083 538 : psWO->pfnProgress = GDALTermProgress;
1084 :
1085 540 : if( dfWarpMemoryLimit != 0.0 )
1086 18 : psWO->dfWarpMemoryLimit = dfWarpMemoryLimit;
1087 :
1088 : /* -------------------------------------------------------------------- */
1089 : /* Setup band mapping. */
1090 : /* -------------------------------------------------------------------- */
1091 540 : if( bEnableSrcAlpha )
1092 2 : psWO->nBandCount = GDALGetRasterCount(hSrcDS) - 1;
1093 : else
1094 538 : psWO->nBandCount = GDALGetRasterCount(hSrcDS);
1095 :
1096 540 : psWO->panSrcBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
1097 540 : psWO->panDstBands = (int *) CPLMalloc(psWO->nBandCount*sizeof(int));
1098 :
1099 1142 : for( i = 0; i < psWO->nBandCount; i++ )
1100 : {
1101 602 : psWO->panSrcBands[i] = i+1;
1102 602 : psWO->panDstBands[i] = i+1;
1103 : }
1104 :
1105 : /* -------------------------------------------------------------------- */
1106 : /* Setup alpha bands used if any. */
1107 : /* -------------------------------------------------------------------- */
1108 540 : if( bEnableSrcAlpha )
1109 2 : psWO->nSrcAlphaBand = GDALGetRasterCount(hSrcDS);
1110 :
1111 540 : if( !bEnableDstAlpha
1112 : && GDALGetRasterCount(hDstDS) == psWO->nBandCount+1
1113 : && GDALGetRasterColorInterpretation(
1114 : GDALGetRasterBand(hDstDS,GDALGetRasterCount(hDstDS)))
1115 : == GCI_AlphaBand )
1116 : {
1117 2 : if( !bQuiet )
1118 : printf( "Using band %d of destination image as alpha.\n",
1119 2 : GDALGetRasterCount(hDstDS) );
1120 :
1121 2 : bEnableDstAlpha = TRUE;
1122 : }
1123 :
1124 540 : if( bEnableDstAlpha )
1125 6 : psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
1126 :
1127 : /* -------------------------------------------------------------------- */
1128 : /* Setup NODATA options. */
1129 : /* -------------------------------------------------------------------- */
1130 540 : if( pszSrcNodata != NULL && !EQUALN(pszSrcNodata,"n",1) )
1131 : {
1132 2 : char **papszTokens = CSLTokenizeString( pszSrcNodata );
1133 2 : int nTokenCount = CSLCount(papszTokens);
1134 :
1135 : psWO->padfSrcNoDataReal = (double *)
1136 2 : CPLMalloc(psWO->nBandCount*sizeof(double));
1137 : psWO->padfSrcNoDataImag = (double *)
1138 2 : CPLMalloc(psWO->nBandCount*sizeof(double));
1139 :
1140 8 : for( i = 0; i < psWO->nBandCount; i++ )
1141 : {
1142 6 : if( i < nTokenCount )
1143 : {
1144 2 : CPLStringToComplex( papszTokens[i],
1145 : psWO->padfSrcNoDataReal + i,
1146 4 : psWO->padfSrcNoDataImag + i );
1147 : }
1148 : else
1149 : {
1150 4 : psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i-1];
1151 4 : psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i-1];
1152 : }
1153 : }
1154 :
1155 2 : CSLDestroy( papszTokens );
1156 :
1157 : psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
1158 2 : "UNIFIED_SRC_NODATA", "YES" );
1159 : }
1160 :
1161 : /* -------------------------------------------------------------------- */
1162 : /* If -srcnodata was not specified, but the data has nodata */
1163 : /* values, use them. */
1164 : /* -------------------------------------------------------------------- */
1165 540 : if( pszSrcNodata == NULL )
1166 : {
1167 538 : int bHaveNodata = FALSE;
1168 538 : double dfReal = 0.0;
1169 :
1170 1116 : for( i = 0; !bHaveNodata && i < psWO->nBandCount; i++ )
1171 : {
1172 578 : GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
1173 578 : dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
1174 : }
1175 :
1176 538 : if( bHaveNodata )
1177 : {
1178 4 : if( !bQuiet )
1179 : {
1180 2 : 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 2 : dfReal, papszSrcFiles[iSrc] );
1186 : }
1187 : psWO->padfSrcNoDataReal = (double *)
1188 4 : CPLMalloc(psWO->nBandCount*sizeof(double));
1189 : psWO->padfSrcNoDataImag = (double *)
1190 4 : CPLMalloc(psWO->nBandCount*sizeof(double));
1191 :
1192 26 : for( i = 0; i < psWO->nBandCount; i++ )
1193 : {
1194 22 : GDALRasterBandH hBand = GDALGetRasterBand( hSrcDS, i+1 );
1195 :
1196 22 : dfReal = GDALGetRasterNoDataValue( hBand, &bHaveNodata );
1197 :
1198 22 : if( bHaveNodata )
1199 : {
1200 22 : psWO->padfSrcNoDataReal[i] = dfReal;
1201 22 : 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 540 : if( pszDstNodata != NULL )
1217 : {
1218 2 : char **papszTokens = CSLTokenizeString( pszDstNodata );
1219 2 : int nTokenCount = CSLCount(papszTokens);
1220 :
1221 : psWO->padfDstNoDataReal = (double *)
1222 2 : CPLMalloc(psWO->nBandCount*sizeof(double));
1223 : psWO->padfDstNoDataImag = (double *)
1224 2 : CPLMalloc(psWO->nBandCount*sizeof(double));
1225 :
1226 4 : for( i = 0; i < psWO->nBandCount; i++ )
1227 : {
1228 2 : if( i < nTokenCount )
1229 : {
1230 2 : CPLStringToComplex( papszTokens[i],
1231 : psWO->padfDstNoDataReal + i,
1232 4 : 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 2 : GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, i+1 );
1241 2 : 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 2 : switch(GDALGetRasterDataType(hBand))
1250 : {
1251 : case GDT_Byte:
1252 2 : CLAMP(psWO->padfDstNoDataReal[i], GByte,
1253 : 0.0, 255.0);
1254 2 : 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 2 : 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 2 : 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 2 : if( bCreateOutput )
1290 : {
1291 : GDALSetRasterNoDataValue(
1292 : GDALGetRasterBand( hDstDS, psWO->panDstBands[i] ),
1293 2 : psWO->padfDstNoDataReal[i] );
1294 : }
1295 : }
1296 :
1297 2 : 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 540 : if( hCutline != NULL )
1305 : {
1306 : TransformCutlineToSource( hSrcDS, hCutline,
1307 : &(psWO->papszWarpOptions),
1308 6 : 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 540 : if( bVRT )
1317 : {
1318 8 : if( GDALInitializeWarpedVRT( hDstDS, psWO ) != CE_None )
1319 0 : GDALExit( 1 );
1320 :
1321 8 : GDALClose( hDstDS );
1322 8 : 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 8 : if (pfnTransformer == GDALApproxTransform)
1328 : {
1329 6 : if( hGenImgProjArg != NULL )
1330 6 : GDALDestroyGenImgProjTransformer( hGenImgProjArg );
1331 : }
1332 :
1333 8 : GDALDestroyWarpOptions( psWO );
1334 :
1335 8 : CPLFree( pszDstFilename );
1336 8 : CSLDestroy( argv );
1337 8 : CSLDestroy( papszSrcFiles );
1338 8 : CSLDestroy( papszWarpOptions );
1339 8 : CSLDestroy( papszTO );
1340 :
1341 8 : GDALDumpOpenDatasets( stderr );
1342 :
1343 8 : GDALDestroyDriverManager();
1344 :
1345 8 : return 0;
1346 : }
1347 :
1348 : /* -------------------------------------------------------------------- */
1349 : /* Initialize and execute the warp. */
1350 : /* -------------------------------------------------------------------- */
1351 532 : GDALWarpOperation oWO;
1352 :
1353 532 : if( oWO.Initialize( psWO ) == CE_None )
1354 : {
1355 : CPLErr eErr;
1356 532 : if( bMulti )
1357 : eErr = oWO.ChunkAndWarpMulti( 0, 0,
1358 : GDALGetRasterXSize( hDstDS ),
1359 2 : GDALGetRasterYSize( hDstDS ) );
1360 : else
1361 : eErr = oWO.ChunkAndWarpImage( 0, 0,
1362 : GDALGetRasterXSize( hDstDS ),
1363 530 : GDALGetRasterYSize( hDstDS ) );
1364 532 : if (eErr != CE_None)
1365 0 : bHasGotErr = TRUE;
1366 : }
1367 :
1368 : /* -------------------------------------------------------------------- */
1369 : /* Cleanup */
1370 : /* -------------------------------------------------------------------- */
1371 532 : if( hApproxArg != NULL )
1372 530 : GDALDestroyApproxTransformer( hApproxArg );
1373 :
1374 532 : if( hGenImgProjArg != NULL )
1375 532 : GDALDestroyGenImgProjTransformer( hGenImgProjArg );
1376 :
1377 532 : GDALDestroyWarpOptions( psWO );
1378 :
1379 532 : GDALClose( hSrcDS );
1380 532 : }
1381 :
1382 : /* -------------------------------------------------------------------- */
1383 : /* Final Cleanup. */
1384 : /* -------------------------------------------------------------------- */
1385 530 : CPLErrorReset();
1386 530 : GDALFlushCache( hDstDS );
1387 530 : if( CPLGetLastErrorType() != CE_None )
1388 0 : bHasGotErr = TRUE;
1389 530 : GDALClose( hDstDS );
1390 :
1391 530 : CPLFree( pszDstFilename );
1392 530 : CSLDestroy( argv );
1393 530 : CSLDestroy( papszSrcFiles );
1394 530 : CSLDestroy( papszWarpOptions );
1395 530 : CSLDestroy( papszTO );
1396 :
1397 530 : GDALDumpOpenDatasets( stderr );
1398 :
1399 530 : GDALDestroyDriverManager();
1400 :
1401 : #ifdef OGR_ENABLED
1402 530 : if( hCutline != NULL )
1403 6 : OGR_G_DestroyGeometry( (OGRGeometryH) hCutline );
1404 530 : OGRCleanupAll();
1405 : #endif
1406 :
1407 530 : 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 536 : 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 536 : GDALColorTableH hCT = NULL;
1434 536 : double dfWrkMinX=0, dfWrkMaxX=0, dfWrkMinY=0, dfWrkMaxY=0;
1435 536 : double dfWrkResX=0, dfWrkResY=0;
1436 536 : int nDstBandCount = 0;
1437 536 : std::vector<GDALColorInterp> apeColorInterpretations;
1438 :
1439 536 : *phTransformArg = NULL;
1440 536 : *phSrcDS = NULL;
1441 :
1442 : /* -------------------------------------------------------------------- */
1443 : /* Find the output driver. */
1444 : /* -------------------------------------------------------------------- */
1445 536 : hDriver = GDALGetDriverByName( pszFormat );
1446 536 : if( hDriver == NULL
1447 : || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
1448 : {
1449 : int iDr;
1450 :
1451 : printf( "Output driver `%s' not recognised or does not support\n",
1452 0 : pszFormat );
1453 : printf( "direct output file creation. The following format drivers are configured\n"
1454 0 : "and support direct output:\n" );
1455 :
1456 0 : for( iDr = 0; iDr < GDALGetDriverCount(); iDr++ )
1457 : {
1458 0 : GDALDriverH hDriver = GDALGetDriver(iDr);
1459 :
1460 0 : if( GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL) != NULL )
1461 : {
1462 : printf( " %s: %s\n",
1463 : GDALGetDriverShortName( hDriver ),
1464 0 : GDALGetDriverLongName( hDriver ) );
1465 : }
1466 : }
1467 0 : printf( "\n" );
1468 0 : GDALExit( 1 );
1469 : }
1470 :
1471 : /* -------------------------------------------------------------------- */
1472 : /* For virtual output files, we have to set a special subclass */
1473 : /* of dataset to create. */
1474 : /* -------------------------------------------------------------------- */
1475 536 : if( bVRT )
1476 : *ppapszCreateOptions =
1477 : CSLSetNameValue( *ppapszCreateOptions, "SUBCLASS",
1478 8 : "VRTWarpedDataset" );
1479 :
1480 : /* -------------------------------------------------------------------- */
1481 : /* Loop over all input files to collect extents. */
1482 : /* -------------------------------------------------------------------- */
1483 : int iSrc;
1484 536 : char *pszThisTargetSRS = (char*)CSLFetchNameValue( papszTO, "DST_SRS" );
1485 536 : if( pszThisTargetSRS != NULL )
1486 24 : pszThisTargetSRS = CPLStrdup( pszThisTargetSRS );
1487 :
1488 1074 : for( iSrc = 0; papszSrcFiles[iSrc] != NULL; iSrc++ )
1489 : {
1490 : GDALDatasetH hSrcDS;
1491 538 : const char *pszThisSourceSRS = CSLFetchNameValue(papszTO,"SRC_SRS");
1492 :
1493 538 : hSrcDS = GDALOpen( papszSrcFiles[iSrc], GA_ReadOnly );
1494 538 : if( hSrcDS == NULL )
1495 0 : GDALExit( 1 );
1496 :
1497 : /* -------------------------------------------------------------------- */
1498 : /* Check that there's at least one raster band */
1499 : /* -------------------------------------------------------------------- */
1500 538 : if ( GDALGetRasterCount(hSrcDS) == 0 )
1501 : {
1502 0 : fprintf(stderr, "Input file %s has no raster bands.\n", papszSrcFiles[iSrc] );
1503 0 : GDALExit( 1 );
1504 : }
1505 :
1506 538 : if( eDT == GDT_Unknown )
1507 534 : eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
1508 :
1509 : /* -------------------------------------------------------------------- */
1510 : /* If we are processing the first file, and it has a color */
1511 : /* table, then we will copy it to the destination file. */
1512 : /* -------------------------------------------------------------------- */
1513 538 : if( iSrc == 0 )
1514 : {
1515 536 : nDstBandCount = GDALGetRasterCount(hSrcDS);
1516 536 : hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
1517 536 : if( hCT != NULL )
1518 : {
1519 0 : hCT = GDALCloneColorTable( hCT );
1520 0 : if( !bQuiet )
1521 : printf( "Copying color table from %s to new file.\n",
1522 0 : papszSrcFiles[iSrc] );
1523 : }
1524 :
1525 1136 : for(int iBand = 0; iBand < nDstBandCount; iBand++)
1526 : {
1527 : apeColorInterpretations.push_back(
1528 600 : GDALGetRasterColorInterpretation(GDALGetRasterBand(hSrcDS,iBand+1)) );
1529 : }
1530 : }
1531 :
1532 : /* -------------------------------------------------------------------- */
1533 : /* Get the sourcesrs from the dataset, if not set already. */
1534 : /* -------------------------------------------------------------------- */
1535 538 : if( pszThisSourceSRS == NULL )
1536 : {
1537 538 : const char *pszMethod = CSLFetchNameValue( papszTO, "METHOD" );
1538 :
1539 538 : if( GDALGetProjectionRef( hSrcDS ) != NULL
1540 : && strlen(GDALGetProjectionRef( hSrcDS )) > 0
1541 : && (pszMethod == NULL || EQUAL(pszMethod,"GEOTRANSFORM")) )
1542 504 : pszThisSourceSRS = GDALGetProjectionRef( hSrcDS );
1543 :
1544 34 : else if( GDALGetGCPProjection( hSrcDS ) != NULL
1545 : && strlen(GDALGetGCPProjection(hSrcDS)) > 0
1546 : && GDALGetGCPCount( hSrcDS ) > 1
1547 : && (pszMethod == NULL || EQUALN(pszMethod,"GCP_",4)) )
1548 30 : pszThisSourceSRS = GDALGetGCPProjection( hSrcDS );
1549 4 : else if( pszMethod != NULL && EQUAL(pszMethod,"RPC") )
1550 0 : pszThisSourceSRS = SRS_WKT_WGS84;
1551 : else
1552 4 : pszThisSourceSRS = "";
1553 : }
1554 :
1555 538 : if( pszThisTargetSRS == NULL )
1556 512 : pszThisTargetSRS = CPLStrdup( pszThisSourceSRS );
1557 :
1558 : /* -------------------------------------------------------------------- */
1559 : /* Create a transformation object from the source to */
1560 : /* destination coordinate system. */
1561 : /* -------------------------------------------------------------------- */
1562 : hTransformArg =
1563 538 : GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
1564 :
1565 538 : if( hTransformArg == NULL )
1566 : {
1567 0 : CPLFree( pszThisTargetSRS );
1568 0 : GDALClose( hSrcDS );
1569 0 : return NULL;
1570 : }
1571 :
1572 538 : GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;
1573 :
1574 : /* -------------------------------------------------------------------- */
1575 : /* Get approximate output definition. */
1576 : /* -------------------------------------------------------------------- */
1577 : double adfThisGeoTransform[6];
1578 : double adfExtent[4];
1579 : int nThisPixels, nThisLines;
1580 :
1581 538 : if( GDALSuggestedWarpOutput2( hSrcDS,
1582 : psInfo->pfnTransform, hTransformArg,
1583 : adfThisGeoTransform,
1584 : &nThisPixels, &nThisLines,
1585 : adfExtent, 0 ) != CE_None )
1586 : {
1587 0 : CPLFree( pszThisTargetSRS );
1588 0 : GDALClose( hSrcDS );
1589 0 : return NULL;
1590 : }
1591 :
1592 538 : if (CPLGetConfigOption( "CHECK_WITH_INVERT_PROJ", NULL ) == NULL)
1593 : {
1594 538 : double MinX = adfExtent[0];
1595 538 : double MaxX = adfExtent[2];
1596 538 : double MaxY = adfExtent[3];
1597 538 : double MinY = adfExtent[1];
1598 538 : int bSuccess = TRUE;
1599 :
1600 : /* Check that the the edges of the target image are in the validity area */
1601 : /* of the target projection */
1602 : #define N_STEPS 20
1603 : int i,j;
1604 11636 : for(i=0;i<=N_STEPS && bSuccess;i++)
1605 : {
1606 243956 : for(j=0;j<=N_STEPS && bSuccess;j++)
1607 : {
1608 232858 : double dfRatioI = i * 1.0 / N_STEPS;
1609 232858 : double dfRatioJ = j * 1.0 / N_STEPS;
1610 232858 : double expected_x = (1 - dfRatioI) * MinX + dfRatioI * MaxX;
1611 232858 : double expected_y = (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
1612 232858 : double x = expected_x;
1613 232858 : double y = expected_y;
1614 232858 : double z = 0;
1615 : /* Target SRS coordinates to source image pixel coordinates */
1616 232858 : if (!psInfo->pfnTransform(hTransformArg, TRUE, 1, &x, &y, &z, &bSuccess) || !bSuccess)
1617 0 : bSuccess = FALSE;
1618 : /* Source image pixel coordinates to target SRS coordinates */
1619 232858 : if (!psInfo->pfnTransform(hTransformArg, FALSE, 1, &x, &y, &z, &bSuccess) || !bSuccess)
1620 0 : bSuccess = FALSE;
1621 232858 : if (fabs(x - expected_x) > (MaxX - MinX) / nThisPixels ||
1622 : fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
1623 10 : bSuccess = FALSE;
1624 : }
1625 : }
1626 :
1627 : /* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces ogrct.cpp */
1628 : /* to check the consistency of each requested projection result with the */
1629 : /* invert projection */
1630 538 : if (!bSuccess)
1631 : {
1632 10 : CPLSetConfigOption( "CHECK_WITH_INVERT_PROJ", "TRUE" );
1633 10 : CPLDebug("WARP", "Recompute out extent with CHECK_WITH_INVERT_PROJ=TRUE");
1634 :
1635 10 : if( GDALSuggestedWarpOutput2( hSrcDS,
1636 : psInfo->pfnTransform, hTransformArg,
1637 : adfThisGeoTransform,
1638 : &nThisPixels, &nThisLines,
1639 : adfExtent, 0 ) != CE_None )
1640 : {
1641 0 : CPLFree( pszThisTargetSRS );
1642 0 : GDALClose( hSrcDS );
1643 0 : return NULL;
1644 : }
1645 : }
1646 : }
1647 :
1648 : /* -------------------------------------------------------------------- */
1649 : /* Expand the working bounds to include this region, ensure the */
1650 : /* working resolution is no more than this resolution. */
1651 : /* -------------------------------------------------------------------- */
1652 1074 : if( dfWrkMaxX == 0.0 && dfWrkMinX == 0.0 )
1653 : {
1654 536 : dfWrkMinX = adfExtent[0];
1655 536 : dfWrkMaxX = adfExtent[2];
1656 536 : dfWrkMaxY = adfExtent[3];
1657 536 : dfWrkMinY = adfExtent[1];
1658 536 : dfWrkResX = adfThisGeoTransform[1];
1659 536 : dfWrkResY = ABS(adfThisGeoTransform[5]);
1660 : }
1661 : else
1662 : {
1663 2 : dfWrkMinX = MIN(dfWrkMinX,adfExtent[0]);
1664 2 : dfWrkMaxX = MAX(dfWrkMaxX,adfExtent[2]);
1665 2 : dfWrkMaxY = MAX(dfWrkMaxY,adfExtent[3]);
1666 2 : dfWrkMinY = MIN(dfWrkMinY,adfExtent[1]);
1667 2 : dfWrkResX = MIN(dfWrkResX,adfThisGeoTransform[1]);
1668 2 : dfWrkResY = MIN(dfWrkResY,ABS(adfThisGeoTransform[5]));
1669 : }
1670 :
1671 1072 : if (iSrc == 0 && papszSrcFiles[1] == NULL)
1672 : {
1673 534 : *phTransformArg = hTransformArg;
1674 534 : *phSrcDS = hSrcDS;
1675 : }
1676 : else
1677 : {
1678 4 : GDALDestroyGenImgProjTransformer( hTransformArg );
1679 4 : GDALClose( hSrcDS );
1680 : }
1681 : }
1682 :
1683 : /* -------------------------------------------------------------------- */
1684 : /* Did we have any usable sources? */
1685 : /* -------------------------------------------------------------------- */
1686 536 : if( nDstBandCount == 0 )
1687 : {
1688 : CPLError( CE_Failure, CPLE_AppDefined,
1689 0 : "No usable source images." );
1690 0 : CPLFree( pszThisTargetSRS );
1691 0 : return NULL;
1692 : }
1693 :
1694 : /* -------------------------------------------------------------------- */
1695 : /* Turn the suggested region into a geotransform and suggested */
1696 : /* number of pixels and lines. */
1697 : /* -------------------------------------------------------------------- */
1698 : double adfDstGeoTransform[6];
1699 : int nPixels, nLines;
1700 :
1701 536 : adfDstGeoTransform[0] = dfWrkMinX;
1702 536 : adfDstGeoTransform[1] = dfWrkResX;
1703 536 : adfDstGeoTransform[2] = 0.0;
1704 536 : adfDstGeoTransform[3] = dfWrkMaxY;
1705 536 : adfDstGeoTransform[4] = 0.0;
1706 536 : adfDstGeoTransform[5] = -1 * dfWrkResY;
1707 :
1708 536 : nPixels = (int) ((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
1709 536 : nLines = (int) ((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
1710 :
1711 : /* -------------------------------------------------------------------- */
1712 : /* Did the user override some parameters? */
1713 : /* -------------------------------------------------------------------- */
1714 542 : if( dfXRes != 0.0 && dfYRes != 0.0 )
1715 : {
1716 6 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1717 : {
1718 6 : dfMinX = adfDstGeoTransform[0];
1719 6 : dfMaxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
1720 6 : dfMaxY = adfDstGeoTransform[3];
1721 6 : dfMinY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
1722 : }
1723 :
1724 6 : if ( bTargetAlignedPixels )
1725 : {
1726 2 : dfMinX = floor(dfMinX / dfXRes) * dfXRes;
1727 2 : dfMaxX = ceil(dfMaxX / dfXRes) * dfXRes;
1728 2 : dfMinY = floor(dfMinY / dfYRes) * dfYRes;
1729 2 : dfMaxY = ceil(dfMaxY / dfYRes) * dfYRes;
1730 : }
1731 :
1732 6 : nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1733 6 : nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1734 6 : adfDstGeoTransform[0] = dfMinX;
1735 6 : adfDstGeoTransform[3] = dfMaxY;
1736 6 : adfDstGeoTransform[1] = dfXRes;
1737 6 : adfDstGeoTransform[5] = -dfYRes;
1738 : }
1739 :
1740 546 : else if( nForcePixels != 0 && nForceLines != 0 )
1741 : {
1742 16 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1743 : {
1744 16 : dfMinX = dfWrkMinX;
1745 16 : dfMaxX = dfWrkMaxX;
1746 16 : dfMaxY = dfWrkMaxY;
1747 16 : dfMinY = dfWrkMinY;
1748 : }
1749 :
1750 16 : dfXRes = (dfMaxX - dfMinX) / nForcePixels;
1751 16 : dfYRes = (dfMaxY - dfMinY) / nForceLines;
1752 :
1753 16 : adfDstGeoTransform[0] = dfMinX;
1754 16 : adfDstGeoTransform[3] = dfMaxY;
1755 16 : adfDstGeoTransform[1] = dfXRes;
1756 16 : adfDstGeoTransform[5] = -dfYRes;
1757 :
1758 16 : nPixels = nForcePixels;
1759 16 : nLines = nForceLines;
1760 : }
1761 :
1762 514 : else if( nForcePixels != 0 )
1763 : {
1764 0 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1765 : {
1766 0 : dfMinX = dfWrkMinX;
1767 0 : dfMaxX = dfWrkMaxX;
1768 0 : dfMaxY = dfWrkMaxY;
1769 0 : dfMinY = dfWrkMinY;
1770 : }
1771 :
1772 0 : dfXRes = (dfMaxX - dfMinX) / nForcePixels;
1773 0 : dfYRes = dfXRes;
1774 :
1775 0 : adfDstGeoTransform[0] = dfMinX;
1776 0 : adfDstGeoTransform[3] = dfMaxY;
1777 0 : adfDstGeoTransform[1] = dfXRes;
1778 0 : adfDstGeoTransform[5] = -dfYRes;
1779 :
1780 0 : nPixels = nForcePixels;
1781 0 : nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1782 : }
1783 :
1784 514 : else if( nForceLines != 0 )
1785 : {
1786 0 : if( dfMinX == 0.0 && dfMinY == 0.0 && dfMaxX == 0.0 && dfMaxY == 0.0 )
1787 : {
1788 0 : dfMinX = dfWrkMinX;
1789 0 : dfMaxX = dfWrkMaxX;
1790 0 : dfMaxY = dfWrkMaxY;
1791 0 : dfMinY = dfWrkMinY;
1792 : }
1793 :
1794 0 : dfYRes = (dfMaxY - dfMinY) / nForceLines;
1795 0 : dfXRes = dfYRes;
1796 :
1797 0 : adfDstGeoTransform[0] = dfMinX;
1798 0 : adfDstGeoTransform[3] = dfMaxY;
1799 0 : adfDstGeoTransform[1] = dfXRes;
1800 0 : adfDstGeoTransform[5] = -dfYRes;
1801 :
1802 0 : nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1803 0 : nLines = nForceLines;
1804 : }
1805 :
1806 514 : else if( dfMinX != 0.0 || dfMinY != 0.0 || dfMaxX != 0.0 || dfMaxY != 0.0 )
1807 : {
1808 2 : dfXRes = adfDstGeoTransform[1];
1809 2 : dfYRes = fabs(adfDstGeoTransform[5]);
1810 :
1811 2 : nPixels = (int) ((dfMaxX - dfMinX + (dfXRes/2.0)) / dfXRes);
1812 2 : nLines = (int) ((dfMaxY - dfMinY + (dfYRes/2.0)) / dfYRes);
1813 :
1814 2 : dfXRes = (dfMaxX - dfMinX) / nPixels;
1815 2 : dfYRes = (dfMaxY - dfMinY) / nLines;
1816 :
1817 2 : adfDstGeoTransform[0] = dfMinX;
1818 2 : adfDstGeoTransform[3] = dfMaxY;
1819 2 : adfDstGeoTransform[1] = dfXRes;
1820 2 : adfDstGeoTransform[5] = -dfYRes;
1821 : }
1822 :
1823 : /* -------------------------------------------------------------------- */
1824 : /* Do we want to generate an alpha band in the output file? */
1825 : /* -------------------------------------------------------------------- */
1826 536 : if( bEnableSrcAlpha )
1827 0 : nDstBandCount--;
1828 :
1829 536 : if( bEnableDstAlpha )
1830 4 : nDstBandCount++;
1831 :
1832 : /* -------------------------------------------------------------------- */
1833 : /* Create the output file. */
1834 : /* -------------------------------------------------------------------- */
1835 536 : if( !bQuiet )
1836 534 : printf( "Creating output file that is %dP x %dL.\n", nPixels, nLines );
1837 :
1838 : hDstDS = GDALCreate( hDriver, pszFilename, nPixels, nLines,
1839 536 : nDstBandCount, eDT, *ppapszCreateOptions );
1840 :
1841 536 : if( hDstDS == NULL )
1842 : {
1843 0 : CPLFree( pszThisTargetSRS );
1844 0 : return NULL;
1845 : }
1846 :
1847 : /* -------------------------------------------------------------------- */
1848 : /* Write out the projection definition. */
1849 : /* -------------------------------------------------------------------- */
1850 536 : GDALSetProjection( hDstDS, pszThisTargetSRS );
1851 536 : GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
1852 :
1853 536 : if (*phTransformArg != NULL)
1854 534 : GDALSetGenImgProjTransformerDstGeoTransform( *phTransformArg, adfDstGeoTransform);
1855 :
1856 : /* -------------------------------------------------------------------- */
1857 : /* Try to set color interpretation of source bands to target */
1858 : /* dataset. */
1859 : /* FIXME? We should likely do that for other drivers than VRT */
1860 : /* but it might create spurious .aux.xml files (at least with HFA, */
1861 : /* and netCDF) */
1862 : /* -------------------------------------------------------------------- */
1863 536 : if( bVRT )
1864 : {
1865 8 : int nBandsToCopy = (int)apeColorInterpretations.size();
1866 8 : if ( bEnableSrcAlpha )
1867 0 : nBandsToCopy --;
1868 30 : for(int iBand = 0; iBand < nBandsToCopy; iBand++)
1869 : {
1870 : GDALSetRasterColorInterpretation(
1871 : GDALGetRasterBand( hDstDS, iBand + 1 ),
1872 22 : apeColorInterpretations[iBand] );
1873 : }
1874 : }
1875 :
1876 : /* -------------------------------------------------------------------- */
1877 : /* Try to set color interpretation of output file alpha band. */
1878 : /* -------------------------------------------------------------------- */
1879 536 : if( bEnableDstAlpha )
1880 : {
1881 : GDALSetRasterColorInterpretation(
1882 : GDALGetRasterBand( hDstDS, nDstBandCount ),
1883 4 : GCI_AlphaBand );
1884 : }
1885 :
1886 : /* -------------------------------------------------------------------- */
1887 : /* Copy the color table, if required. */
1888 : /* -------------------------------------------------------------------- */
1889 536 : if( hCT != NULL )
1890 : {
1891 0 : GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
1892 0 : GDALDestroyColorTable( hCT );
1893 : }
1894 :
1895 536 : CPLFree( pszThisTargetSRS );
1896 536 : return hDstDS;
1897 : }
1898 :
1899 : /************************************************************************/
1900 : /* GeoTransform_Transformer() */
1901 : /* */
1902 : /* Convert points from georef coordinates to pixel/line based */
1903 : /* on a geotransform. */
1904 : /************************************************************************/
1905 :
1906 : class CutlineTransformer : public OGRCoordinateTransformation
1907 12 : {
1908 : public:
1909 :
1910 : void *hSrcImageTransformer;
1911 :
1912 0 : virtual OGRSpatialReference *GetSourceCS() { return NULL; }
1913 18 : virtual OGRSpatialReference *GetTargetCS() { return NULL; }
1914 :
1915 0 : virtual int Transform( int nCount,
1916 : double *x, double *y, double *z = NULL ) {
1917 : int nResult;
1918 :
1919 0 : int *pabSuccess = (int *) CPLCalloc(sizeof(int),nCount);
1920 0 : nResult = TransformEx( nCount, x, y, z, pabSuccess );
1921 0 : CPLFree( pabSuccess );
1922 :
1923 0 : return nResult;
1924 : }
1925 :
1926 6 : virtual int TransformEx( int nCount,
1927 : double *x, double *y, double *z = NULL,
1928 : int *pabSuccess = NULL ) {
1929 : return GDALGenImgProjTransform( hSrcImageTransformer, TRUE,
1930 6 : nCount, x, y, z, pabSuccess );
1931 : }
1932 : };
1933 :
1934 :
1935 : /************************************************************************/
1936 : /* LoadCutline() */
1937 : /* */
1938 : /* Load blend cutline from OGR datasource. */
1939 : /************************************************************************/
1940 :
1941 : static void
1942 6 : LoadCutline( const char *pszCutlineDSName, const char *pszCLayer,
1943 : const char *pszCWHERE, const char *pszCSQL,
1944 : void **phCutlineRet )
1945 :
1946 : {
1947 : #ifndef OGR_ENABLED
1948 : CPLError( CE_Failure, CPLE_AppDefined,
1949 : "Request to load a cutline failed, this build does not support OGR features.\n" );
1950 : GDALExit( 1 );
1951 : #else // def OGR_ENABLED
1952 6 : OGRRegisterAll();
1953 :
1954 : /* -------------------------------------------------------------------- */
1955 : /* Open source vector dataset. */
1956 : /* -------------------------------------------------------------------- */
1957 : OGRDataSourceH hSrcDS;
1958 :
1959 6 : hSrcDS = OGROpen( pszCutlineDSName, FALSE, NULL );
1960 6 : if( hSrcDS == NULL )
1961 0 : GDALExit( 1 );
1962 :
1963 : /* -------------------------------------------------------------------- */
1964 : /* Get the source layer */
1965 : /* -------------------------------------------------------------------- */
1966 6 : OGRLayerH hLayer = NULL;
1967 :
1968 6 : if( pszCSQL != NULL )
1969 0 : hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszCSQL, NULL, NULL );
1970 6 : else if( pszCLayer != NULL )
1971 6 : hLayer = OGR_DS_GetLayerByName( hSrcDS, pszCLayer );
1972 : else
1973 0 : hLayer = OGR_DS_GetLayer( hSrcDS, 0 );
1974 :
1975 6 : if( hLayer == NULL )
1976 : {
1977 0 : fprintf( stderr, "Failed to identify source layer from datasource.\n" );
1978 0 : GDALExit( 1 );
1979 : }
1980 :
1981 : /* -------------------------------------------------------------------- */
1982 : /* Apply WHERE clause if there is one. */
1983 : /* -------------------------------------------------------------------- */
1984 6 : if( pszCWHERE != NULL )
1985 0 : OGR_L_SetAttributeFilter( hLayer, pszCWHERE );
1986 :
1987 : /* -------------------------------------------------------------------- */
1988 : /* Collect the geometries from this layer, and build list of */
1989 : /* burn values. */
1990 : /* -------------------------------------------------------------------- */
1991 : OGRFeatureH hFeat;
1992 6 : OGRGeometryH hMultiPolygon = OGR_G_CreateGeometry( wkbMultiPolygon );
1993 :
1994 6 : OGR_L_ResetReading( hLayer );
1995 :
1996 18 : while( (hFeat = OGR_L_GetNextFeature( hLayer )) != NULL )
1997 : {
1998 6 : OGRGeometryH hGeom = OGR_F_GetGeometryRef(hFeat);
1999 :
2000 6 : if( hGeom == NULL )
2001 : {
2002 0 : fprintf( stderr, "ERROR: Cutline feature without a geometry.\n" );
2003 0 : GDALExit( 1 );
2004 : }
2005 :
2006 6 : OGRwkbGeometryType eType = wkbFlatten(OGR_G_GetGeometryType( hGeom ));
2007 :
2008 6 : if( eType == wkbPolygon )
2009 6 : OGR_G_AddGeometry( hMultiPolygon, hGeom );
2010 0 : else if( eType == wkbMultiPolygon )
2011 : {
2012 : int iGeom;
2013 :
2014 0 : for( iGeom = 0; iGeom < OGR_G_GetGeometryCount( hGeom ); iGeom++ )
2015 : {
2016 : OGR_G_AddGeometry( hMultiPolygon,
2017 0 : OGR_G_GetGeometryRef(hGeom,iGeom) );
2018 : }
2019 : }
2020 : else
2021 : {
2022 0 : fprintf( stderr, "ERROR: Cutline not of polygon type.\n" );
2023 0 : GDALExit( 1 );
2024 : }
2025 :
2026 6 : OGR_F_Destroy( hFeat );
2027 : }
2028 :
2029 6 : if( OGR_G_GetGeometryCount( hMultiPolygon ) == 0 )
2030 : {
2031 0 : fprintf( stderr, "ERROR: Did not get any cutline features.\n" );
2032 0 : GDALExit( 1 );
2033 : }
2034 :
2035 : /* -------------------------------------------------------------------- */
2036 : /* Ensure the coordinate system gets set on the geometry. */
2037 : /* -------------------------------------------------------------------- */
2038 : OGR_G_AssignSpatialReference(
2039 6 : hMultiPolygon, OGR_L_GetSpatialRef(hLayer) );
2040 :
2041 6 : *phCutlineRet = (void *) hMultiPolygon;
2042 :
2043 : /* -------------------------------------------------------------------- */
2044 : /* Cleanup */
2045 : /* -------------------------------------------------------------------- */
2046 6 : if( pszCSQL != NULL )
2047 0 : OGR_DS_ReleaseResultSet( hSrcDS, hLayer );
2048 :
2049 6 : OGR_DS_Destroy( hSrcDS );
2050 : #endif
2051 6 : }
2052 :
2053 : /************************************************************************/
2054 : /* TransformCutlineToSource() */
2055 : /* */
2056 : /* Transform cutline from its SRS to source pixel/line coordinates.*/
2057 : /************************************************************************/
2058 : static void
2059 6 : TransformCutlineToSource( GDALDatasetH hSrcDS, void *hCutline,
2060 : char ***ppapszWarpOptions, char **papszTO_In )
2061 :
2062 : {
2063 : #ifdef OGR_ENABLED
2064 6 : OGRGeometryH hMultiPolygon = OGR_G_Clone( (OGRGeometryH) hCutline );
2065 6 : char **papszTO = CSLDuplicate( papszTO_In );
2066 :
2067 : /* -------------------------------------------------------------------- */
2068 : /* Checkout that SRS are the same. */
2069 : /* -------------------------------------------------------------------- */
2070 6 : OGRSpatialReferenceH hRasterSRS = NULL;
2071 6 : const char *pszProjection = NULL;
2072 :
2073 6 : if( GDALGetProjectionRef( hSrcDS ) != NULL
2074 : && strlen(GDALGetProjectionRef( hSrcDS )) > 0 )
2075 6 : pszProjection = GDALGetProjectionRef( hSrcDS );
2076 0 : else if( GDALGetGCPProjection( hSrcDS ) != NULL )
2077 0 : pszProjection = GDALGetGCPProjection( hSrcDS );
2078 :
2079 6 : if( pszProjection != NULL )
2080 : {
2081 6 : hRasterSRS = OSRNewSpatialReference(NULL);
2082 6 : if( OSRImportFromWkt( hRasterSRS, (char **)&pszProjection ) != CE_None )
2083 : {
2084 0 : OSRDestroySpatialReference(hRasterSRS);
2085 0 : hRasterSRS = NULL;
2086 : }
2087 : }
2088 :
2089 6 : OGRSpatialReferenceH hCutlineSRS = OGR_G_GetSpatialReference( hMultiPolygon );
2090 6 : if( hRasterSRS != NULL && hCutlineSRS != NULL )
2091 : {
2092 : /* ok, we will reproject */
2093 : }
2094 0 : else if( hRasterSRS != NULL && hCutlineSRS == NULL )
2095 : {
2096 : fprintf(stderr,
2097 : "Warning : the source raster dataset has a SRS, but the cutline features\n"
2098 : "not. We assume that the cutline coordinates are expressed in the destination SRS.\n"
2099 0 : "If not, cutline results may be incorrect.\n");
2100 : }
2101 0 : else if( hRasterSRS == NULL && hCutlineSRS != NULL )
2102 : {
2103 : fprintf(stderr,
2104 : "Warning : the input vector layer has a SRS, but the source raster dataset does not.\n"
2105 0 : "Cutline results may be incorrect.\n");
2106 : }
2107 :
2108 6 : if( hRasterSRS != NULL )
2109 6 : OSRDestroySpatialReference(hRasterSRS);
2110 :
2111 : /* -------------------------------------------------------------------- */
2112 : /* Extract the cutline SRS WKT. */
2113 : /* -------------------------------------------------------------------- */
2114 6 : if( hCutlineSRS != NULL )
2115 : {
2116 6 : char *pszCutlineSRS_WKT = NULL;
2117 :
2118 6 : OSRExportToWkt( hCutlineSRS, &pszCutlineSRS_WKT );
2119 6 : papszTO = CSLSetNameValue( papszTO, "DST_SRS", pszCutlineSRS_WKT );
2120 6 : CPLFree( pszCutlineSRS_WKT );
2121 : }
2122 :
2123 : /* -------------------------------------------------------------------- */
2124 : /* It may be unwise to let the mask geometry be re-wrapped by */
2125 : /* the CENTER_LONG machinery as this can easily screw up world */
2126 : /* spanning masks and invert the mask topology. */
2127 : /* -------------------------------------------------------------------- */
2128 6 : papszTO = CSLSetNameValue( papszTO, "INSERT_CENTER_LONG", "FALSE" );
2129 :
2130 : /* -------------------------------------------------------------------- */
2131 : /* Transform the geometry to pixel/line coordinates. */
2132 : /* -------------------------------------------------------------------- */
2133 6 : CutlineTransformer oTransformer;
2134 :
2135 : /* The cutline transformer will *invert* the hSrcImageTransformer */
2136 : /* so it will convert from the cutline SRS to the source pixel/line */
2137 : /* coordinates */
2138 : oTransformer.hSrcImageTransformer =
2139 6 : GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszTO );
2140 :
2141 6 : CSLDestroy( papszTO );
2142 :
2143 6 : if( oTransformer.hSrcImageTransformer == NULL )
2144 0 : GDALExit( 1 );
2145 :
2146 : OGR_G_Transform( hMultiPolygon,
2147 6 : (OGRCoordinateTransformationH) &oTransformer );
2148 :
2149 6 : GDALDestroyGenImgProjTransformer( oTransformer.hSrcImageTransformer );
2150 :
2151 : /* -------------------------------------------------------------------- */
2152 : /* Convert aggregate geometry into WKT. */
2153 : /* -------------------------------------------------------------------- */
2154 6 : char *pszWKT = NULL;
2155 :
2156 6 : OGR_G_ExportToWkt( hMultiPolygon, &pszWKT );
2157 6 : OGR_G_DestroyGeometry( hMultiPolygon );
2158 :
2159 : *ppapszWarpOptions = CSLSetNameValue( *ppapszWarpOptions,
2160 6 : "CUTLINE", pszWKT );
2161 6 : CPLFree( pszWKT );
2162 : #endif
2163 6 : }
2164 :
2165 : static void
2166 4 : RemoveConflictingMetadata( GDALMajorObjectH hObj, char **papszMetadata,
2167 : const char *pszValueConflict )
2168 : {
2169 4 : if ( hObj == NULL ) return;
2170 :
2171 4 : char *pszKey = NULL;
2172 : const char *pszValueRef;
2173 : const char *pszValueComp;
2174 4 : char ** papszMetadataRef = CSLDuplicate( papszMetadata );
2175 4 : int nCount = CSLCount( papszMetadataRef );
2176 :
2177 4 : for( int i = 0; i < nCount; i++ )
2178 : {
2179 0 : pszValueRef = CPLParseNameValue( papszMetadataRef[i], &pszKey );
2180 0 : pszValueComp = GDALGetMetadataItem( hObj, pszKey, NULL );
2181 0 : if ( ( pszValueRef == NULL || pszValueComp == NULL ||
2182 : ! EQUAL( pszValueRef, pszValueComp ) ) &&
2183 : ! EQUAL( pszValueComp, pszValueConflict ) ) {
2184 0 : GDALSetMetadataItem( hObj, pszKey, pszValueConflict, NULL );
2185 : }
2186 0 : CPLFree( pszKey );
2187 : }
2188 :
2189 4 : CSLDestroy( papszMetadataRef );
2190 : }
|
