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