Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: High Performance Image Reprojector
4 : * Purpose: Test program for high performance warper API.
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2002, i3 - information integration and imaging
9 : * Fort Collins, CO
10 : * Copyright (c) 2007-2015, Even Rouault <even dot rouault at spatialys.com>
11 : * Copyright (c) 2015, Faza Mahamood
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include "cpl_port.h"
17 : #include "gdal_utils.h"
18 : #include "gdal_utils_priv.h"
19 : #include "gdalargumentparser.h"
20 :
21 : #include <cctype>
22 : #include <cmath>
23 : #include <cstdio>
24 : #include <cstdlib>
25 : #include <cstring>
26 :
27 : #include <algorithm>
28 : #include <array>
29 : #include <limits>
30 : #include <set>
31 : #include <utility>
32 : #include <vector>
33 :
34 : // Suppress deprecation warning for GDALOpenVerticalShiftGrid and
35 : // GDALApplyVerticalShiftGrid
36 : #ifndef CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid
37 : #define CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid(x)
38 : #define CPL_WARN_DEPRECATED_GDALApplyVerticalShiftGrid(x)
39 : #endif
40 :
41 : #include "commonutils.h"
42 : #include "cpl_conv.h"
43 : #include "cpl_error.h"
44 : #include "cpl_progress.h"
45 : #include "cpl_string.h"
46 : #include "gdal.h"
47 : #include "gdal_alg.h"
48 : #include "gdal_alg_priv.h"
49 : #include "gdal_priv.h"
50 : #include "gdalwarper.h"
51 : #include "ogr_api.h"
52 : #include "ogr_core.h"
53 : #include "ogr_geometry.h"
54 : #include "ogr_spatialref.h"
55 : #include "ogr_srs_api.h"
56 : #include "ogr_proj_p.h"
57 : #include "ogrct_priv.h"
58 : #include "ogrsf_frmts.h"
59 : #include "vrtdataset.h"
60 : #include "../frmts/gtiff/cogdriver.h"
61 :
62 : /************************************************************************/
63 : /* GDALWarpAppOptions */
64 : /************************************************************************/
65 :
66 : /** Options for use with GDALWarp(). GDALWarpAppOptions* must be allocated and
67 : * freed with GDALWarpAppOptionsNew() and GDALWarpAppOptionsFree() respectively.
68 : */
69 : struct GDALWarpAppOptions
70 : {
71 : /*! Raw program arguments */
72 : CPLStringList aosArgs{};
73 :
74 : /*! set georeferenced extents of output file to be created (in target SRS by
75 : default, or in the SRS specified with pszTE_SRS) */
76 : double dfMinX = 0;
77 : double dfMinY = 0;
78 : double dfMaxX = 0;
79 : double dfMaxY = 0;
80 :
81 : /*! the SRS in which to interpret the coordinates given in
82 : GDALWarpAppOptions::dfMinX, GDALWarpAppOptions::dfMinY,
83 : GDALWarpAppOptions::dfMaxX and GDALWarpAppOptions::dfMaxY. The SRS may be
84 : any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n or a file
85 : containing the WKT. It is a convenience e.g. when knowing the output
86 : coordinates in a geodetic long/lat SRS, but still wanting a result in a
87 : projected coordinate system. */
88 : std::string osTE_SRS{};
89 :
90 : /*! set output file resolution (in target georeferenced units) */
91 : double dfXRes = 0;
92 : double dfYRes = 0;
93 :
94 : /*! whether target pixels should have dfXRes == dfYRes */
95 : bool bSquarePixels = false;
96 :
97 : /*! align the coordinates of the extent of the output file to the values of
98 : the GDALWarpAppOptions::dfXRes and GDALWarpAppOptions::dfYRes, such that
99 : the aligned extent includes the minimum extent. */
100 : bool bTargetAlignedPixels = false;
101 :
102 : /*! set output file size in pixels and lines. If
103 : GDALWarpAppOptions::nForcePixels or GDALWarpAppOptions::nForceLines is
104 : set to 0, the other dimension will be guessed from the computed
105 : resolution. Note that GDALWarpAppOptions::nForcePixels and
106 : GDALWarpAppOptions::nForceLines cannot be used with
107 : GDALWarpAppOptions::dfXRes and GDALWarpAppOptions::dfYRes. */
108 : int nForcePixels = 0;
109 : int nForceLines = 0;
110 :
111 : /*! allow or suppress progress monitor and other non-error output */
112 : bool bQuiet = true;
113 :
114 : /*! the progress function to use */
115 : GDALProgressFunc pfnProgress = GDALDummyProgress;
116 :
117 : /*! pointer to the progress data variable */
118 : void *pProgressData = nullptr;
119 :
120 : /*! creates an output alpha band to identify nodata (unset/transparent)
121 : pixels when set to true */
122 : bool bEnableDstAlpha = false;
123 :
124 : /*! forces the last band of an input file to be considered as alpha band. */
125 : bool bEnableSrcAlpha = false;
126 :
127 : /*! Prevent a source alpha band from being considered as such */
128 : bool bDisableSrcAlpha = false;
129 :
130 : /*! output format. Use the short format name. */
131 : std::string osFormat{};
132 :
133 : bool bCreateOutput = false;
134 :
135 : /*! list of warp options. ("NAME1=VALUE1","NAME2=VALUE2",...). The
136 : GDALWarpOptions::aosWarpOptions docs show all options. */
137 : CPLStringList aosWarpOptions{};
138 :
139 : double dfErrorThreshold = -1;
140 :
141 : /*! the amount of memory (in megabytes) that the warp API is allowed
142 : to use for caching. */
143 : double dfWarpMemoryLimit = 0;
144 :
145 : /*! list of create options for the output format driver. See format
146 : specific documentation for legal creation options for each format. */
147 : CPLStringList aosCreateOptions{};
148 :
149 : /*! the data type of the output bands */
150 : GDALDataType eOutputType = GDT_Unknown;
151 :
152 : /*! working pixel data type. The data type of pixels in the source
153 : image and destination image buffers. */
154 : GDALDataType eWorkingType = GDT_Unknown;
155 :
156 : /*! the resampling method. Available methods are: near, bilinear,
157 : cubic, cubicspline, lanczos, average, mode, max, min, med,
158 : q1, q3, sum */
159 : GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;
160 :
161 : /*! whether -r was specified */
162 : bool bResampleAlgSpecifiedByUser = false;
163 :
164 : /*! nodata masking values for input bands (different values can be supplied
165 : for each band). ("value1 value2 ..."). Masked values will not be used
166 : in interpolation. Use a value of "None" to ignore intrinsic nodata
167 : settings on the source dataset. */
168 : std::string osSrcNodata{};
169 :
170 : /*! nodata values for output bands (different values can be supplied for
171 : each band). ("value1 value2 ..."). New files will be initialized to
172 : this value and if possible the nodata value will be recorded in the
173 : output file. Use a value of "None" to ensure that nodata is not defined.
174 : If this argument is not used then nodata values will be copied from
175 : the source dataset. */
176 : std::string osDstNodata{};
177 :
178 : /*! use multithreaded warping implementation. Multiple threads will be used
179 : to process chunks of image and perform input/output operation
180 : simultaneously. */
181 : bool bMulti = false;
182 :
183 : /*! list of transformer options suitable to pass to
184 : GDALCreateGenImgProjTransformer2().
185 : ("NAME1=VALUE1","NAME2=VALUE2",...) */
186 : CPLStringList aosTransformerOptions{};
187 :
188 : /*! enable use of a blend cutline from a vector dataset name or a WKT
189 : * geometry
190 : */
191 : std::string osCutlineDSNameOrWKT{};
192 :
193 : /*! cutline SRS */
194 : std::string osCutlineSRS{};
195 :
196 : /*! the named layer to be selected from the cutline datasource */
197 : std::string osCLayer{};
198 :
199 : /*! restrict desired cutline features based on attribute query */
200 : std::string osCWHERE{};
201 :
202 : /*! SQL query to select the cutline features instead of from a layer
203 : with osCLayer */
204 : std::string osCSQL{};
205 :
206 : /*! crop the extent of the target dataset to the extent of the cutline */
207 : bool bCropToCutline = false;
208 :
209 : /*! copy dataset and band metadata will be copied from the first source
210 : dataset. Items that differ between source datasets will be set "*" (see
211 : GDALWarpAppOptions::pszMDConflictValue) */
212 : bool bCopyMetadata = true;
213 :
214 : /*! copy band information from the first source dataset */
215 : bool bCopyBandInfo = true;
216 :
217 : /*! value to set metadata items that conflict between source datasets
218 : (default is "*"). Use "" to remove conflicting items. */
219 : std::string osMDConflictValue = "*";
220 :
221 : /*! set the color interpretation of the bands of the target dataset from the
222 : * source dataset */
223 : bool bSetColorInterpretation = false;
224 :
225 : /*! overview level of source files to be used */
226 : int nOvLevel = OVR_LEVEL_AUTO;
227 :
228 : /*! Whether to enable vertical shift adjustment */
229 : bool bVShift = false;
230 :
231 : /*! Whether to disable vertical shift adjustment */
232 : bool bNoVShift = false;
233 :
234 : /*! Source bands */
235 : std::vector<int> anSrcBands{};
236 :
237 : /*! Destination bands */
238 : std::vector<int> anDstBands{};
239 :
240 : /*! Used when using a temporary TIFF file while warping */
241 : bool bDeleteOutputFileOnceCreated = false;
242 : };
243 :
244 : static CPLErr
245 : LoadCutline(const std::string &osCutlineDSNameOrWKT, const std::string &osSRS,
246 : const std::string &oszCLayer, const std::string &osCWHERE,
247 : const std::string &osCSQL, OGRGeometryH *phCutlineRet);
248 : static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,
249 : OGRGeometry *poCutline,
250 : char ***ppapszWarpOptions,
251 : CSLConstList papszTO);
252 :
253 : static GDALDatasetH GDALWarpCreateOutput(
254 : int nSrcCount, GDALDatasetH *pahSrcDS, const char *pszFilename,
255 : const char *pszFormat, char **papszTO, CSLConstList papszCreateOptions,
256 : GDALDataType eDT, GDALTransformerArgUniquePtr &hTransformArg,
257 : bool bSetColorInterpretation, GDALWarpAppOptions *psOptions,
258 : bool bUpdateTransformerWithDestGT);
259 :
260 : static void RemoveConflictingMetadata(GDALMajorObjectH hObj,
261 : CSLConstList papszMetadata,
262 : const char *pszValueConflict);
263 :
264 3 : static double GetAverageSegmentLength(const OGRGeometry *poGeom)
265 : {
266 3 : if (!poGeom)
267 0 : return 0;
268 3 : switch (wkbFlatten(poGeom->getGeometryType()))
269 : {
270 1 : case wkbLineString:
271 : {
272 1 : const auto *poLS = poGeom->toLineString();
273 1 : double dfSum = 0;
274 1 : const int nPoints = poLS->getNumPoints();
275 1 : if (nPoints == 0)
276 0 : return 0;
277 5 : for (int i = 0; i < nPoints - 1; i++)
278 : {
279 4 : double dfX1 = poLS->getX(i);
280 4 : double dfY1 = poLS->getY(i);
281 4 : double dfX2 = poLS->getX(i + 1);
282 4 : double dfY2 = poLS->getY(i + 1);
283 4 : double dfDX = dfX2 - dfX1;
284 4 : double dfDY = dfY2 - dfY1;
285 4 : dfSum += sqrt(dfDX * dfDX + dfDY * dfDY);
286 : }
287 1 : return dfSum / nPoints;
288 : }
289 :
290 1 : case wkbPolygon:
291 : {
292 1 : if (poGeom->IsEmpty())
293 0 : return 0;
294 1 : double dfSum = 0;
295 2 : for (const auto *poLS : poGeom->toPolygon())
296 : {
297 1 : dfSum += GetAverageSegmentLength(poLS);
298 : }
299 1 : return dfSum / (1 + poGeom->toPolygon()->getNumInteriorRings());
300 : }
301 :
302 1 : case wkbMultiPolygon:
303 : case wkbMultiLineString:
304 : case wkbGeometryCollection:
305 : {
306 1 : if (poGeom->IsEmpty())
307 0 : return 0;
308 1 : double dfSum = 0;
309 2 : for (const auto *poSubGeom : poGeom->toGeometryCollection())
310 : {
311 1 : dfSum += GetAverageSegmentLength(poSubGeom);
312 : }
313 1 : return dfSum / poGeom->toGeometryCollection()->getNumGeometries();
314 : }
315 :
316 0 : default:
317 0 : return 0;
318 : }
319 : }
320 :
321 : /************************************************************************/
322 : /* FetchSrcMethod() */
323 : /************************************************************************/
324 :
325 3136 : static const char *FetchSrcMethod(CSLConstList papszTO,
326 : const char *pszDefault = nullptr)
327 : {
328 3136 : const char *pszMethod = CSLFetchNameValue(papszTO, "SRC_METHOD");
329 3136 : if (!pszMethod)
330 3060 : pszMethod = CSLFetchNameValueDef(papszTO, "METHOD", pszDefault);
331 3136 : return pszMethod;
332 : }
333 :
334 1185 : static const char *FetchSrcMethod(const CPLStringList &aosTO,
335 : const char *pszDefault = nullptr)
336 : {
337 1185 : const char *pszMethod = aosTO.FetchNameValue("SRC_METHOD");
338 1185 : if (!pszMethod)
339 1151 : pszMethod = aosTO.FetchNameValueDef("METHOD", pszDefault);
340 1185 : return pszMethod;
341 : }
342 :
343 : /************************************************************************/
344 : /* GetSrcDSProjection() */
345 : /* */
346 : /* Takes into account SRC_SRS transformer option in priority, and then */
347 : /* dataset characteristics as well as the METHOD transformer */
348 : /* option to determine the source SRS. */
349 : /************************************************************************/
350 :
351 1327 : static CPLString GetSrcDSProjection(GDALDatasetH hDS, CSLConstList papszTO)
352 : {
353 1327 : const char *pszProjection = CSLFetchNameValue(papszTO, "SRC_SRS");
354 1327 : if (pszProjection != nullptr || hDS == nullptr)
355 : {
356 65 : return pszProjection ? pszProjection : "";
357 : }
358 :
359 1262 : const char *pszMethod = FetchSrcMethod(papszTO);
360 1262 : CSLConstList papszMD = nullptr;
361 1262 : const OGRSpatialReferenceH hSRS = GDALGetSpatialRef(hDS);
362 : const char *pszGeolocationDataset =
363 1262 : CSLFetchNameValueDef(papszTO, "SRC_GEOLOC_ARRAY",
364 : CSLFetchNameValue(papszTO, "GEOLOC_ARRAY"));
365 1262 : if (pszGeolocationDataset != nullptr &&
366 2 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")))
367 : {
368 : auto aosMD =
369 3 : GDALCreateGeolocationMetadata(hDS, pszGeolocationDataset, true);
370 3 : pszProjection = aosMD.FetchNameValue("SRS");
371 3 : if (pszProjection)
372 3 : return pszProjection; // return in this scope so that aosMD is
373 : // still valid
374 : }
375 1259 : else if (hSRS && (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")))
376 : {
377 983 : char *pszWKT = nullptr;
378 : {
379 1966 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
380 983 : if (OSRExportToWkt(hSRS, &pszWKT) != OGRERR_NONE)
381 : {
382 0 : CPLFree(pszWKT);
383 0 : pszWKT = nullptr;
384 0 : const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};
385 0 : OSRExportToWktEx(hSRS, &pszWKT, apszOptions);
386 : }
387 : }
388 1966 : CPLString osWKT = pszWKT ? pszWKT : "";
389 983 : CPLFree(pszWKT);
390 983 : return osWKT;
391 : }
392 276 : else if (GDALGetGCPProjection(hDS) != nullptr &&
393 276 : strlen(GDALGetGCPProjection(hDS)) > 0 &&
394 557 : GDALGetGCPCount(hDS) > 1 &&
395 5 : (pszMethod == nullptr || STARTS_WITH_CI(pszMethod, "GCP_")))
396 : {
397 27 : pszProjection = GDALGetGCPProjection(hDS);
398 : }
399 252 : else if (GDALGetMetadata(hDS, "RPC") != nullptr &&
400 3 : (pszMethod == nullptr || EQUAL(pszMethod, "RPC")))
401 : {
402 14 : pszProjection = SRS_WKT_WGS84_LAT_LONG;
403 : }
404 248 : else if ((papszMD = GDALGetMetadata(hDS, "GEOLOCATION")) != nullptr &&
405 13 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")))
406 : {
407 22 : pszProjection = CSLFetchNameValue(papszMD, "SRS");
408 : }
409 276 : return pszProjection ? pszProjection : "";
410 : }
411 :
412 : /************************************************************************/
413 : /* CreateCTCutlineToSrc() */
414 : /************************************************************************/
415 :
416 58 : static std::unique_ptr<OGRCoordinateTransformation> CreateCTCutlineToSrc(
417 : const OGRSpatialReference *poRasterSRS, const OGRSpatialReference *poDstSRS,
418 : const OGRSpatialReference *poCutlineSRS, CSLConstList papszTO)
419 : {
420 58 : const OGRSpatialReference *poCutlineOrTargetSRS =
421 58 : poCutlineSRS ? poCutlineSRS : poDstSRS;
422 58 : std::unique_ptr<OGRCoordinateTransformation> poCTCutlineToSrc;
423 112 : if (poCutlineOrTargetSRS && poRasterSRS &&
424 54 : !poCutlineOrTargetSRS->IsSame(poRasterSRS))
425 : {
426 16 : OGRCoordinateTransformationOptions oOptions;
427 : // If the cutline SRS is the same as the target SRS and there is
428 : // an explicit -ct between the source SRS and the target SRS, then
429 : // use it in the reverse way to transform from the cutline SRS to
430 : // the source SRS.
431 8 : if (poDstSRS && poCutlineOrTargetSRS->IsSame(poDstSRS))
432 : {
433 : const char *pszCT =
434 2 : CSLFetchNameValue(papszTO, "COORDINATE_OPERATION");
435 2 : if (pszCT)
436 : {
437 1 : oOptions.SetCoordinateOperation(pszCT, /* bInverse = */ true);
438 : }
439 : }
440 8 : poCTCutlineToSrc.reset(OGRCreateCoordinateTransformation(
441 : poCutlineOrTargetSRS, poRasterSRS, oOptions));
442 : }
443 58 : return poCTCutlineToSrc;
444 : }
445 :
446 : /************************************************************************/
447 : /* CropToCutline() */
448 : /************************************************************************/
449 :
450 18 : static CPLErr CropToCutline(const OGRGeometry *poCutline, CSLConstList papszTO,
451 : CSLConstList papszWarpOptions, int nSrcCount,
452 : GDALDatasetH *pahSrcDS, double &dfMinX,
453 : double &dfMinY, double &dfMaxX, double &dfMaxY,
454 : const GDALWarpAppOptions *psOptions)
455 : {
456 : // We could possibly directly reproject from cutline SRS to target SRS,
457 : // but when applying the cutline, it is reprojected to source raster image
458 : // space using the source SRS. To be consistent, we reproject
459 : // the cutline from cutline SRS to source SRS and then from source SRS to
460 : // target SRS.
461 18 : const OGRSpatialReference *poCutlineSRS = poCutline->getSpatialReference();
462 18 : const char *pszThisTargetSRS = CSLFetchNameValue(papszTO, "DST_SRS");
463 18 : std::unique_ptr<OGRSpatialReference> poSrcSRS;
464 18 : std::unique_ptr<OGRSpatialReference> poDstSRS;
465 :
466 : const CPLString osThisSourceSRS =
467 36 : GetSrcDSProjection(nSrcCount > 0 ? pahSrcDS[0] : nullptr, papszTO);
468 18 : if (!osThisSourceSRS.empty())
469 : {
470 14 : poSrcSRS = std::make_unique<OGRSpatialReference>();
471 14 : poSrcSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
472 14 : if (poSrcSRS->SetFromUserInput(osThisSourceSRS) != OGRERR_NONE)
473 : {
474 1 : CPLError(CE_Failure, CPLE_AppDefined,
475 : "Cannot compute bounding box of cutline.");
476 1 : return CE_Failure;
477 : }
478 : }
479 4 : else if (!pszThisTargetSRS && !poCutlineSRS)
480 : {
481 1 : OGREnvelope sEnvelope;
482 1 : poCutline->getEnvelope(&sEnvelope);
483 :
484 1 : dfMinX = sEnvelope.MinX;
485 1 : dfMinY = sEnvelope.MinY;
486 1 : dfMaxX = sEnvelope.MaxX;
487 1 : dfMaxY = sEnvelope.MaxY;
488 :
489 1 : return CE_None;
490 : }
491 : else
492 : {
493 3 : CPLError(CE_Failure, CPLE_AppDefined,
494 : "Cannot compute bounding box of cutline. Cannot find "
495 : "source SRS");
496 3 : return CE_Failure;
497 : }
498 :
499 13 : if (pszThisTargetSRS)
500 : {
501 3 : poDstSRS = std::make_unique<OGRSpatialReference>();
502 3 : poDstSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
503 3 : if (poDstSRS->SetFromUserInput(pszThisTargetSRS) != OGRERR_NONE)
504 : {
505 1 : CPLError(CE_Failure, CPLE_AppDefined,
506 : "Cannot compute bounding box of cutline.");
507 1 : return CE_Failure;
508 : }
509 : }
510 : else
511 : {
512 10 : poDstSRS.reset(poSrcSRS->Clone());
513 : }
514 :
515 24 : auto poCutlineGeom = std::unique_ptr<OGRGeometry>(poCutline->clone());
516 12 : auto poCTCutlineToSrc = CreateCTCutlineToSrc(poSrcSRS.get(), poDstSRS.get(),
517 24 : poCutlineSRS, papszTO);
518 :
519 12 : std::unique_ptr<OGRCoordinateTransformation> poCTSrcToDst;
520 12 : if (!poSrcSRS->IsSame(poDstSRS.get()))
521 : {
522 1 : poCTSrcToDst.reset(
523 1 : OGRCreateCoordinateTransformation(poSrcSRS.get(), poDstSRS.get()));
524 : }
525 :
526 : // Reproject cutline to target SRS, by doing intermediate vertex
527 : // densification in source SRS.
528 12 : if (poCTSrcToDst || poCTCutlineToSrc)
529 : {
530 2 : OGREnvelope sLastEnvelope, sCurEnvelope;
531 0 : std::unique_ptr<OGRGeometry> poTransformedGeom;
532 : auto poGeomInSrcSRS =
533 2 : std::unique_ptr<OGRGeometry>(poCutlineGeom->clone());
534 2 : if (poCTCutlineToSrc)
535 : {
536 4 : poGeomInSrcSRS.reset(OGRGeometryFactory::transformWithOptions(
537 2 : poGeomInSrcSRS.get(), poCTCutlineToSrc.get(), nullptr));
538 2 : if (!poGeomInSrcSRS)
539 0 : return CE_Failure;
540 : }
541 :
542 : // Do not use a smaller epsilon, otherwise it could cause useless
543 : // segmentization (https://github.com/OSGeo/gdal/issues/4826)
544 2 : constexpr double epsilon = 1e-10;
545 3 : for (int nIter = 0; nIter < 10; nIter++)
546 : {
547 3 : poTransformedGeom.reset(poGeomInSrcSRS->clone());
548 3 : if (poCTSrcToDst)
549 : {
550 4 : poTransformedGeom.reset(
551 : OGRGeometryFactory::transformWithOptions(
552 2 : poTransformedGeom.get(), poCTSrcToDst.get(), nullptr));
553 2 : if (!poTransformedGeom)
554 0 : return CE_Failure;
555 : }
556 3 : poTransformedGeom->getEnvelope(&sCurEnvelope);
557 3 : if (nIter > 0 || !poCTSrcToDst)
558 : {
559 2 : if (std::abs(sCurEnvelope.MinX - sLastEnvelope.MinX) <=
560 2 : epsilon *
561 2 : std::abs(sCurEnvelope.MinX + sLastEnvelope.MinX) &&
562 2 : std::abs(sCurEnvelope.MinY - sLastEnvelope.MinY) <=
563 2 : epsilon *
564 2 : std::abs(sCurEnvelope.MinY + sLastEnvelope.MinY) &&
565 2 : std::abs(sCurEnvelope.MaxX - sLastEnvelope.MaxX) <=
566 2 : epsilon *
567 6 : std::abs(sCurEnvelope.MaxX + sLastEnvelope.MaxX) &&
568 2 : std::abs(sCurEnvelope.MaxY - sLastEnvelope.MaxY) <=
569 2 : epsilon *
570 2 : std::abs(sCurEnvelope.MaxY + sLastEnvelope.MaxY))
571 : {
572 2 : break;
573 : }
574 : }
575 : double dfAverageSegmentLength =
576 1 : GetAverageSegmentLength(poGeomInSrcSRS.get());
577 1 : poGeomInSrcSRS->segmentize(dfAverageSegmentLength / 4);
578 :
579 1 : sLastEnvelope = sCurEnvelope;
580 : }
581 :
582 2 : poCutlineGeom = std::move(poTransformedGeom);
583 : }
584 :
585 12 : OGREnvelope sEnvelope;
586 12 : poCutlineGeom->getEnvelope(&sEnvelope);
587 :
588 12 : dfMinX = sEnvelope.MinX;
589 12 : dfMinY = sEnvelope.MinY;
590 12 : dfMaxX = sEnvelope.MaxX;
591 12 : dfMaxY = sEnvelope.MaxY;
592 22 : if (!poCTSrcToDst && nSrcCount > 0 && psOptions->dfXRes == 0.0 &&
593 10 : psOptions->dfYRes == 0.0)
594 : {
595 : // No raster reprojection: stick on exact pixel boundaries of the source
596 : // to preserve resolution and avoid resampling
597 : double adfGT[6];
598 10 : if (GDALGetGeoTransform(pahSrcDS[0], adfGT) == CE_None)
599 : {
600 : // We allow for a relative error in coordinates up to 0.1% of the
601 : // pixel size for rounding purposes.
602 10 : constexpr double REL_EPS_PIXEL = 1e-3;
603 10 : if (CPLFetchBool(papszWarpOptions, "CUTLINE_ALL_TOUCHED", false))
604 : {
605 : // All touched ? Then make the extent a bit larger than the
606 : // cutline envelope
607 3 : dfMinX = adfGT[0] +
608 3 : floor((dfMinX - adfGT[0]) / adfGT[1] + REL_EPS_PIXEL) *
609 3 : adfGT[1];
610 3 : dfMinY = adfGT[3] +
611 3 : ceil((dfMinY - adfGT[3]) / adfGT[5] - REL_EPS_PIXEL) *
612 3 : adfGT[5];
613 3 : dfMaxX = adfGT[0] +
614 3 : ceil((dfMaxX - adfGT[0]) / adfGT[1] - REL_EPS_PIXEL) *
615 3 : adfGT[1];
616 3 : dfMaxY = adfGT[3] +
617 3 : floor((dfMaxY - adfGT[3]) / adfGT[5] + REL_EPS_PIXEL) *
618 3 : adfGT[5];
619 : }
620 : else
621 : {
622 : // Otherwise, make it a bit smaller
623 7 : dfMinX = adfGT[0] +
624 7 : ceil((dfMinX - adfGT[0]) / adfGT[1] - REL_EPS_PIXEL) *
625 7 : adfGT[1];
626 7 : dfMinY = adfGT[3] +
627 7 : floor((dfMinY - adfGT[3]) / adfGT[5] + REL_EPS_PIXEL) *
628 7 : adfGT[5];
629 7 : dfMaxX = adfGT[0] +
630 7 : floor((dfMaxX - adfGT[0]) / adfGT[1] + REL_EPS_PIXEL) *
631 7 : adfGT[1];
632 7 : dfMaxY = adfGT[3] +
633 7 : ceil((dfMaxY - adfGT[3]) / adfGT[5] - REL_EPS_PIXEL) *
634 7 : adfGT[5];
635 : }
636 : }
637 : }
638 :
639 12 : return CE_None;
640 : }
641 :
642 2314 : static bool MustApplyVerticalShift(GDALDatasetH hWrkSrcDS,
643 : const GDALWarpAppOptions *psOptions,
644 : OGRSpatialReference &oSRSSrc,
645 : OGRSpatialReference &oSRSDst,
646 : bool &bSrcHasVertAxis, bool &bDstHasVertAxis)
647 : {
648 2314 : bool bApplyVShift = psOptions->bVShift;
649 :
650 : // Check if we must do vertical shift grid transform
651 : const char *pszSrcWKT =
652 2314 : psOptions->aosTransformerOptions.FetchNameValue("SRC_SRS");
653 2314 : if (pszSrcWKT)
654 62 : oSRSSrc.SetFromUserInput(pszSrcWKT);
655 : else
656 : {
657 2252 : auto hSRS = GDALGetSpatialRef(hWrkSrcDS);
658 2252 : if (hSRS)
659 1638 : oSRSSrc = *(OGRSpatialReference::FromHandle(hSRS));
660 : else
661 614 : return false;
662 : }
663 :
664 : const char *pszDstWKT =
665 1700 : psOptions->aosTransformerOptions.FetchNameValue("DST_SRS");
666 1700 : if (pszDstWKT)
667 545 : oSRSDst.SetFromUserInput(pszDstWKT);
668 : else
669 1155 : return false;
670 :
671 545 : if (oSRSSrc.IsSame(&oSRSDst))
672 20 : return false;
673 :
674 1020 : bSrcHasVertAxis = oSRSSrc.IsCompound() ||
675 495 : ((oSRSSrc.IsProjected() || oSRSSrc.IsGeographic()) &&
676 495 : oSRSSrc.GetAxesCount() == 3);
677 :
678 1028 : bDstHasVertAxis = oSRSDst.IsCompound() ||
679 503 : ((oSRSDst.IsProjected() || oSRSDst.IsGeographic()) &&
680 502 : oSRSDst.GetAxesCount() == 3);
681 :
682 877 : if ((GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&
683 352 : (bSrcHasVertAxis || bDstHasVertAxis))
684 : {
685 40 : bApplyVShift = true;
686 : }
687 525 : return bApplyVShift;
688 : }
689 :
690 : /************************************************************************/
691 : /* ApplyVerticalShift() */
692 : /************************************************************************/
693 :
694 1152 : static bool ApplyVerticalShift(GDALDatasetH hWrkSrcDS,
695 : const GDALWarpAppOptions *psOptions,
696 : GDALWarpOptions *psWO)
697 : {
698 1152 : if (psOptions->bVShift)
699 : {
700 0 : psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
701 : "APPLY_VERTICAL_SHIFT", "YES");
702 : }
703 :
704 2304 : OGRSpatialReference oSRSSrc;
705 1152 : OGRSpatialReference oSRSDst;
706 1152 : bool bSrcHasVertAxis = false;
707 1152 : bool bDstHasVertAxis = false;
708 : bool bApplyVShift =
709 1152 : MustApplyVerticalShift(hWrkSrcDS, psOptions, oSRSSrc, oSRSDst,
710 : bSrcHasVertAxis, bDstHasVertAxis);
711 :
712 2008 : if ((GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&
713 856 : (bSrcHasVertAxis || bDstHasVertAxis))
714 : {
715 20 : bApplyVShift = true;
716 20 : psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
717 : "APPLY_VERTICAL_SHIFT", "YES");
718 :
719 20 : if (CSLFetchNameValue(psWO->papszWarpOptions,
720 20 : "MULT_FACTOR_VERTICAL_SHIFT") == nullptr)
721 : {
722 : // Select how to go from input dataset units to meters
723 20 : double dfToMeterSrc = 1.0;
724 : const char *pszUnit =
725 20 : GDALGetRasterUnitType(GDALGetRasterBand(hWrkSrcDS, 1));
726 :
727 20 : double dfToMeterSrcAxis = 1.0;
728 20 : if (bSrcHasVertAxis)
729 : {
730 16 : oSRSSrc.GetAxis(nullptr, 2, nullptr, &dfToMeterSrcAxis);
731 : }
732 :
733 20 : if (pszUnit && (EQUAL(pszUnit, "m") || EQUAL(pszUnit, "meter") ||
734 19 : EQUAL(pszUnit, "metre")))
735 : {
736 : }
737 19 : else if (pszUnit &&
738 19 : (EQUAL(pszUnit, "ft") || EQUAL(pszUnit, "foot")))
739 : {
740 1 : dfToMeterSrc = CPLAtof(SRS_UL_FOOT_CONV);
741 : }
742 18 : else if (pszUnit && (EQUAL(pszUnit, "US survey foot")))
743 : {
744 2 : dfToMeterSrc = CPLAtof(SRS_UL_US_FOOT_CONV);
745 : }
746 16 : else if (pszUnit && !EQUAL(pszUnit, ""))
747 : {
748 1 : if (bSrcHasVertAxis)
749 : {
750 1 : dfToMeterSrc = dfToMeterSrcAxis;
751 : }
752 : else
753 : {
754 0 : CPLError(CE_Warning, CPLE_AppDefined,
755 : "Unknown units=%s. Assuming metre.", pszUnit);
756 : }
757 : }
758 : else
759 : {
760 15 : if (bSrcHasVertAxis)
761 11 : oSRSSrc.GetAxis(nullptr, 2, nullptr, &dfToMeterSrc);
762 : }
763 :
764 20 : double dfToMeterDst = 1.0;
765 20 : if (bDstHasVertAxis)
766 20 : oSRSDst.GetAxis(nullptr, 2, nullptr, &dfToMeterDst);
767 :
768 20 : if (dfToMeterSrc > 0 && dfToMeterDst > 0)
769 : {
770 20 : const double dfMultFactorVerticalShift =
771 20 : dfToMeterSrc / dfToMeterDst;
772 20 : CPLDebug("WARP", "Applying MULT_FACTOR_VERTICAL_SHIFT=%.18g",
773 : dfMultFactorVerticalShift);
774 20 : psWO->papszWarpOptions = CSLSetNameValue(
775 : psWO->papszWarpOptions, "MULT_FACTOR_VERTICAL_SHIFT",
776 : CPLSPrintf("%.18g", dfMultFactorVerticalShift));
777 :
778 20 : const double dfMultFactorVerticalShiftPipeline =
779 20 : dfToMeterSrcAxis / dfToMeterDst;
780 20 : CPLDebug("WARP",
781 : "Applying MULT_FACTOR_VERTICAL_SHIFT_PIPELINE=%.18g",
782 : dfMultFactorVerticalShiftPipeline);
783 20 : psWO->papszWarpOptions = CSLSetNameValue(
784 : psWO->papszWarpOptions,
785 : "MULT_FACTOR_VERTICAL_SHIFT_PIPELINE",
786 : CPLSPrintf("%.18g", dfMultFactorVerticalShiftPipeline));
787 : }
788 : }
789 : }
790 :
791 2304 : return bApplyVShift;
792 : }
793 :
794 : /************************************************************************/
795 : /* CanUseBuildVRT() */
796 : /************************************************************************/
797 :
798 5 : static bool CanUseBuildVRT(int nSrcCount, GDALDatasetH *pahSrcDS)
799 : {
800 :
801 5 : bool bCanUseBuildVRT = true;
802 10 : std::vector<std::array<double, 4>> aoExtents;
803 5 : bool bSrcHasAlpha = false;
804 5 : int nPrevBandCount = 0;
805 5 : OGRSpatialReference oSRSPrev;
806 5 : double dfLastResX = 0;
807 5 : double dfLastResY = 0;
808 14 : for (int i = 0; i < nSrcCount; i++)
809 : {
810 : double adfGT[6];
811 9 : auto hSrcDS = pahSrcDS[i];
812 9 : if (EQUAL(GDALGetDescription(hSrcDS), ""))
813 : {
814 0 : bCanUseBuildVRT = false;
815 0 : break;
816 : }
817 18 : if (GDALGetGeoTransform(hSrcDS, adfGT) != CE_None || adfGT[2] != 0 ||
818 18 : adfGT[4] != 0 || adfGT[5] > 0)
819 : {
820 0 : bCanUseBuildVRT = false;
821 0 : break;
822 : }
823 9 : const double dfMinX = adfGT[0];
824 9 : const double dfMinY = adfGT[3] + GDALGetRasterYSize(hSrcDS) * adfGT[5];
825 9 : const double dfMaxX = adfGT[0] + GDALGetRasterXSize(hSrcDS) * adfGT[1];
826 9 : const double dfMaxY = adfGT[3];
827 9 : const int nBands = GDALGetRasterCount(hSrcDS);
828 9 : if (nBands > 1 && GDALGetRasterColorInterpretation(GDALGetRasterBand(
829 : hSrcDS, nBands)) == GCI_AlphaBand)
830 : {
831 4 : bSrcHasAlpha = true;
832 : }
833 : aoExtents.emplace_back(
834 9 : std::array<double, 4>{{dfMinX, dfMinY, dfMaxX, dfMaxY}});
835 9 : const auto poSRS = GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
836 9 : if (i == 0)
837 : {
838 5 : nPrevBandCount = nBands;
839 5 : if (poSRS)
840 5 : oSRSPrev = *poSRS;
841 5 : dfLastResX = adfGT[1];
842 5 : dfLastResY = adfGT[5];
843 : }
844 : else
845 : {
846 4 : if (nPrevBandCount != nBands)
847 : {
848 0 : bCanUseBuildVRT = false;
849 0 : break;
850 : }
851 4 : if (poSRS == nullptr && !oSRSPrev.IsEmpty())
852 : {
853 0 : bCanUseBuildVRT = false;
854 0 : break;
855 : }
856 8 : if (poSRS != nullptr &&
857 4 : (oSRSPrev.IsEmpty() || !poSRS->IsSame(&oSRSPrev)))
858 : {
859 0 : bCanUseBuildVRT = false;
860 0 : break;
861 : }
862 4 : if (dfLastResX != adfGT[1] || dfLastResY != adfGT[5])
863 : {
864 0 : bCanUseBuildVRT = false;
865 0 : break;
866 : }
867 : }
868 : }
869 5 : if (bSrcHasAlpha && bCanUseBuildVRT)
870 : {
871 : // Quadratic performance loop. If that happens to be an issue,
872 : // we might need to build a quad tree
873 2 : for (size_t i = 0; i < aoExtents.size(); i++)
874 : {
875 2 : const double dfMinX = aoExtents[i][0];
876 2 : const double dfMinY = aoExtents[i][1];
877 2 : const double dfMaxX = aoExtents[i][2];
878 2 : const double dfMaxY = aoExtents[i][3];
879 2 : for (size_t j = i + 1; j < aoExtents.size(); j++)
880 : {
881 2 : const double dfOtherMinX = aoExtents[j][0];
882 2 : const double dfOtherMinY = aoExtents[j][1];
883 2 : const double dfOtherMaxX = aoExtents[j][2];
884 2 : const double dfOtherMaxY = aoExtents[j][3];
885 2 : if (dfMinX < dfOtherMaxX && dfOtherMinX < dfMaxX &&
886 2 : dfMinY < dfOtherMaxY && dfOtherMinY < dfMaxY)
887 : {
888 2 : bCanUseBuildVRT = false;
889 2 : break;
890 : }
891 : }
892 2 : if (!bCanUseBuildVRT)
893 2 : break;
894 : }
895 : }
896 10 : return bCanUseBuildVRT;
897 : }
898 :
899 : #ifdef HAVE_TIFF
900 :
901 : /************************************************************************/
902 : /* DealWithCOGOptions() */
903 : /************************************************************************/
904 :
905 8 : static bool DealWithCOGOptions(CPLStringList &aosCreateOptions, int nSrcCount,
906 : GDALDatasetH *pahSrcDS,
907 : GDALWarpAppOptions *psOptions,
908 : GDALTransformerArgUniquePtr &hUniqueTransformArg)
909 : {
910 8 : const auto SetDstSRS = [psOptions](const std::string &osTargetSRS)
911 : {
912 : const char *pszExistingDstSRS =
913 5 : psOptions->aosTransformerOptions.FetchNameValue("DST_SRS");
914 5 : if (pszExistingDstSRS)
915 : {
916 2 : OGRSpatialReference oSRS1;
917 2 : oSRS1.SetFromUserInput(pszExistingDstSRS);
918 2 : OGRSpatialReference oSRS2;
919 2 : oSRS2.SetFromUserInput(osTargetSRS.c_str());
920 2 : if (!oSRS1.IsSame(&oSRS2))
921 : {
922 2 : CPLError(CE_Failure, CPLE_AppDefined,
923 : "Target SRS implied by COG creation options is not "
924 : "the same as the one specified by -t_srs");
925 2 : return false;
926 : }
927 : }
928 : psOptions->aosTransformerOptions.SetNameValue("DST_SRS",
929 3 : osTargetSRS.c_str());
930 3 : return true;
931 8 : };
932 :
933 16 : CPLString osTargetSRS;
934 8 : if (COGGetTargetSRS(aosCreateOptions.List(), osTargetSRS))
935 : {
936 5 : if (!SetDstSRS(osTargetSRS))
937 2 : return false;
938 : }
939 :
940 6 : if (!(psOptions->dfMinX == 0 && psOptions->dfMinY == 0 &&
941 5 : psOptions->dfMaxX == 0 && psOptions->dfMaxY == 0 &&
942 5 : psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
943 5 : psOptions->nForcePixels == 0 && psOptions->nForceLines == 0))
944 : {
945 2 : if (!psOptions->bResampleAlgSpecifiedByUser && nSrcCount > 0)
946 : {
947 : try
948 : {
949 1 : GDALGetWarpResampleAlg(
950 2 : COGGetResampling(GDALDataset::FromHandle(pahSrcDS[0]),
951 1 : aosCreateOptions.List())
952 : .c_str(),
953 1 : psOptions->eResampleAlg);
954 : }
955 0 : catch (const std::invalid_argument &)
956 : {
957 : // Cannot happen actually. Coverity Scan false positive...
958 0 : CPLAssert(false);
959 : }
960 : }
961 2 : return true;
962 : }
963 :
964 8 : GDALWarpAppOptions oClonedOptions(*psOptions);
965 4 : oClonedOptions.bQuiet = true;
966 : const CPLString osTmpFilename(
967 8 : VSIMemGenerateHiddenFilename("gdalwarp_tmp.tif"));
968 8 : CPLStringList aosTmpGTiffCreateOptions;
969 4 : aosTmpGTiffCreateOptions.SetNameValue("SPARSE_OK", "YES");
970 4 : aosTmpGTiffCreateOptions.SetNameValue("TILED", "YES");
971 4 : aosTmpGTiffCreateOptions.SetNameValue("BLOCKXSIZE", "4096");
972 4 : aosTmpGTiffCreateOptions.SetNameValue("BLOCKYSIZE", "4096");
973 4 : auto hTmpDS = GDALWarpCreateOutput(
974 : nSrcCount, pahSrcDS, osTmpFilename, "GTiff",
975 : oClonedOptions.aosTransformerOptions.List(),
976 4 : aosTmpGTiffCreateOptions.List(), oClonedOptions.eOutputType,
977 : hUniqueTransformArg, false, &oClonedOptions,
978 : /* bUpdateTransformerWithDestGT = */ false);
979 4 : if (hTmpDS == nullptr)
980 : {
981 0 : return false;
982 : }
983 :
984 4 : CPLString osResampling;
985 4 : int nXSize = 0;
986 4 : int nYSize = 0;
987 4 : double dfMinX = 0;
988 4 : double dfMinY = 0;
989 4 : double dfMaxX = 0;
990 4 : double dfMaxY = 0;
991 4 : bool bRet = true;
992 4 : if (COGGetWarpingCharacteristics(GDALDataset::FromHandle(hTmpDS),
993 4 : aosCreateOptions.List(), osResampling,
994 : osTargetSRS, nXSize, nYSize, dfMinX,
995 : dfMinY, dfMaxX, dfMaxY))
996 : {
997 2 : if (!psOptions->bResampleAlgSpecifiedByUser)
998 : {
999 : try
1000 : {
1001 2 : GDALGetWarpResampleAlg(osResampling, psOptions->eResampleAlg);
1002 : }
1003 0 : catch (const std::invalid_argument &)
1004 : {
1005 : // Cannot happen actually. Coverity Scan false positive...
1006 0 : CPLAssert(false);
1007 : }
1008 : }
1009 :
1010 2 : psOptions->dfMinX = dfMinX;
1011 2 : psOptions->dfMinY = dfMinY;
1012 2 : psOptions->dfMaxX = dfMaxX;
1013 2 : psOptions->dfMaxY = dfMaxY;
1014 2 : psOptions->nForcePixels = nXSize;
1015 2 : psOptions->nForceLines = nYSize;
1016 2 : COGRemoveWarpingOptions(aosCreateOptions);
1017 : }
1018 4 : GDALClose(hTmpDS);
1019 4 : VSIUnlink(osTmpFilename);
1020 4 : return bRet;
1021 : }
1022 :
1023 : #endif
1024 :
1025 : /************************************************************************/
1026 : /* GDALWarpIndirect() */
1027 : /************************************************************************/
1028 :
1029 : static GDALDatasetH
1030 : GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS, int nSrcCount,
1031 : GDALDatasetH *pahSrcDS,
1032 : GDALTransformerArgUniquePtr hUniqueTransformArg,
1033 : GDALWarpAppOptions *psOptions, int *pbUsageError);
1034 :
1035 1028 : static int CPL_STDCALL myScaledProgress(double dfProgress, const char *,
1036 : void *pProgressData)
1037 : {
1038 1028 : return GDALScaledProgress(dfProgress, nullptr, pProgressData);
1039 : }
1040 :
1041 17 : static GDALDatasetH GDALWarpIndirect(const char *pszDest, GDALDriverH hDriver,
1042 : int nSrcCount, GDALDatasetH *pahSrcDS,
1043 : GDALWarpAppOptions *psOptions,
1044 : int *pbUsageError)
1045 : {
1046 34 : CPLStringList aosCreateOptions(psOptions->aosCreateOptions);
1047 17 : psOptions->aosCreateOptions.Clear();
1048 :
1049 : // Do not use a warped VRT input for COG output, because that would cause
1050 : // warping to be done both during overview computation and creation of
1051 : // full resolution image. Better materialize a temporary GTiff a bit later
1052 : // in that method.
1053 17 : if (nSrcCount == 1 && !EQUAL(psOptions->osFormat.c_str(), "COG"))
1054 : {
1055 6 : psOptions->osFormat = "VRT";
1056 6 : auto pfnProgress = psOptions->pfnProgress;
1057 6 : psOptions->pfnProgress = GDALDummyProgress;
1058 6 : auto pProgressData = psOptions->pProgressData;
1059 6 : psOptions->pProgressData = nullptr;
1060 :
1061 6 : auto hTmpDS = GDALWarpDirect("", nullptr, nSrcCount, pahSrcDS, nullptr,
1062 : psOptions, pbUsageError);
1063 6 : if (hTmpDS)
1064 : {
1065 6 : auto hRet = GDALCreateCopy(hDriver, pszDest, hTmpDS, FALSE,
1066 6 : aosCreateOptions.List(), pfnProgress,
1067 : pProgressData);
1068 6 : GDALClose(hTmpDS);
1069 6 : return hRet;
1070 : }
1071 0 : return nullptr;
1072 : }
1073 :
1074 : // Detect a pure mosaicing situation where a BuildVRT approach is
1075 : // sufficient.
1076 11 : GDALDatasetH hTmpDS = nullptr;
1077 11 : if (psOptions->aosTransformerOptions.empty() &&
1078 5 : psOptions->eOutputType == GDT_Unknown && psOptions->dfMinX == 0 &&
1079 5 : psOptions->dfMinY == 0 && psOptions->dfMaxX == 0 &&
1080 5 : psOptions->dfMaxY == 0 && psOptions->dfXRes == 0 &&
1081 5 : psOptions->dfYRes == 0 && psOptions->nForcePixels == 0 &&
1082 10 : psOptions->nForceLines == 0 &&
1083 21 : psOptions->osCutlineDSNameOrWKT.empty() &&
1084 5 : CanUseBuildVRT(nSrcCount, pahSrcDS))
1085 : {
1086 6 : CPLStringList aosArgv;
1087 3 : const int nBands = GDALGetRasterCount(pahSrcDS[0]);
1088 0 : if ((nBands == 1 ||
1089 0 : (nBands > 1 && GDALGetRasterColorInterpretation(GDALGetRasterBand(
1090 6 : pahSrcDS[0], nBands)) != GCI_AlphaBand)) &&
1091 3 : (psOptions->bEnableDstAlpha
1092 : #ifdef HAVE_TIFF
1093 3 : || (EQUAL(psOptions->osFormat.c_str(), "COG") &&
1094 3 : COGHasWarpingOptions(aosCreateOptions.List()) &&
1095 2 : CPLTestBool(
1096 : aosCreateOptions.FetchNameValueDef("ADD_ALPHA", "YES")))
1097 : #endif
1098 : ))
1099 : {
1100 2 : aosArgv.AddString("-addalpha");
1101 : }
1102 : auto psBuildVRTOptions =
1103 3 : GDALBuildVRTOptionsNew(aosArgv.List(), nullptr);
1104 3 : hTmpDS = GDALBuildVRT("", nSrcCount, pahSrcDS, nullptr,
1105 : psBuildVRTOptions, nullptr);
1106 3 : GDALBuildVRTOptionsFree(psBuildVRTOptions);
1107 : }
1108 11 : auto pfnProgress = psOptions->pfnProgress;
1109 11 : auto pProgressData = psOptions->pProgressData;
1110 22 : CPLString osTmpFilename;
1111 11 : double dfStartPctCreateCopy = 0.0;
1112 11 : if (hTmpDS == nullptr)
1113 : {
1114 0 : GDALTransformerArgUniquePtr hUniqueTransformArg;
1115 : #ifdef HAVE_TIFF
1116 : // Special processing for COG output. As some of its options do
1117 : // on-the-fly reprojection, take them into account now, and remove them
1118 : // from the COG creation stage.
1119 16 : if (EQUAL(psOptions->osFormat.c_str(), "COG") &&
1120 8 : !DealWithCOGOptions(aosCreateOptions, nSrcCount, pahSrcDS,
1121 : psOptions, hUniqueTransformArg))
1122 : {
1123 2 : return nullptr;
1124 : }
1125 : #endif
1126 :
1127 : // Materialize a temporary GeoTIFF with the result of the warp
1128 6 : psOptions->osFormat = "GTiff";
1129 6 : psOptions->aosCreateOptions.AddString("SPARSE_OK=YES");
1130 6 : psOptions->aosCreateOptions.AddString("COMPRESS=LZW");
1131 6 : psOptions->aosCreateOptions.AddString("TILED=YES");
1132 6 : psOptions->aosCreateOptions.AddString("BIGTIFF=YES");
1133 6 : psOptions->pfnProgress = myScaledProgress;
1134 6 : dfStartPctCreateCopy = 2. / 3;
1135 6 : psOptions->pProgressData = GDALCreateScaledProgress(
1136 : 0, dfStartPctCreateCopy, pfnProgress, pProgressData);
1137 6 : psOptions->bDeleteOutputFileOnceCreated = true;
1138 6 : osTmpFilename = CPLGenerateTempFilenameSafe(CPLGetFilename(pszDest));
1139 6 : hTmpDS = GDALWarpDirect(osTmpFilename, nullptr, nSrcCount, pahSrcDS,
1140 6 : std::move(hUniqueTransformArg), psOptions,
1141 : pbUsageError);
1142 6 : GDALDestroyScaledProgress(psOptions->pProgressData);
1143 6 : psOptions->pfnProgress = nullptr;
1144 6 : psOptions->pProgressData = nullptr;
1145 : }
1146 9 : if (hTmpDS)
1147 : {
1148 9 : auto pScaledProgressData = GDALCreateScaledProgress(
1149 : dfStartPctCreateCopy, 1.0, pfnProgress, pProgressData);
1150 9 : auto hRet = GDALCreateCopy(hDriver, pszDest, hTmpDS, FALSE,
1151 9 : aosCreateOptions.List(), myScaledProgress,
1152 : pScaledProgressData);
1153 9 : GDALDestroyScaledProgress(pScaledProgressData);
1154 9 : GDALClose(hTmpDS);
1155 : VSIStatBufL sStat;
1156 15 : if (!osTmpFilename.empty() &&
1157 6 : VSIStatL(osTmpFilename.c_str(), &sStat) == 0)
1158 : {
1159 6 : GDALDeleteDataset(GDALGetDriverByName("GTiff"), osTmpFilename);
1160 : }
1161 9 : return hRet;
1162 : }
1163 0 : return nullptr;
1164 : }
1165 :
1166 : /************************************************************************/
1167 : /* GDALWarp() */
1168 : /************************************************************************/
1169 :
1170 : /**
1171 : * Image reprojection and warping function.
1172 : *
1173 : * This is the equivalent of the <a href="/programs/gdalwarp.html">gdalwarp</a>
1174 : * utility.
1175 : *
1176 : * GDALWarpAppOptions* must be allocated and freed with GDALWarpAppOptionsNew()
1177 : * and GDALWarpAppOptionsFree() respectively.
1178 : * pszDest and hDstDS cannot be used at the same time.
1179 : *
1180 : * @param pszDest the destination dataset path or NULL.
1181 : * @param hDstDS the destination dataset or NULL.
1182 : * @param nSrcCount the number of input datasets.
1183 : * @param pahSrcDS the list of input datasets. For practical purposes, the type
1184 : * of this argument should be considered as "const GDALDatasetH* const*", that
1185 : * is neither the array nor its values are mutated by this function.
1186 : * @param psOptionsIn the options struct returned by GDALWarpAppOptionsNew() or
1187 : * NULL.
1188 : * @param pbUsageError pointer to a integer output variable to store if any
1189 : * usage error has occurred, or NULL.
1190 : * @return the output dataset (new dataset that must be closed using
1191 : * GDALClose(), or hDstDS if not NULL) or NULL in case of error. If the output
1192 : * format is a VRT dataset, then the returned VRT dataset has a reference to
1193 : * pahSrcDS[0]. Hence pahSrcDS[0] should be closed after the returned dataset
1194 : * if using GDALClose().
1195 : * A safer alternative is to use GDALReleaseDataset() instead of using
1196 : * GDALClose(), in which case you can close datasets in any order.
1197 : *
1198 : * @since GDAL 2.1
1199 : */
1200 :
1201 1174 : GDALDatasetH GDALWarp(const char *pszDest, GDALDatasetH hDstDS, int nSrcCount,
1202 : GDALDatasetH *pahSrcDS,
1203 : const GDALWarpAppOptions *psOptionsIn, int *pbUsageError)
1204 : {
1205 1174 : CPLErrorReset();
1206 :
1207 2371 : for (int i = 0; i < nSrcCount; i++)
1208 : {
1209 1197 : if (!pahSrcDS[i])
1210 0 : return nullptr;
1211 : }
1212 :
1213 2348 : GDALWarpAppOptions oOptionsTmp;
1214 1174 : if (psOptionsIn)
1215 1172 : oOptionsTmp = *psOptionsIn;
1216 1174 : GDALWarpAppOptions *psOptions = &oOptionsTmp;
1217 :
1218 1174 : if (hDstDS == nullptr)
1219 : {
1220 1081 : if (psOptions->osFormat.empty())
1221 : {
1222 417 : psOptions->osFormat = GetOutputDriverForRaster(pszDest);
1223 417 : if (psOptions->osFormat.empty())
1224 : {
1225 0 : return nullptr;
1226 : }
1227 : }
1228 :
1229 1081 : auto hDriver = GDALGetDriverByName(psOptions->osFormat.c_str());
1230 2162 : if (hDriver != nullptr &&
1231 1081 : (EQUAL(psOptions->osFormat.c_str(), "COG") ||
1232 1070 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) ==
1233 6 : nullptr &&
1234 6 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) !=
1235 : nullptr)))
1236 : {
1237 17 : auto ret = GDALWarpIndirect(pszDest, hDriver, nSrcCount, pahSrcDS,
1238 : psOptions, pbUsageError);
1239 17 : return ret;
1240 : }
1241 : }
1242 :
1243 1157 : auto ret = GDALWarpDirect(pszDest, hDstDS, nSrcCount, pahSrcDS, nullptr,
1244 : psOptions, pbUsageError);
1245 :
1246 1157 : return ret;
1247 : }
1248 :
1249 : /************************************************************************/
1250 : /* UseTEAndTSAndTRConsistently() */
1251 : /************************************************************************/
1252 :
1253 1081 : static bool UseTEAndTSAndTRConsistently(const GDALWarpAppOptions *psOptions)
1254 : {
1255 : // We normally don't allow -te, -ts and -tr together, unless they are all
1256 : // consistent. The interest of this is to use the -tr values to produce
1257 : // exact pixel size, rather than inferring it from -te and -ts
1258 :
1259 : // Constant and logic to be kept in sync with cogdriver.cpp
1260 1081 : constexpr double RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP = 1e-8;
1261 199 : return psOptions->nForcePixels != 0 && psOptions->nForceLines != 0 &&
1262 194 : psOptions->dfXRes != 0 && psOptions->dfYRes != 0 &&
1263 52 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
1264 0 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0) &&
1265 52 : fabs((psOptions->dfMaxX - psOptions->dfMinX) / psOptions->dfXRes -
1266 52 : psOptions->nForcePixels) <=
1267 1280 : RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP &&
1268 52 : fabs((psOptions->dfMaxY - psOptions->dfMinY) / psOptions->dfYRes -
1269 52 : psOptions->nForceLines) <=
1270 1081 : RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP;
1271 : }
1272 :
1273 : /************************************************************************/
1274 : /* CheckOptions() */
1275 : /************************************************************************/
1276 :
1277 1169 : static bool CheckOptions(const char *pszDest, GDALDatasetH hDstDS,
1278 : int nSrcCount, GDALDatasetH *pahSrcDS,
1279 : GDALWarpAppOptions *psOptions, bool &bVRT,
1280 : int *pbUsageError)
1281 : {
1282 :
1283 1169 : if (hDstDS)
1284 : {
1285 93 : if (psOptions->bCreateOutput == true)
1286 : {
1287 0 : CPLError(CE_Warning, CPLE_AppDefined,
1288 : "All options related to creation ignored in update mode");
1289 0 : psOptions->bCreateOutput = false;
1290 : }
1291 : }
1292 :
1293 2524 : if ((psOptions->osFormat.empty() &&
1294 2338 : EQUAL(CPLGetExtensionSafe(pszDest).c_str(), "VRT")) ||
1295 1169 : (EQUAL(psOptions->osFormat.c_str(), "VRT")))
1296 : {
1297 181 : if (hDstDS != nullptr)
1298 : {
1299 0 : CPLError(CE_Warning, CPLE_NotSupported,
1300 : "VRT output not compatible with existing dataset.");
1301 0 : return false;
1302 : }
1303 :
1304 181 : bVRT = true;
1305 :
1306 181 : if (nSrcCount > 1)
1307 : {
1308 0 : CPLError(CE_Warning, CPLE_AppDefined,
1309 : "gdalwarp -of VRT just takes into account "
1310 : "the first source dataset.\nIf all source datasets "
1311 : "are in the same projection, try making a mosaic of\n"
1312 : "them with gdalbuildvrt, and use the resulting "
1313 : "VRT file as the input of\ngdalwarp -of VRT.");
1314 : }
1315 : }
1316 :
1317 : /* -------------------------------------------------------------------- */
1318 : /* Check that incompatible options are not used */
1319 : /* -------------------------------------------------------------------- */
1320 :
1321 996 : if ((psOptions->nForcePixels != 0 || psOptions->nForceLines != 0) &&
1322 2191 : (psOptions->dfXRes != 0 && psOptions->dfYRes != 0) &&
1323 26 : !UseTEAndTSAndTRConsistently(psOptions))
1324 : {
1325 0 : CPLError(CE_Failure, CPLE_IllegalArg,
1326 : "-tr and -ts options cannot be used at the same time.");
1327 0 : if (pbUsageError)
1328 0 : *pbUsageError = TRUE;
1329 0 : return false;
1330 : }
1331 :
1332 1169 : if (psOptions->bTargetAlignedPixels && psOptions->dfXRes == 0 &&
1333 1 : psOptions->dfYRes == 0)
1334 : {
1335 1 : CPLError(CE_Failure, CPLE_IllegalArg,
1336 : "-tap option cannot be used without using -tr.");
1337 1 : if (pbUsageError)
1338 1 : *pbUsageError = TRUE;
1339 1 : return false;
1340 : }
1341 :
1342 1168 : if (!psOptions->bQuiet &&
1343 81 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
1344 64 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0))
1345 : {
1346 17 : if (psOptions->dfMinX >= psOptions->dfMaxX)
1347 0 : CPLError(CE_Warning, CPLE_AppDefined,
1348 : "-te values have minx >= maxx. This will result in a "
1349 : "horizontally flipped image.");
1350 17 : if (psOptions->dfMinY >= psOptions->dfMaxY)
1351 0 : CPLError(CE_Warning, CPLE_AppDefined,
1352 : "-te values have miny >= maxy. This will result in a "
1353 : "vertically flipped image.");
1354 : }
1355 :
1356 1168 : if (psOptions->dfErrorThreshold < 0)
1357 : {
1358 : // By default, use approximate transformer unless RPC_DEM is specified
1359 1155 : if (psOptions->aosTransformerOptions.FetchNameValue("RPC_DEM") !=
1360 : nullptr)
1361 4 : psOptions->dfErrorThreshold = 0.0;
1362 : else
1363 1151 : psOptions->dfErrorThreshold = 0.125;
1364 : }
1365 :
1366 : /* -------------------------------------------------------------------- */
1367 : /* -te_srs option */
1368 : /* -------------------------------------------------------------------- */
1369 1168 : if (!psOptions->osTE_SRS.empty())
1370 : {
1371 6 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
1372 0 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
1373 : {
1374 0 : CPLError(CE_Warning, CPLE_AppDefined,
1375 : "-te_srs ignored since -te is not specified.");
1376 : }
1377 : else
1378 : {
1379 6 : OGRSpatialReference oSRSIn;
1380 6 : oSRSIn.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1381 6 : oSRSIn.SetFromUserInput(psOptions->osTE_SRS.c_str());
1382 6 : OGRSpatialReference oSRSDS;
1383 6 : oSRSDS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1384 6 : bool bOK = false;
1385 6 : if (psOptions->aosTransformerOptions.FetchNameValue("DST_SRS") !=
1386 : nullptr)
1387 : {
1388 3 : oSRSDS.SetFromUserInput(
1389 : psOptions->aosTransformerOptions.FetchNameValue("DST_SRS"));
1390 3 : bOK = true;
1391 : }
1392 3 : else if (psOptions->aosTransformerOptions.FetchNameValue(
1393 3 : "SRC_SRS") != nullptr)
1394 : {
1395 0 : oSRSDS.SetFromUserInput(
1396 : psOptions->aosTransformerOptions.FetchNameValue("SRC_SRS"));
1397 0 : bOK = true;
1398 : }
1399 3 : else if (nSrcCount)
1400 : {
1401 : const auto poSrcSRS =
1402 3 : GDALDataset::FromHandle(pahSrcDS[0])->GetSpatialRef();
1403 3 : if (poSrcSRS)
1404 : {
1405 3 : oSRSDS = *poSrcSRS;
1406 3 : bOK = true;
1407 : }
1408 : }
1409 6 : if (!bOK)
1410 : {
1411 0 : CPLError(CE_Failure, CPLE_AppDefined,
1412 : "-te_srs ignored since none of -t_srs, -s_srs is "
1413 : "specified or the input dataset has no projection.");
1414 0 : return false;
1415 : }
1416 6 : if (!oSRSIn.IsSame(&oSRSDS))
1417 : {
1418 5 : double dfWestLongitudeDeg = 0.0;
1419 5 : double dfSouthLatitudeDeg = 0.0;
1420 5 : double dfEastLongitudeDeg = 0.0;
1421 5 : double dfNorthLatitudeDeg = 0.0;
1422 :
1423 5 : OGRCoordinateTransformationOptions options;
1424 5 : if (GDALComputeAreaOfInterest(
1425 : &oSRSIn, psOptions->dfMinX, psOptions->dfMinY,
1426 : psOptions->dfMaxX, psOptions->dfMaxY,
1427 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
1428 : dfEastLongitudeDeg, dfNorthLatitudeDeg))
1429 : {
1430 5 : options.SetAreaOfInterest(
1431 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
1432 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
1433 : }
1434 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
1435 : OGRCreateCoordinateTransformation(&oSRSIn, &oSRSDS,
1436 5 : options));
1437 5 : constexpr int DENSIFY_PTS = 21;
1438 10 : if (!(poCT && poCT->TransformBounds(
1439 : psOptions->dfMinX, psOptions->dfMinY,
1440 : psOptions->dfMaxX, psOptions->dfMaxY,
1441 : &(psOptions->dfMinX), &(psOptions->dfMinY),
1442 : &(psOptions->dfMaxX), &(psOptions->dfMaxY),
1443 10 : DENSIFY_PTS)))
1444 : {
1445 0 : CPLError(CE_Failure, CPLE_AppDefined,
1446 : "-te_srs ignored since coordinate transformation "
1447 : "failed.");
1448 0 : return false;
1449 : }
1450 : }
1451 : }
1452 : }
1453 1168 : return true;
1454 : }
1455 :
1456 : /************************************************************************/
1457 : /* ProcessCutlineOptions() */
1458 : /************************************************************************/
1459 :
1460 1168 : static bool ProcessCutlineOptions(int nSrcCount, GDALDatasetH *pahSrcDS,
1461 : GDALWarpAppOptions *psOptions,
1462 : std::unique_ptr<OGRGeometry> &poCutline)
1463 : {
1464 1168 : if (!psOptions->osCutlineDSNameOrWKT.empty())
1465 : {
1466 : CPLErr eError;
1467 58 : OGRGeometryH hCutline = nullptr;
1468 116 : eError = LoadCutline(psOptions->osCutlineDSNameOrWKT,
1469 58 : psOptions->osCutlineSRS, psOptions->osCLayer,
1470 58 : psOptions->osCWHERE, psOptions->osCSQL, &hCutline);
1471 58 : poCutline.reset(OGRGeometry::FromHandle(hCutline));
1472 58 : if (eError == CE_Failure)
1473 : {
1474 6 : return false;
1475 : }
1476 : }
1477 :
1478 1162 : if (psOptions->bCropToCutline && poCutline)
1479 : {
1480 : CPLErr eError;
1481 18 : eError = CropToCutline(poCutline.get(),
1482 18 : psOptions->aosTransformerOptions.List(),
1483 18 : psOptions->aosWarpOptions.List(), nSrcCount,
1484 18 : pahSrcDS, psOptions->dfMinX, psOptions->dfMinY,
1485 18 : psOptions->dfMaxX, psOptions->dfMaxY, psOptions);
1486 18 : if (eError == CE_Failure)
1487 : {
1488 5 : return false;
1489 : }
1490 : }
1491 :
1492 : const char *pszWarpThreads =
1493 1157 : psOptions->aosWarpOptions.FetchNameValue("NUM_THREADS");
1494 1157 : if (pszWarpThreads != nullptr)
1495 : {
1496 : /* Used by TPS transformer to parallelize direct and inverse matrix
1497 : * computation */
1498 : psOptions->aosTransformerOptions.SetNameValue("NUM_THREADS",
1499 28 : pszWarpThreads);
1500 : }
1501 :
1502 1157 : return true;
1503 : }
1504 :
1505 : /************************************************************************/
1506 : /* CreateOutput() */
1507 : /************************************************************************/
1508 :
1509 : static GDALDatasetH
1510 1064 : CreateOutput(const char *pszDest, int nSrcCount, GDALDatasetH *pahSrcDS,
1511 : GDALWarpAppOptions *psOptions, const bool bInitDestSetByUser,
1512 : GDALTransformerArgUniquePtr &hUniqueTransformArg)
1513 : {
1514 1064 : if (nSrcCount == 1 && !psOptions->bDisableSrcAlpha)
1515 : {
1516 2066 : if (GDALGetRasterCount(pahSrcDS[0]) > 0 &&
1517 1033 : GDALGetRasterColorInterpretation(GDALGetRasterBand(
1518 : pahSrcDS[0], GDALGetRasterCount(pahSrcDS[0]))) == GCI_AlphaBand)
1519 : {
1520 39 : psOptions->bEnableSrcAlpha = true;
1521 39 : psOptions->bEnableDstAlpha = true;
1522 39 : if (!psOptions->bQuiet)
1523 0 : printf("Using band %d of source image as alpha.\n",
1524 : GDALGetRasterCount(pahSrcDS[0]));
1525 : }
1526 : }
1527 :
1528 1064 : auto hDstDS = GDALWarpCreateOutput(
1529 : nSrcCount, pahSrcDS, pszDest, psOptions->osFormat.c_str(),
1530 : psOptions->aosTransformerOptions.List(),
1531 1064 : psOptions->aosCreateOptions.List(), psOptions->eOutputType,
1532 1064 : hUniqueTransformArg, psOptions->bSetColorInterpretation, psOptions,
1533 : /* bUpdateTransformerWithDestGT = */ true);
1534 1064 : if (hDstDS == nullptr)
1535 : {
1536 19 : return nullptr;
1537 : }
1538 1045 : psOptions->bCreateOutput = true;
1539 :
1540 1045 : if (!bInitDestSetByUser)
1541 : {
1542 1025 : if (psOptions->osDstNodata.empty())
1543 : {
1544 893 : psOptions->aosWarpOptions.SetNameValue("INIT_DEST", "0");
1545 : }
1546 : else
1547 : {
1548 132 : psOptions->aosWarpOptions.SetNameValue("INIT_DEST", "NO_DATA");
1549 : }
1550 : }
1551 :
1552 1045 : return hDstDS;
1553 : }
1554 :
1555 : /************************************************************************/
1556 : /* EditISIS3ForMetadataChanges() */
1557 : /************************************************************************/
1558 :
1559 : static std::string
1560 2 : EditISIS3ForMetadataChanges(const char *pszJSON,
1561 : const GDALWarpAppOptions *psOptions)
1562 : {
1563 4 : CPLJSONDocument oJSONDocument;
1564 2 : if (!oJSONDocument.LoadMemory(pszJSON))
1565 : {
1566 0 : return std::string();
1567 : }
1568 :
1569 4 : auto oRoot = oJSONDocument.GetRoot();
1570 2 : if (!oRoot.IsValid())
1571 : {
1572 0 : return std::string();
1573 : }
1574 :
1575 6 : auto oGDALHistory = oRoot.GetObj("GDALHistory");
1576 2 : if (!oGDALHistory.IsValid())
1577 : {
1578 2 : oGDALHistory = CPLJSONObject();
1579 2 : oRoot.Add("GDALHistory", oGDALHistory);
1580 : }
1581 2 : oGDALHistory["_type"] = "object";
1582 :
1583 2 : char szFullFilename[2048] = {0};
1584 2 : if (!CPLGetExecPath(szFullFilename, sizeof(szFullFilename) - 1))
1585 0 : strcpy(szFullFilename, "unknown_program");
1586 4 : const CPLString osProgram(CPLGetBasenameSafe(szFullFilename));
1587 4 : const CPLString osPath(CPLGetPathSafe(szFullFilename));
1588 :
1589 2 : oGDALHistory["GdalVersion"] = GDALVersionInfo("RELEASE_NAME");
1590 2 : oGDALHistory["Program"] = osProgram;
1591 2 : if (osPath != ".")
1592 2 : oGDALHistory["ProgramPath"] = osPath;
1593 :
1594 4 : std::string osArgs;
1595 6 : for (const char *pszArg : psOptions->aosArgs)
1596 : {
1597 4 : if (!osArgs.empty())
1598 2 : osArgs += ' ';
1599 4 : osArgs += pszArg;
1600 : }
1601 2 : oGDALHistory["ProgramArguments"] = osArgs;
1602 :
1603 2 : oGDALHistory.Add(
1604 : "Comment",
1605 : "Part of that metadata might be invalid due to a reprojection "
1606 : "operation having been performed by GDAL tools");
1607 :
1608 2 : return oRoot.Format(CPLJSONObject::PrettyFormat::Pretty);
1609 : }
1610 :
1611 : /************************************************************************/
1612 : /* ProcessMetadata() */
1613 : /************************************************************************/
1614 :
1615 1159 : static void ProcessMetadata(int iSrc, GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
1616 : const GDALWarpAppOptions *psOptions,
1617 : const bool bEnableDstAlpha,
1618 : bool bVerticalShiftApplied)
1619 : {
1620 1159 : if (psOptions->bCopyMetadata)
1621 : {
1622 1159 : const char *pszSrcInfo = nullptr;
1623 1159 : GDALRasterBandH hSrcBand = nullptr;
1624 1159 : GDALRasterBandH hDstBand = nullptr;
1625 :
1626 : /* copy metadata from first dataset */
1627 1159 : if (iSrc == 0)
1628 : {
1629 1135 : CPLDebug(
1630 : "WARP",
1631 : "Copying metadata from first source to destination dataset");
1632 : /* copy dataset-level metadata */
1633 1135 : CSLConstList papszMetadata = GDALGetMetadata(hSrcDS, nullptr);
1634 :
1635 1135 : char **papszMetadataNew = nullptr;
1636 2947 : for (int i = 0;
1637 2947 : papszMetadata != nullptr && papszMetadata[i] != nullptr; i++)
1638 : {
1639 : // Do not preserve NODATA_VALUES when the output includes an
1640 : // alpha band
1641 1812 : if (bEnableDstAlpha &&
1642 198 : STARTS_WITH_CI(papszMetadata[i], "NODATA_VALUES="))
1643 : {
1644 1 : continue;
1645 : }
1646 : // Do not preserve the CACHE_PATH from the WMS driver
1647 1811 : if (STARTS_WITH_CI(papszMetadata[i], "CACHE_PATH="))
1648 : {
1649 0 : continue;
1650 : }
1651 :
1652 : papszMetadataNew =
1653 1811 : CSLAddString(papszMetadataNew, papszMetadata[i]);
1654 : }
1655 :
1656 1135 : if (CSLCount(papszMetadataNew) > 0)
1657 : {
1658 759 : if (GDALSetMetadata(hDstDS, papszMetadataNew, nullptr) !=
1659 : CE_None)
1660 0 : CPLError(CE_Warning, CPLE_AppDefined,
1661 : "error copying metadata to destination dataset.");
1662 : }
1663 :
1664 1135 : CSLDestroy(papszMetadataNew);
1665 :
1666 : /* ISIS3 -> ISIS3 special case */
1667 1135 : if (EQUAL(psOptions->osFormat.c_str(), "ISIS3") ||
1668 1134 : EQUAL(psOptions->osFormat.c_str(), "PDS4") ||
1669 2953 : EQUAL(psOptions->osFormat.c_str(), "GTIFF") ||
1670 684 : EQUAL(psOptions->osFormat.c_str(), "COG"))
1671 : {
1672 : CSLConstList papszMD_ISIS3 =
1673 451 : GDALGetMetadata(hSrcDS, "json:ISIS3");
1674 451 : if (papszMD_ISIS3 != nullptr && papszMD_ISIS3[0])
1675 : {
1676 4 : std::string osJSON = papszMD_ISIS3[0];
1677 : osJSON =
1678 2 : EditISIS3ForMetadataChanges(osJSON.c_str(), psOptions);
1679 2 : if (!osJSON.empty())
1680 : {
1681 2 : char *apszMD[] = {osJSON.data(), nullptr};
1682 2 : GDALSetMetadata(hDstDS, apszMD, "json:ISIS3");
1683 : }
1684 : }
1685 : }
1686 684 : else if (EQUAL(psOptions->osFormat.c_str(), "PDS4"))
1687 : {
1688 0 : CSLConstList papszMD_PDS4 = GDALGetMetadata(hSrcDS, "xml:PDS4");
1689 0 : if (papszMD_PDS4 != nullptr)
1690 0 : GDALSetMetadata(hDstDS, papszMD_PDS4, "xml:PDS4");
1691 : }
1692 684 : else if (EQUAL(psOptions->osFormat.c_str(), "VICAR"))
1693 : {
1694 : CSLConstList papszMD_VICAR =
1695 0 : GDALGetMetadata(hSrcDS, "json:VICAR");
1696 0 : if (papszMD_VICAR != nullptr)
1697 0 : GDALSetMetadata(hDstDS, papszMD_VICAR, "json:VICAR");
1698 : }
1699 :
1700 : /* copy band-level metadata and other info */
1701 1135 : if (GDALGetRasterCount(hSrcDS) == GDALGetRasterCount(hDstDS))
1702 : {
1703 2458 : for (int iBand = 0; iBand < GDALGetRasterCount(hSrcDS); iBand++)
1704 : {
1705 1406 : hSrcBand = GDALGetRasterBand(hSrcDS, iBand + 1);
1706 1406 : hDstBand = GDALGetRasterBand(hDstDS, iBand + 1);
1707 : /* copy metadata, except stats (#5319) */
1708 1406 : papszMetadata = GDALGetMetadata(hSrcBand, nullptr);
1709 1406 : if (CSLCount(papszMetadata) > 0)
1710 : {
1711 : // GDALSetMetadata( hDstBand, papszMetadata, NULL );
1712 22 : papszMetadataNew = nullptr;
1713 120 : for (int i = 0; papszMetadata != nullptr &&
1714 120 : papszMetadata[i] != nullptr;
1715 : i++)
1716 : {
1717 98 : if (!STARTS_WITH(papszMetadata[i], "STATISTICS_"))
1718 83 : papszMetadataNew = CSLAddString(
1719 83 : papszMetadataNew, papszMetadata[i]);
1720 : }
1721 22 : GDALSetMetadata(hDstBand, papszMetadataNew, nullptr);
1722 22 : CSLDestroy(papszMetadataNew);
1723 : }
1724 : /* copy other info (Description, Unit Type) - what else? */
1725 1406 : if (psOptions->bCopyBandInfo && !bVerticalShiftApplied)
1726 : {
1727 1386 : pszSrcInfo = GDALGetDescription(hSrcBand);
1728 1386 : if (pszSrcInfo != nullptr && strlen(pszSrcInfo) > 0)
1729 3 : GDALSetDescription(hDstBand, pszSrcInfo);
1730 1386 : pszSrcInfo = GDALGetRasterUnitType(hSrcBand);
1731 1386 : if (pszSrcInfo != nullptr && strlen(pszSrcInfo) > 0)
1732 34 : GDALSetRasterUnitType(hDstBand, pszSrcInfo);
1733 : }
1734 : }
1735 : }
1736 : }
1737 : /* remove metadata that conflicts between datasets */
1738 : else
1739 : {
1740 24 : CPLDebug("WARP",
1741 : "Removing conflicting metadata from destination dataset "
1742 : "(source #%d)",
1743 : iSrc);
1744 : /* remove conflicting dataset-level metadata */
1745 24 : RemoveConflictingMetadata(hDstDS, GDALGetMetadata(hSrcDS, nullptr),
1746 : psOptions->osMDConflictValue.c_str());
1747 :
1748 : /* remove conflicting copy band-level metadata and other info */
1749 24 : if (GDALGetRasterCount(hSrcDS) == GDALGetRasterCount(hDstDS))
1750 : {
1751 44 : for (int iBand = 0; iBand < GDALGetRasterCount(hSrcDS); iBand++)
1752 : {
1753 23 : hSrcBand = GDALGetRasterBand(hSrcDS, iBand + 1);
1754 23 : hDstBand = GDALGetRasterBand(hDstDS, iBand + 1);
1755 : /* remove conflicting metadata */
1756 23 : RemoveConflictingMetadata(
1757 : hDstBand, GDALGetMetadata(hSrcBand, nullptr),
1758 : psOptions->osMDConflictValue.c_str());
1759 : /* remove conflicting info */
1760 23 : if (psOptions->bCopyBandInfo)
1761 : {
1762 23 : pszSrcInfo = GDALGetDescription(hSrcBand);
1763 23 : const char *pszDstInfo = GDALGetDescription(hDstBand);
1764 23 : if (!(pszSrcInfo != nullptr && strlen(pszSrcInfo) > 0 &&
1765 0 : pszDstInfo != nullptr && strlen(pszDstInfo) > 0 &&
1766 0 : EQUAL(pszSrcInfo, pszDstInfo)))
1767 23 : GDALSetDescription(hDstBand, "");
1768 23 : pszSrcInfo = GDALGetRasterUnitType(hSrcBand);
1769 23 : pszDstInfo = GDALGetRasterUnitType(hDstBand);
1770 23 : if (!(pszSrcInfo != nullptr && strlen(pszSrcInfo) > 0 &&
1771 0 : pszDstInfo != nullptr && strlen(pszDstInfo) > 0 &&
1772 0 : EQUAL(pszSrcInfo, pszDstInfo)))
1773 23 : GDALSetRasterUnitType(hDstBand, "");
1774 : }
1775 : }
1776 : }
1777 : }
1778 : }
1779 1159 : }
1780 :
1781 : /************************************************************************/
1782 : /* SetupNoData() */
1783 : /************************************************************************/
1784 :
1785 1156 : static CPLErr SetupNoData(const char *pszDest, int iSrc, GDALDatasetH hSrcDS,
1786 : GDALDatasetH hWrkSrcDS, GDALDatasetH hDstDS,
1787 : GDALWarpOptions *psWO, GDALWarpAppOptions *psOptions,
1788 : const bool bEnableDstAlpha,
1789 : const bool bInitDestSetByUser)
1790 : {
1791 1206 : if (!psOptions->osSrcNodata.empty() &&
1792 50 : !EQUAL(psOptions->osSrcNodata.c_str(), "none"))
1793 : {
1794 : CPLStringList aosTokens(
1795 46 : CSLTokenizeString(psOptions->osSrcNodata.c_str()));
1796 46 : const int nTokenCount = aosTokens.Count();
1797 :
1798 46 : psWO->padfSrcNoDataReal =
1799 46 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
1800 46 : psWO->padfSrcNoDataImag = nullptr;
1801 :
1802 123 : for (int i = 0; i < psWO->nBandCount; i++)
1803 : {
1804 78 : if (i < nTokenCount)
1805 : {
1806 : double dfNoDataReal;
1807 : double dfNoDataImag;
1808 :
1809 48 : if (CPLStringToComplex(aosTokens[i], &dfNoDataReal,
1810 48 : &dfNoDataImag) != CE_None)
1811 : {
1812 1 : CPLError(CE_Failure, CPLE_AppDefined,
1813 : "Error parsing srcnodata for band %d", i + 1);
1814 1 : return CE_Failure;
1815 : }
1816 :
1817 94 : psWO->padfSrcNoDataReal[i] =
1818 47 : GDALAdjustNoDataCloseToFloatMax(dfNoDataReal);
1819 :
1820 47 : if (strchr(aosTokens[i], 'i') != nullptr)
1821 : {
1822 1 : if (psWO->padfSrcNoDataImag == nullptr)
1823 : {
1824 1 : psWO->padfSrcNoDataImag = static_cast<double *>(
1825 1 : CPLCalloc(psWO->nBandCount, sizeof(double)));
1826 : }
1827 2 : psWO->padfSrcNoDataImag[i] =
1828 1 : GDALAdjustNoDataCloseToFloatMax(dfNoDataImag);
1829 : }
1830 : }
1831 : else
1832 : {
1833 30 : psWO->padfSrcNoDataReal[i] = psWO->padfSrcNoDataReal[i - 1];
1834 30 : if (psWO->padfSrcNoDataImag != nullptr)
1835 : {
1836 0 : psWO->padfSrcNoDataImag[i] = psWO->padfSrcNoDataImag[i - 1];
1837 : }
1838 : }
1839 : }
1840 :
1841 71 : if (psWO->nBandCount > 1 &&
1842 26 : CSLFetchNameValue(psWO->papszWarpOptions, "UNIFIED_SRC_NODATA") ==
1843 : nullptr)
1844 : {
1845 6 : CPLDebug("WARP", "Set UNIFIED_SRC_NODATA=YES");
1846 6 : psWO->papszWarpOptions = CSLSetNameValue(
1847 : psWO->papszWarpOptions, "UNIFIED_SRC_NODATA", "YES");
1848 : }
1849 : }
1850 :
1851 : /* -------------------------------------------------------------------- */
1852 : /* If -srcnodata was not specified, but the data has nodata */
1853 : /* values, use them. */
1854 : /* -------------------------------------------------------------------- */
1855 1155 : if (psOptions->osSrcNodata.empty())
1856 : {
1857 1106 : int bHaveNodata = FALSE;
1858 1106 : double dfReal = 0.0;
1859 :
1860 2523 : for (int i = 0; !bHaveNodata && i < psWO->nBandCount; i++)
1861 : {
1862 : GDALRasterBandH hBand =
1863 1417 : GDALGetRasterBand(hWrkSrcDS, psWO->panSrcBands[i]);
1864 1417 : dfReal = GDALGetRasterNoDataValue(hBand, &bHaveNodata);
1865 : }
1866 :
1867 1106 : if (bHaveNodata)
1868 : {
1869 96 : if (!psOptions->bQuiet)
1870 : {
1871 10 : if (std::isnan(dfReal))
1872 0 : printf("Using internal nodata values (e.g. nan) for image "
1873 : "%s.\n",
1874 : GDALGetDescription(hSrcDS));
1875 : else
1876 10 : printf("Using internal nodata values (e.g. %g) for image "
1877 : "%s.\n",
1878 : dfReal, GDALGetDescription(hSrcDS));
1879 : }
1880 96 : psWO->padfSrcNoDataReal = static_cast<double *>(
1881 96 : CPLMalloc(psWO->nBandCount * sizeof(double)));
1882 :
1883 229 : for (int i = 0; i < psWO->nBandCount; i++)
1884 : {
1885 : GDALRasterBandH hBand =
1886 133 : GDALGetRasterBand(hWrkSrcDS, psWO->panSrcBands[i]);
1887 :
1888 133 : dfReal = GDALGetRasterNoDataValue(hBand, &bHaveNodata);
1889 :
1890 133 : if (bHaveNodata)
1891 : {
1892 133 : psWO->padfSrcNoDataReal[i] = dfReal;
1893 : }
1894 : else
1895 : {
1896 0 : psWO->padfSrcNoDataReal[i] = -123456.789;
1897 : }
1898 : }
1899 : }
1900 : }
1901 :
1902 : /* -------------------------------------------------------------------- */
1903 : /* If the output dataset was created, and we have a destination */
1904 : /* nodata value, go through marking the bands with the information.*/
1905 : /* -------------------------------------------------------------------- */
1906 1299 : if (!psOptions->osDstNodata.empty() &&
1907 144 : !EQUAL(psOptions->osDstNodata.c_str(), "none"))
1908 : {
1909 : CPLStringList aosTokens(
1910 144 : CSLTokenizeString(psOptions->osDstNodata.c_str()));
1911 144 : const int nTokenCount = aosTokens.Count();
1912 144 : bool bDstNoDataNone = true;
1913 :
1914 144 : psWO->padfDstNoDataReal =
1915 144 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
1916 144 : psWO->padfDstNoDataImag =
1917 144 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
1918 :
1919 370 : for (int i = 0; i < psWO->nBandCount; i++)
1920 : {
1921 227 : psWO->padfDstNoDataReal[i] = -1.1e20;
1922 227 : psWO->padfDstNoDataImag[i] = 0.0;
1923 :
1924 227 : if (i < nTokenCount)
1925 : {
1926 145 : if (aosTokens[i] != nullptr && EQUAL(aosTokens[i], "none"))
1927 : {
1928 0 : CPLDebug("WARP", "dstnodata of band %d not set", i);
1929 0 : bDstNoDataNone = true;
1930 0 : continue;
1931 : }
1932 145 : else if (aosTokens[i] ==
1933 : nullptr) // this should not happen, but just in case
1934 : {
1935 0 : CPLError(CE_Failure, CPLE_AppDefined,
1936 : "Error parsing dstnodata arg #%d", i);
1937 0 : bDstNoDataNone = true;
1938 0 : continue;
1939 : }
1940 :
1941 145 : if (CPLStringToComplex(aosTokens[i],
1942 145 : psWO->padfDstNoDataReal + i,
1943 290 : psWO->padfDstNoDataImag + i) != CE_None)
1944 : {
1945 :
1946 1 : CPLError(CE_Failure, CPLE_AppDefined,
1947 : "Error parsing dstnodata for band %d", i + 1);
1948 1 : return CE_Failure;
1949 : }
1950 :
1951 288 : psWO->padfDstNoDataReal[i] =
1952 144 : GDALAdjustNoDataCloseToFloatMax(psWO->padfDstNoDataReal[i]);
1953 288 : psWO->padfDstNoDataImag[i] =
1954 144 : GDALAdjustNoDataCloseToFloatMax(psWO->padfDstNoDataImag[i]);
1955 144 : bDstNoDataNone = false;
1956 144 : CPLDebug("WARP", "dstnodata of band %d set to %f", i,
1957 144 : psWO->padfDstNoDataReal[i]);
1958 : }
1959 : else
1960 : {
1961 82 : if (!bDstNoDataNone)
1962 : {
1963 82 : psWO->padfDstNoDataReal[i] = psWO->padfDstNoDataReal[i - 1];
1964 82 : psWO->padfDstNoDataImag[i] = psWO->padfDstNoDataImag[i - 1];
1965 82 : CPLDebug("WARP",
1966 : "dstnodata of band %d set from previous band", i);
1967 : }
1968 : else
1969 : {
1970 0 : CPLDebug("WARP", "dstnodata value of band %d not set", i);
1971 0 : continue;
1972 : }
1973 : }
1974 :
1975 : GDALRasterBandH hBand =
1976 226 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]);
1977 226 : int bClamped = FALSE;
1978 226 : int bRounded = FALSE;
1979 226 : psWO->padfDstNoDataReal[i] = GDALAdjustValueToDataType(
1980 226 : GDALGetRasterDataType(hBand), psWO->padfDstNoDataReal[i],
1981 : &bClamped, &bRounded);
1982 :
1983 226 : if (bClamped)
1984 : {
1985 0 : CPLError(
1986 : CE_Warning, CPLE_AppDefined,
1987 : "for band %d, destination nodata value has been clamped "
1988 : "to %.0f, the original value being out of range.",
1989 0 : psWO->panDstBands[i], psWO->padfDstNoDataReal[i]);
1990 : }
1991 226 : else if (bRounded)
1992 : {
1993 0 : CPLError(
1994 : CE_Warning, CPLE_AppDefined,
1995 : "for band %d, destination nodata value has been rounded "
1996 : "to %.0f, %s being an integer datatype.",
1997 0 : psWO->panDstBands[i], psWO->padfDstNoDataReal[i],
1998 : GDALGetDataTypeName(GDALGetRasterDataType(hBand)));
1999 : }
2000 :
2001 226 : if (psOptions->bCreateOutput && iSrc == 0)
2002 : {
2003 222 : GDALSetRasterNoDataValue(
2004 222 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]),
2005 222 : psWO->padfDstNoDataReal[i]);
2006 : }
2007 : }
2008 : }
2009 :
2010 : /* check if the output dataset has already nodata */
2011 1154 : if (psOptions->osDstNodata.empty())
2012 : {
2013 1011 : int bHaveNodataAll = TRUE;
2014 2327 : for (int i = 0; i < psWO->nBandCount; i++)
2015 : {
2016 : GDALRasterBandH hBand =
2017 1316 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]);
2018 1316 : int bHaveNodata = FALSE;
2019 1316 : GDALGetRasterNoDataValue(hBand, &bHaveNodata);
2020 1316 : bHaveNodataAll &= bHaveNodata;
2021 : }
2022 1011 : if (bHaveNodataAll)
2023 : {
2024 6 : psWO->padfDstNoDataReal = static_cast<double *>(
2025 6 : CPLMalloc(psWO->nBandCount * sizeof(double)));
2026 13 : for (int i = 0; i < psWO->nBandCount; i++)
2027 : {
2028 : GDALRasterBandH hBand =
2029 7 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]);
2030 7 : int bHaveNodata = FALSE;
2031 14 : psWO->padfDstNoDataReal[i] =
2032 7 : GDALGetRasterNoDataValue(hBand, &bHaveNodata);
2033 7 : CPLDebug("WARP", "band=%d dstNoData=%f", i,
2034 7 : psWO->padfDstNoDataReal[i]);
2035 : }
2036 : }
2037 : }
2038 :
2039 : // If creating a new file that has default nodata value,
2040 : // try to override the default output nodata values with the source ones.
2041 2165 : if (psOptions->osDstNodata.empty() && psWO->padfSrcNoDataReal != nullptr &&
2042 99 : psWO->padfDstNoDataReal != nullptr && psOptions->bCreateOutput &&
2043 2165 : iSrc == 0 && !bEnableDstAlpha)
2044 : {
2045 5 : for (int i = 0; i < psWO->nBandCount; i++)
2046 : {
2047 : GDALRasterBandH hBand =
2048 3 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]);
2049 3 : int bHaveNodata = FALSE;
2050 3 : CPLPushErrorHandler(CPLQuietErrorHandler);
2051 : bool bRedefinedOK =
2052 3 : (GDALSetRasterNoDataValue(hBand, psWO->padfSrcNoDataReal[i]) ==
2053 3 : CE_None &&
2054 3 : GDALGetRasterNoDataValue(hBand, &bHaveNodata) ==
2055 9 : psWO->padfSrcNoDataReal[i] &&
2056 3 : bHaveNodata);
2057 3 : CPLPopErrorHandler();
2058 3 : if (bRedefinedOK)
2059 : {
2060 3 : if (i == 0 && !psOptions->bQuiet)
2061 0 : printf("Copying nodata values from source %s "
2062 : "to destination %s.\n",
2063 : GDALGetDescription(hSrcDS), pszDest);
2064 3 : psWO->padfDstNoDataReal[i] = psWO->padfSrcNoDataReal[i];
2065 :
2066 3 : if (i == 0 && !bInitDestSetByUser)
2067 : {
2068 : /* As we didn't know at the beginning if there was source
2069 : * nodata */
2070 : /* we have initialized INIT_DEST=0. Override this with
2071 : * NO_DATA now */
2072 2 : psWO->papszWarpOptions = CSLSetNameValue(
2073 : psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
2074 : }
2075 : }
2076 : else
2077 : {
2078 0 : break;
2079 : }
2080 : }
2081 : }
2082 :
2083 : /* else try to fill dstNoData from source bands, unless -dstalpha is
2084 : * specified */
2085 1152 : else if (psOptions->osDstNodata.empty() &&
2086 1009 : psWO->padfSrcNoDataReal != nullptr &&
2087 2161 : psWO->padfDstNoDataReal == nullptr && !bEnableDstAlpha)
2088 : {
2089 82 : psWO->padfDstNoDataReal =
2090 82 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
2091 :
2092 82 : if (psWO->padfSrcNoDataImag != nullptr)
2093 : {
2094 1 : psWO->padfDstNoDataImag = static_cast<double *>(
2095 1 : CPLMalloc(psWO->nBandCount * sizeof(double)));
2096 : }
2097 :
2098 82 : if (!psOptions->bQuiet)
2099 3 : printf("Copying nodata values from source %s to destination %s.\n",
2100 : GDALGetDescription(hSrcDS), pszDest);
2101 :
2102 183 : for (int i = 0; i < psWO->nBandCount; i++)
2103 : {
2104 101 : psWO->padfDstNoDataReal[i] = psWO->padfSrcNoDataReal[i];
2105 101 : if (psWO->padfSrcNoDataImag != nullptr)
2106 : {
2107 1 : psWO->padfDstNoDataImag[i] = psWO->padfSrcNoDataImag[i];
2108 : }
2109 101 : CPLDebug("WARP", "srcNoData=%f dstNoData=%f",
2110 101 : psWO->padfSrcNoDataReal[i], psWO->padfDstNoDataReal[i]);
2111 :
2112 101 : if (psOptions->bCreateOutput && iSrc == 0)
2113 : {
2114 101 : CPLDebug("WARP",
2115 : "calling GDALSetRasterNoDataValue() for band#%d", i);
2116 101 : GDALSetRasterNoDataValue(
2117 101 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]),
2118 101 : psWO->padfDstNoDataReal[i]);
2119 : }
2120 : }
2121 :
2122 82 : if (psOptions->bCreateOutput && !bInitDestSetByUser && iSrc == 0)
2123 : {
2124 : /* As we didn't know at the beginning if there was source nodata */
2125 : /* we have initialized INIT_DEST=0. Override this with NO_DATA now
2126 : */
2127 79 : psWO->papszWarpOptions =
2128 79 : CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
2129 : }
2130 : }
2131 :
2132 1154 : return CE_None;
2133 : }
2134 :
2135 : /************************************************************************/
2136 : /* SetupSkipNoSource() */
2137 : /************************************************************************/
2138 :
2139 1154 : static void SetupSkipNoSource(int iSrc, GDALDatasetH hDstDS,
2140 : GDALWarpOptions *psWO,
2141 : GDALWarpAppOptions *psOptions)
2142 : {
2143 1063 : if (psOptions->bCreateOutput && iSrc == 0 &&
2144 1041 : CSLFetchNameValue(psWO->papszWarpOptions, "SKIP_NOSOURCE") == nullptr &&
2145 1036 : CSLFetchNameValue(psWO->papszWarpOptions, "STREAMABLE_OUTPUT") ==
2146 2217 : nullptr &&
2147 : // This white list of drivers could potentially be extended.
2148 1035 : (EQUAL(psOptions->osFormat.c_str(), "MEM") ||
2149 629 : EQUAL(psOptions->osFormat.c_str(), "GTiff") ||
2150 185 : EQUAL(psOptions->osFormat.c_str(), "GPKG")))
2151 : {
2152 : // We can enable the optimization only if the user didn't specify
2153 : // a INIT_DEST value that would contradict the destination nodata.
2154 :
2155 850 : bool bOKRegardingInitDest = false;
2156 : const char *pszInitDest =
2157 850 : CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST");
2158 850 : if (pszInitDest == nullptr || EQUAL(pszInitDest, "NO_DATA"))
2159 : {
2160 176 : bOKRegardingInitDest = true;
2161 :
2162 : // The MEM driver will return non-initialized blocks at 0
2163 : // so make sure that the nodata value is 0.
2164 176 : if (EQUAL(psOptions->osFormat.c_str(), "MEM"))
2165 : {
2166 347 : for (int i = 0; i < GDALGetRasterCount(hDstDS); i++)
2167 : {
2168 235 : int bHasNoData = false;
2169 235 : double dfDstNoDataVal = GDALGetRasterNoDataValue(
2170 : GDALGetRasterBand(hDstDS, i + 1), &bHasNoData);
2171 235 : if (bHasNoData && dfDstNoDataVal != 0.0)
2172 : {
2173 39 : bOKRegardingInitDest = false;
2174 39 : break;
2175 : }
2176 : }
2177 176 : }
2178 : }
2179 : else
2180 : {
2181 674 : char **papszTokensInitDest = CSLTokenizeString(pszInitDest);
2182 674 : const int nTokenCountInitDest = CSLCount(papszTokensInitDest);
2183 674 : if (nTokenCountInitDest == 1 ||
2184 0 : nTokenCountInitDest == GDALGetRasterCount(hDstDS))
2185 : {
2186 674 : bOKRegardingInitDest = true;
2187 1567 : for (int i = 0; i < GDALGetRasterCount(hDstDS); i++)
2188 : {
2189 901 : double dfInitVal = GDALAdjustNoDataCloseToFloatMax(
2190 901 : CPLAtofM(papszTokensInitDest[std::min(
2191 901 : i, nTokenCountInitDest - 1)]));
2192 901 : int bHasNoData = false;
2193 901 : double dfDstNoDataVal = GDALGetRasterNoDataValue(
2194 : GDALGetRasterBand(hDstDS, i + 1), &bHasNoData);
2195 901 : if (!((bHasNoData && dfInitVal == dfDstNoDataVal) ||
2196 899 : (!bHasNoData && dfInitVal == 0.0)))
2197 : {
2198 7 : bOKRegardingInitDest = false;
2199 8 : break;
2200 : }
2201 894 : if (EQUAL(psOptions->osFormat.c_str(), "MEM") &&
2202 894 : bHasNoData && dfDstNoDataVal != 0.0)
2203 : {
2204 1 : bOKRegardingInitDest = false;
2205 1 : break;
2206 : }
2207 : }
2208 : }
2209 674 : CSLDestroy(papszTokensInitDest);
2210 : }
2211 :
2212 850 : if (bOKRegardingInitDest)
2213 : {
2214 803 : CPLDebug("GDALWARP", "Defining SKIP_NOSOURCE=YES");
2215 803 : psWO->papszWarpOptions =
2216 803 : CSLSetNameValue(psWO->papszWarpOptions, "SKIP_NOSOURCE", "YES");
2217 : }
2218 : }
2219 1154 : }
2220 :
2221 : /************************************************************************/
2222 : /* AdjustOutputExtentForRPC() */
2223 : /************************************************************************/
2224 :
2225 : /** Returns false if there's no intersection between source extent defined
2226 : * by RPC and target extent.
2227 : */
2228 1154 : static bool AdjustOutputExtentForRPC(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
2229 : GDALTransformerFunc pfnTransformer,
2230 : void *hTransformArg, GDALWarpOptions *psWO,
2231 : GDALWarpAppOptions *psOptions,
2232 : int &nWarpDstXOff, int &nWarpDstYOff,
2233 : int &nWarpDstXSize, int &nWarpDstYSize)
2234 : {
2235 1154 : if (CPLTestBool(CSLFetchNameValueDef(psWO->papszWarpOptions,
2236 919 : "SKIP_NOSOURCE", "NO")) &&
2237 919 : GDALGetMetadata(hSrcDS, "RPC") != nullptr &&
2238 2085 : EQUAL(FetchSrcMethod(psOptions->aosTransformerOptions, "RPC"), "RPC") &&
2239 12 : CPLTestBool(
2240 : CPLGetConfigOption("RESTRICT_OUTPUT_DATASET_UPDATE", "YES")))
2241 : {
2242 : double adfSuggestedGeoTransform[6];
2243 : double adfExtent[4];
2244 : int nPixels, nLines;
2245 10 : if (GDALSuggestedWarpOutput2(hSrcDS, pfnTransformer, hTransformArg,
2246 : adfSuggestedGeoTransform, &nPixels,
2247 10 : &nLines, adfExtent, 0) == CE_None)
2248 : {
2249 6 : const double dfMinX = adfExtent[0];
2250 6 : const double dfMinY = adfExtent[1];
2251 6 : const double dfMaxX = adfExtent[2];
2252 6 : const double dfMaxY = adfExtent[3];
2253 6 : const double dfThreshold = static_cast<double>(INT_MAX) / 2;
2254 6 : if (std::fabs(dfMinX) < dfThreshold &&
2255 6 : std::fabs(dfMinY) < dfThreshold &&
2256 6 : std::fabs(dfMaxX) < dfThreshold &&
2257 6 : std::fabs(dfMaxY) < dfThreshold)
2258 : {
2259 6 : const int nPadding = 5;
2260 6 : nWarpDstXOff =
2261 6 : std::max(nWarpDstXOff,
2262 6 : static_cast<int>(std::floor(dfMinX)) - nPadding);
2263 6 : nWarpDstYOff =
2264 6 : std::max(nWarpDstYOff,
2265 6 : static_cast<int>(std::floor(dfMinY)) - nPadding);
2266 12 : nWarpDstXSize = std::min(nWarpDstXSize - nWarpDstXOff,
2267 6 : static_cast<int>(std::ceil(dfMaxX)) +
2268 6 : nPadding - nWarpDstXOff);
2269 12 : nWarpDstYSize = std::min(nWarpDstYSize - nWarpDstYOff,
2270 6 : static_cast<int>(std::ceil(dfMaxY)) +
2271 6 : nPadding - nWarpDstYOff);
2272 6 : if (nWarpDstXSize <= 0 || nWarpDstYSize <= 0)
2273 : {
2274 1 : CPLDebug("WARP",
2275 : "No intersection between source extent defined "
2276 : "by RPC and target extent");
2277 1 : return false;
2278 : }
2279 2 : if (nWarpDstXOff != 0 || nWarpDstYOff != 0 ||
2280 9 : nWarpDstXSize != GDALGetRasterXSize(hDstDS) ||
2281 2 : nWarpDstYSize != GDALGetRasterYSize(hDstDS))
2282 : {
2283 3 : CPLDebug("WARP",
2284 : "Restricting warping to output dataset window "
2285 : "%d,%d,%dx%d",
2286 : nWarpDstXOff, nWarpDstYOff, nWarpDstXSize,
2287 : nWarpDstYSize);
2288 : }
2289 : }
2290 : }
2291 : }
2292 1153 : return true;
2293 : }
2294 :
2295 : /************************************************************************/
2296 : /* GDALWarpDirect() */
2297 : /************************************************************************/
2298 :
2299 : static GDALDatasetH
2300 1169 : GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS, int nSrcCount,
2301 : GDALDatasetH *pahSrcDS,
2302 : GDALTransformerArgUniquePtr hUniqueTransformArg,
2303 : GDALWarpAppOptions *psOptions, int *pbUsageError)
2304 : {
2305 1169 : CPLErrorReset();
2306 1169 : if (pszDest == nullptr && hDstDS == nullptr)
2307 : {
2308 0 : CPLError(CE_Failure, CPLE_AppDefined,
2309 : "pszDest == NULL && hDstDS == NULL");
2310 :
2311 0 : if (pbUsageError)
2312 0 : *pbUsageError = TRUE;
2313 0 : return nullptr;
2314 : }
2315 1169 : if (pszDest == nullptr)
2316 90 : pszDest = GDALGetDescription(hDstDS);
2317 :
2318 : #ifdef DEBUG
2319 1169 : GDALDataset *poDstDS = GDALDataset::FromHandle(hDstDS);
2320 : const int nExpectedRefCountAtEnd =
2321 1169 : (poDstDS != nullptr) ? poDstDS->GetRefCount() : 1;
2322 : (void)nExpectedRefCountAtEnd;
2323 : #endif
2324 1169 : const bool bDropDstDSRef = (hDstDS != nullptr);
2325 1169 : if (hDstDS != nullptr)
2326 93 : GDALReferenceDataset(hDstDS);
2327 :
2328 1169 : bool bVerticalShiftApplied = false;
2329 :
2330 1169 : if (psOptions->bNoVShift)
2331 : {
2332 1 : psOptions->aosTransformerOptions.SetNameValue("@STRIP_VERT_CS", "YES");
2333 : }
2334 1168 : else if (nSrcCount)
2335 : {
2336 1162 : bool bSrcHasVertAxis = false;
2337 1162 : bool bDstHasVertAxis = false;
2338 2324 : OGRSpatialReference oSRSSrc;
2339 2324 : OGRSpatialReference oSRSDst;
2340 :
2341 : bVerticalShiftApplied =
2342 1162 : MustApplyVerticalShift(pahSrcDS[0], psOptions, oSRSSrc, oSRSDst,
2343 : bSrcHasVertAxis, bDstHasVertAxis);
2344 1162 : if (bVerticalShiftApplied)
2345 : {
2346 : psOptions->aosTransformerOptions.SetNameValue("PROMOTE_TO_3D",
2347 20 : "YES");
2348 : }
2349 : }
2350 :
2351 1169 : bool bVRT = false;
2352 1169 : if (!CheckOptions(pszDest, hDstDS, nSrcCount, pahSrcDS, psOptions, bVRT,
2353 : pbUsageError))
2354 : {
2355 1 : return nullptr;
2356 : }
2357 :
2358 : /* -------------------------------------------------------------------- */
2359 : /* If we have a cutline datasource read it and attach it in the */
2360 : /* warp options. */
2361 : /* -------------------------------------------------------------------- */
2362 1168 : std::unique_ptr<OGRGeometry> poCutline;
2363 1168 : if (!ProcessCutlineOptions(nSrcCount, pahSrcDS, psOptions, poCutline))
2364 : {
2365 11 : return nullptr;
2366 : }
2367 :
2368 : /* -------------------------------------------------------------------- */
2369 : /* If the target dataset does not exist, we need to create it. */
2370 : /* -------------------------------------------------------------------- */
2371 : const bool bInitDestSetByUser =
2372 1157 : (psOptions->aosWarpOptions.FetchNameValue("INIT_DEST") != nullptr);
2373 :
2374 1157 : const bool bFigureoutCorrespondingWindow =
2375 2221 : (hDstDS != nullptr) ||
2376 1064 : (((psOptions->nForcePixels != 0 && psOptions->nForceLines != 0) ||
2377 896 : (psOptions->dfXRes != 0 && psOptions->dfYRes != 0)) &&
2378 261 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
2379 133 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0));
2380 :
2381 1157 : const char *pszMethod = FetchSrcMethod(psOptions->aosTransformerOptions);
2382 40 : if (pszMethod && EQUAL(pszMethod, "GCP_TPS") &&
2383 1202 : psOptions->dfErrorThreshold > 0 &&
2384 5 : !psOptions->aosTransformerOptions.FetchNameValue(
2385 : "SRC_APPROX_ERROR_IN_PIXEL"))
2386 : {
2387 : psOptions->aosTransformerOptions.SetNameValue(
2388 : "SRC_APPROX_ERROR_IN_PIXEL",
2389 5 : CPLSPrintf("%g", psOptions->dfErrorThreshold));
2390 : }
2391 :
2392 1157 : if (hDstDS == nullptr)
2393 : {
2394 1064 : hDstDS = CreateOutput(pszDest, nSrcCount, pahSrcDS, psOptions,
2395 : bInitDestSetByUser, hUniqueTransformArg);
2396 1064 : if (!hDstDS)
2397 : {
2398 19 : return nullptr;
2399 : }
2400 : #ifdef DEBUG
2401 : // Do not remove this if the #ifdef DEBUG before is still there !
2402 1045 : poDstDS = GDALDataset::FromHandle(hDstDS);
2403 1045 : CPL_IGNORE_RET_VAL(poDstDS);
2404 : #endif
2405 : }
2406 : else
2407 : {
2408 93 : if (psOptions->aosWarpOptions.FetchNameValue("SKIP_NOSOURCE") ==
2409 : nullptr)
2410 : {
2411 92 : CPLDebug("GDALWARP", "Defining SKIP_NOSOURCE=YES");
2412 92 : psOptions->aosWarpOptions.SetNameValue("SKIP_NOSOURCE", "YES");
2413 : }
2414 : }
2415 :
2416 : /* -------------------------------------------------------------------- */
2417 : /* Detect if output has alpha channel. */
2418 : /* -------------------------------------------------------------------- */
2419 1138 : bool bEnableDstAlpha = psOptions->bEnableDstAlpha;
2420 1039 : if (!bEnableDstAlpha && GDALGetRasterCount(hDstDS) &&
2421 1038 : GDALGetRasterColorInterpretation(GDALGetRasterBand(
2422 2177 : hDstDS, GDALGetRasterCount(hDstDS))) == GCI_AlphaBand &&
2423 34 : !psOptions->bDisableSrcAlpha)
2424 : {
2425 33 : if (!psOptions->bQuiet)
2426 1 : printf("Using band %d of destination image as alpha.\n",
2427 : GDALGetRasterCount(hDstDS));
2428 :
2429 33 : bEnableDstAlpha = true;
2430 : }
2431 :
2432 : /* -------------------------------------------------------------------- */
2433 : /* Create global progress function. */
2434 : /* -------------------------------------------------------------------- */
2435 : struct Progress
2436 : {
2437 : GDALProgressFunc pfnExternalProgress;
2438 : void *pExternalProgressData;
2439 : int iSrc;
2440 : int nSrcCount;
2441 : GDALDatasetH *pahSrcDS;
2442 :
2443 18993 : int Do(double dfComplete)
2444 : {
2445 37986 : CPLString osMsg;
2446 : osMsg.Printf("Processing %s [%d/%d]",
2447 18993 : CPLGetFilename(GDALGetDescription(pahSrcDS[iSrc])),
2448 18993 : iSrc + 1, nSrcCount);
2449 18993 : return pfnExternalProgress((iSrc + dfComplete) / nSrcCount,
2450 37986 : osMsg.c_str(), pExternalProgressData);
2451 : }
2452 :
2453 17664 : static int CPL_STDCALL ProgressFunc(double dfComplete, const char *,
2454 : void *pThis)
2455 : {
2456 17664 : return static_cast<Progress *>(pThis)->Do(dfComplete);
2457 : }
2458 : };
2459 :
2460 : Progress oProgress;
2461 1138 : oProgress.pfnExternalProgress = psOptions->pfnProgress;
2462 1138 : oProgress.pExternalProgressData = psOptions->pProgressData;
2463 1138 : oProgress.nSrcCount = nSrcCount;
2464 1138 : oProgress.pahSrcDS = pahSrcDS;
2465 :
2466 : /* -------------------------------------------------------------------- */
2467 : /* Loop over all source files, processing each in turn. */
2468 : /* -------------------------------------------------------------------- */
2469 1138 : bool bHasGotErr = false;
2470 2114 : for (int iSrc = 0; iSrc < nSrcCount; iSrc++)
2471 : {
2472 : GDALDatasetH hSrcDS;
2473 :
2474 : /* --------------------------------------------------------------------
2475 : */
2476 : /* Open this file. */
2477 : /* --------------------------------------------------------------------
2478 : */
2479 1160 : hSrcDS = pahSrcDS[iSrc];
2480 1160 : oProgress.iSrc = iSrc;
2481 :
2482 : /* --------------------------------------------------------------------
2483 : */
2484 : /* Check that there's at least one raster band */
2485 : /* --------------------------------------------------------------------
2486 : */
2487 1160 : if (GDALGetRasterCount(hSrcDS) == 0)
2488 : {
2489 1 : CPLError(CE_Failure, CPLE_AppDefined,
2490 : "Input file %s has no raster bands.",
2491 : GDALGetDescription(hSrcDS));
2492 1 : GDALReleaseDataset(hDstDS);
2493 184 : return nullptr;
2494 : }
2495 :
2496 : /* --------------------------------------------------------------------
2497 : */
2498 : /* Do we have a source alpha band? */
2499 : /* --------------------------------------------------------------------
2500 : */
2501 1159 : bool bEnableSrcAlpha = psOptions->bEnableSrcAlpha;
2502 1159 : if (GDALGetRasterColorInterpretation(GDALGetRasterBand(
2503 73 : hSrcDS, GDALGetRasterCount(hSrcDS))) == GCI_AlphaBand &&
2504 1159 : !bEnableSrcAlpha && !psOptions->bDisableSrcAlpha)
2505 : {
2506 31 : bEnableSrcAlpha = true;
2507 31 : if (!psOptions->bQuiet)
2508 0 : printf("Using band %d of source image as alpha.\n",
2509 : GDALGetRasterCount(hSrcDS));
2510 : }
2511 :
2512 : /* --------------------------------------------------------------------
2513 : */
2514 : /* Get the metadata of the first source DS and copy it to the */
2515 : /* destination DS. Copy Band-level metadata and other info, only */
2516 : /* if source and destination band count are equal. Any values that
2517 : */
2518 : /* conflict between source datasets are set to pszMDConflictValue.
2519 : */
2520 : /* --------------------------------------------------------------------
2521 : */
2522 1159 : ProcessMetadata(iSrc, hSrcDS, hDstDS, psOptions, bEnableDstAlpha,
2523 : bVerticalShiftApplied);
2524 :
2525 : /* --------------------------------------------------------------------
2526 : */
2527 : /* Warns if the file has a color table and something more */
2528 : /* complicated than nearest neighbour resampling is asked */
2529 : /* --------------------------------------------------------------------
2530 : */
2531 :
2532 2845 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
2533 1645 : psOptions->eResampleAlg != GRA_Mode &&
2534 486 : GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1)) != nullptr)
2535 : {
2536 0 : if (!psOptions->bQuiet)
2537 0 : CPLError(
2538 : CE_Warning, CPLE_AppDefined,
2539 : "Input file %s has a color table, which will likely lead "
2540 : "to "
2541 : "bad results when using a resampling method other than "
2542 : "nearest neighbour or mode. Converting the dataset prior "
2543 : "to 24/32 bit "
2544 : "is advised.",
2545 : GDALGetDescription(hSrcDS));
2546 : }
2547 :
2548 : // For RPC warping add a few extra source pixels by default
2549 : // (probably mostly needed in the RPC DEM case)
2550 1162 : if (iSrc == 0 && (GDALGetMetadata(hSrcDS, "RPC") != nullptr &&
2551 3 : (pszMethod == nullptr || EQUAL(pszMethod, "RPC"))))
2552 : {
2553 13 : if (!psOptions->aosWarpOptions.FetchNameValue("SOURCE_EXTRA"))
2554 : {
2555 13 : CPLDebug(
2556 : "WARP",
2557 : "Set SOURCE_EXTRA=5 warping options due to RPC warping");
2558 13 : psOptions->aosWarpOptions.SetNameValue("SOURCE_EXTRA", "5");
2559 : }
2560 :
2561 13 : if (!psOptions->aosWarpOptions.FetchNameValue("SAMPLE_STEPS") &&
2562 26 : !psOptions->aosWarpOptions.FetchNameValue("SAMPLE_GRID") &&
2563 13 : psOptions->aosTransformerOptions.FetchNameValue("RPC_DEM"))
2564 : {
2565 10 : CPLDebug("WARP", "Set SAMPLE_STEPS=ALL warping options due to "
2566 : "RPC DEM warping");
2567 10 : psOptions->aosWarpOptions.SetNameValue("SAMPLE_STEPS", "ALL");
2568 : }
2569 : }
2570 : // Also do the same for GCP TPS warping, e.g. to solve use case of
2571 : // https://github.com/OSGeo/gdal/issues/12736
2572 2268 : else if (iSrc == 0 &&
2573 1122 : (GDALGetGCPCount(hSrcDS) > 0 &&
2574 6 : (pszMethod == nullptr || EQUAL(pszMethod, "TPS"))))
2575 : {
2576 21 : if (!psOptions->aosWarpOptions.FetchNameValue("SOURCE_EXTRA"))
2577 : {
2578 21 : CPLDebug(
2579 : "WARP",
2580 : "Set SOURCE_EXTRA=5 warping options due to TPS warping");
2581 : }
2582 : }
2583 :
2584 1159 : if (iSrc > 0)
2585 : psOptions->aosWarpOptions.SetNameValue("RESET_DEST_PIXELS",
2586 24 : nullptr);
2587 :
2588 : /* --------------------------------------------------------------------
2589 : */
2590 : /* Create a transformation object from the source to */
2591 : /* destination coordinate system. */
2592 : /* --------------------------------------------------------------------
2593 : */
2594 0 : GDALTransformerArgUniquePtr hTransformArg;
2595 1159 : if (hUniqueTransformArg)
2596 1023 : hTransformArg = std::move(hUniqueTransformArg);
2597 : else
2598 : {
2599 136 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
2600 136 : hSrcDS, hDstDS, psOptions->aosTransformerOptions.List()));
2601 136 : if (hTransformArg == nullptr)
2602 : {
2603 0 : GDALReleaseDataset(hDstDS);
2604 0 : return nullptr;
2605 : }
2606 : }
2607 :
2608 1159 : GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;
2609 :
2610 : // Check if transformation is inversible
2611 : {
2612 1159 : double dfX = GDALGetRasterXSize(hDstDS) / 2.0;
2613 1159 : double dfY = GDALGetRasterYSize(hDstDS) / 2.0;
2614 1159 : double dfZ = 0;
2615 1159 : int bSuccess = false;
2616 1159 : const auto nErrorCounterBefore = CPLGetErrorCounter();
2617 1159 : pfnTransformer(hTransformArg.get(), TRUE, 1, &dfX, &dfY, &dfZ,
2618 : &bSuccess);
2619 1159 : if (!bSuccess && CPLGetErrorCounter() > nErrorCounterBefore &&
2620 0 : strstr(CPLGetLastErrorMsg(), "No inverse operation"))
2621 : {
2622 0 : GDALReleaseDataset(hDstDS);
2623 0 : return nullptr;
2624 : }
2625 : }
2626 :
2627 : /* --------------------------------------------------------------------
2628 : */
2629 : /* Determine if we must work with the full-resolution source */
2630 : /* dataset, or one of its overview level. */
2631 : /* --------------------------------------------------------------------
2632 : */
2633 1159 : GDALDataset *poSrcDS = static_cast<GDALDataset *>(hSrcDS);
2634 1159 : GDALDataset *poSrcOvrDS = nullptr;
2635 1159 : int nOvCount = poSrcDS->GetRasterBand(1)->GetOverviewCount();
2636 1159 : if (psOptions->nOvLevel <= OVR_LEVEL_AUTO && nOvCount > 0)
2637 : {
2638 22 : double dfTargetRatio = 0;
2639 22 : double dfTargetRatioX = 0;
2640 22 : double dfTargetRatioY = 0;
2641 :
2642 22 : if (bFigureoutCorrespondingWindow)
2643 : {
2644 : // If the user has explicitly set the target bounds and
2645 : // resolution, or we're updating an existing file, then figure
2646 : // out which source window corresponds to the target raster.
2647 4 : constexpr int nPointsOneDim = 10;
2648 4 : constexpr int nPoints = nPointsOneDim * nPointsOneDim;
2649 8 : std::vector<double> adfX(nPoints);
2650 8 : std::vector<double> adfY(nPoints);
2651 8 : std::vector<double> adfZ(nPoints);
2652 4 : const int nDstXSize = GDALGetRasterXSize(hDstDS);
2653 4 : const int nDstYSize = GDALGetRasterYSize(hDstDS);
2654 4 : int iPoint = 0;
2655 44 : for (int iX = 0; iX < nPointsOneDim; ++iX)
2656 : {
2657 440 : for (int iY = 0; iY < nPointsOneDim; ++iY)
2658 : {
2659 400 : adfX[iPoint] = nDstXSize * static_cast<double>(iX) /
2660 : (nPointsOneDim - 1);
2661 400 : adfY[iPoint] = nDstYSize * static_cast<double>(iY) /
2662 : (nPointsOneDim - 1);
2663 400 : iPoint++;
2664 : }
2665 : }
2666 4 : std::vector<int> abSuccess(nPoints);
2667 4 : pfnTransformer(hTransformArg.get(), TRUE, nPoints, &adfX[0],
2668 4 : &adfY[0], &adfZ[0], &abSuccess[0]);
2669 :
2670 4 : double dfMinSrcX = std::numeric_limits<double>::infinity();
2671 4 : double dfMaxSrcX = -std::numeric_limits<double>::infinity();
2672 4 : double dfMinSrcY = std::numeric_limits<double>::infinity();
2673 4 : double dfMaxSrcY = -std::numeric_limits<double>::infinity();
2674 404 : for (int i = 0; i < nPoints; i++)
2675 : {
2676 400 : if (abSuccess[i])
2677 : {
2678 400 : dfMinSrcX = std::min(dfMinSrcX, adfX[i]);
2679 400 : dfMaxSrcX = std::max(dfMaxSrcX, adfX[i]);
2680 400 : dfMinSrcY = std::min(dfMinSrcY, adfY[i]);
2681 400 : dfMaxSrcY = std::max(dfMaxSrcY, adfY[i]);
2682 : }
2683 : }
2684 4 : if (dfMaxSrcX > dfMinSrcX)
2685 : {
2686 4 : dfTargetRatioX =
2687 4 : (dfMaxSrcX - dfMinSrcX) / GDALGetRasterXSize(hDstDS);
2688 : }
2689 4 : if (dfMaxSrcY > dfMinSrcY)
2690 : {
2691 4 : dfTargetRatioY =
2692 4 : (dfMaxSrcY - dfMinSrcY) / GDALGetRasterYSize(hDstDS);
2693 : }
2694 : // take the minimum of these ratios #7019
2695 4 : dfTargetRatio = std::min(dfTargetRatioX, dfTargetRatioY);
2696 : }
2697 : else
2698 : {
2699 : /* Compute what the "natural" output resolution (in pixels)
2700 : * would be for this */
2701 : /* input dataset */
2702 : double adfSuggestedGeoTransform[6];
2703 : int nPixels, nLines;
2704 18 : if (GDALSuggestedWarpOutput(
2705 : hSrcDS, pfnTransformer, hTransformArg.get(),
2706 18 : adfSuggestedGeoTransform, &nPixels, &nLines) == CE_None)
2707 : {
2708 18 : dfTargetRatio = 1.0 / adfSuggestedGeoTransform[1];
2709 : }
2710 : }
2711 :
2712 22 : if (dfTargetRatio > 1.0)
2713 : {
2714 15 : const char *pszOversampligThreshold = CPLGetConfigOption(
2715 : "GDALWARP_OVERSAMPLING_THRESHOLD", nullptr);
2716 : const double dfOversamplingThreshold =
2717 15 : pszOversampligThreshold ? CPLAtof(pszOversampligThreshold)
2718 15 : : 1.0;
2719 :
2720 15 : const int iBestOvr = GDALBandGetBestOverviewLevel(
2721 : poSrcDS->GetRasterBand(1), dfTargetRatio,
2722 : dfOversamplingThreshold);
2723 15 : const int iOvr =
2724 15 : iBestOvr + (psOptions->nOvLevel - OVR_LEVEL_AUTO);
2725 15 : if (iOvr >= 0)
2726 : {
2727 9 : CPLDebug("WARP", "Selecting overview level %d for %s", iOvr,
2728 : GDALGetDescription(hSrcDS));
2729 : poSrcOvrDS =
2730 9 : GDALCreateOverviewDataset(poSrcDS, iOvr,
2731 : /* bThisLevelOnly = */ false);
2732 : }
2733 22 : }
2734 : }
2735 1137 : else if (psOptions->nOvLevel >= 0)
2736 : {
2737 6 : poSrcOvrDS = GDALCreateOverviewDataset(poSrcDS, psOptions->nOvLevel,
2738 : /* bThisLevelOnly = */ true);
2739 6 : if (poSrcOvrDS == nullptr)
2740 : {
2741 1 : if (!psOptions->bQuiet)
2742 : {
2743 1 : CPLError(CE_Warning, CPLE_AppDefined,
2744 : "cannot get overview level %d for "
2745 : "dataset %s. Defaulting to level %d",
2746 : psOptions->nOvLevel, GDALGetDescription(hSrcDS),
2747 : nOvCount - 1);
2748 : }
2749 1 : if (nOvCount > 0)
2750 : poSrcOvrDS =
2751 1 : GDALCreateOverviewDataset(poSrcDS, nOvCount - 1,
2752 : /* bThisLevelOnly = */ false);
2753 : }
2754 : else
2755 : {
2756 5 : CPLDebug("WARP", "Selecting overview level %d for %s",
2757 : psOptions->nOvLevel, GDALGetDescription(hSrcDS));
2758 : }
2759 : }
2760 :
2761 1159 : if (poSrcOvrDS == nullptr)
2762 1144 : GDALReferenceDataset(hSrcDS);
2763 :
2764 1159 : GDALDatasetH hWrkSrcDS =
2765 1159 : poSrcOvrDS ? static_cast<GDALDatasetH>(poSrcOvrDS) : hSrcDS;
2766 :
2767 : /* --------------------------------------------------------------------
2768 : */
2769 : /* Clear temporary INIT_DEST settings after the first image. */
2770 : /* --------------------------------------------------------------------
2771 : */
2772 1159 : if (psOptions->bCreateOutput && iSrc == 1)
2773 22 : psOptions->aosWarpOptions.SetNameValue("INIT_DEST", nullptr);
2774 :
2775 : /* --------------------------------------------------------------------
2776 : */
2777 : /* Define SKIP_NOSOURCE after the first image (since
2778 : * initialization*/
2779 : /* has already be done). */
2780 : /* --------------------------------------------------------------------
2781 : */
2782 1159 : if (iSrc == 1 && psOptions->aosWarpOptions.FetchNameValue(
2783 : "SKIP_NOSOURCE") == nullptr)
2784 : {
2785 22 : CPLDebug("GDALWARP", "Defining SKIP_NOSOURCE=YES");
2786 22 : psOptions->aosWarpOptions.SetNameValue("SKIP_NOSOURCE", "YES");
2787 : }
2788 :
2789 : /* --------------------------------------------------------------------
2790 : */
2791 : /* Setup warp options. */
2792 : /* --------------------------------------------------------------------
2793 : */
2794 : std::unique_ptr<GDALWarpOptions, decltype(&GDALDestroyWarpOptions)>
2795 1159 : psWO(GDALCreateWarpOptions(), GDALDestroyWarpOptions);
2796 :
2797 1159 : psWO->papszWarpOptions = CSLDuplicate(psOptions->aosWarpOptions.List());
2798 1159 : psWO->eWorkingDataType = psOptions->eWorkingType;
2799 :
2800 1159 : psWO->eResampleAlg = psOptions->eResampleAlg;
2801 :
2802 1159 : psWO->hSrcDS = hWrkSrcDS;
2803 1159 : psWO->hDstDS = hDstDS;
2804 :
2805 1159 : if (!bVRT)
2806 : {
2807 982 : if (psOptions->pfnProgress == GDALDummyProgress)
2808 : {
2809 860 : psWO->pfnProgress = GDALDummyProgress;
2810 860 : psWO->pProgressArg = nullptr;
2811 : }
2812 : else
2813 : {
2814 122 : psWO->pfnProgress = Progress::ProgressFunc;
2815 122 : psWO->pProgressArg = &oProgress;
2816 : }
2817 : }
2818 :
2819 1159 : if (psOptions->dfWarpMemoryLimit != 0.0)
2820 35 : psWO->dfWarpMemoryLimit = psOptions->dfWarpMemoryLimit;
2821 :
2822 : /* --------------------------------------------------------------------
2823 : */
2824 : /* Setup band mapping. */
2825 : /* --------------------------------------------------------------------
2826 : */
2827 1159 : if (psOptions->anSrcBands.empty())
2828 : {
2829 1132 : if (bEnableSrcAlpha)
2830 62 : psWO->nBandCount = GDALGetRasterCount(hWrkSrcDS) - 1;
2831 : else
2832 1070 : psWO->nBandCount = GDALGetRasterCount(hWrkSrcDS);
2833 : }
2834 : else
2835 : {
2836 27 : psWO->nBandCount = static_cast<int>(psOptions->anSrcBands.size());
2837 : }
2838 :
2839 : const int nNeededDstBands =
2840 1159 : psWO->nBandCount + (bEnableDstAlpha ? 1 : 0);
2841 1159 : if (nNeededDstBands > GDALGetRasterCount(hDstDS))
2842 : {
2843 1 : CPLError(CE_Failure, CPLE_AppDefined,
2844 : "Destination dataset has %d bands, but at least %d "
2845 : "are needed",
2846 : GDALGetRasterCount(hDstDS), nNeededDstBands);
2847 1 : GDALReleaseDataset(hWrkSrcDS);
2848 1 : GDALReleaseDataset(hDstDS);
2849 1 : return nullptr;
2850 : }
2851 :
2852 2316 : psWO->panSrcBands =
2853 1158 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
2854 2316 : psWO->panDstBands =
2855 1158 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
2856 1158 : if (psOptions->anSrcBands.empty())
2857 : {
2858 2634 : for (int i = 0; i < psWO->nBandCount; i++)
2859 : {
2860 1503 : psWO->panSrcBands[i] = i + 1;
2861 1503 : psWO->panDstBands[i] = i + 1;
2862 : }
2863 : }
2864 : else
2865 : {
2866 68 : for (int i = 0; i < psWO->nBandCount; i++)
2867 : {
2868 86 : if (psOptions->anSrcBands[i] <= 0 ||
2869 43 : psOptions->anSrcBands[i] > GDALGetRasterCount(hSrcDS))
2870 : {
2871 1 : CPLError(CE_Failure, CPLE_AppDefined,
2872 : "-srcband[%d] = %d is invalid", i,
2873 1 : psOptions->anSrcBands[i]);
2874 1 : GDALReleaseDataset(hWrkSrcDS);
2875 1 : GDALReleaseDataset(hDstDS);
2876 1 : return nullptr;
2877 : }
2878 84 : if (psOptions->anDstBands[i] <= 0 ||
2879 42 : psOptions->anDstBands[i] > GDALGetRasterCount(hDstDS))
2880 : {
2881 1 : CPLError(CE_Failure, CPLE_AppDefined,
2882 : "-dstband[%d] = %d is invalid", i,
2883 1 : psOptions->anDstBands[i]);
2884 1 : GDALReleaseDataset(hWrkSrcDS);
2885 1 : GDALReleaseDataset(hDstDS);
2886 1 : return nullptr;
2887 : }
2888 41 : psWO->panSrcBands[i] = psOptions->anSrcBands[i];
2889 41 : psWO->panDstBands[i] = psOptions->anDstBands[i];
2890 : }
2891 : }
2892 :
2893 : /* --------------------------------------------------------------------
2894 : */
2895 : /* Setup alpha bands used if any. */
2896 : /* --------------------------------------------------------------------
2897 : */
2898 1156 : if (bEnableSrcAlpha)
2899 70 : psWO->nSrcAlphaBand = GDALGetRasterCount(hWrkSrcDS);
2900 :
2901 1156 : if (bEnableDstAlpha)
2902 : {
2903 137 : if (psOptions->anSrcBands.empty())
2904 128 : psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
2905 : else
2906 9 : psWO->nDstAlphaBand =
2907 9 : static_cast<int>(psOptions->anDstBands.size()) + 1;
2908 : }
2909 :
2910 : /* ------------------------------------------------------------------ */
2911 : /* Setup NODATA options. */
2912 : /* ------------------------------------------------------------------ */
2913 1156 : if (SetupNoData(pszDest, iSrc, hSrcDS, hWrkSrcDS, hDstDS, psWO.get(),
2914 : psOptions, bEnableDstAlpha,
2915 1156 : bInitDestSetByUser) != CE_None)
2916 : {
2917 2 : GDALReleaseDataset(hWrkSrcDS);
2918 2 : GDALReleaseDataset(hDstDS);
2919 2 : return nullptr;
2920 : }
2921 :
2922 1154 : oProgress.Do(0);
2923 :
2924 : /* --------------------------------------------------------------------
2925 : */
2926 : /* For the first source image of a newly created dataset, decide */
2927 : /* if we can safely enable SKIP_NOSOURCE optimization. */
2928 : /* --------------------------------------------------------------------
2929 : */
2930 1154 : SetupSkipNoSource(iSrc, hDstDS, psWO.get(), psOptions);
2931 :
2932 : /* --------------------------------------------------------------------
2933 : */
2934 : /* In some cases, RPC evaluation can find valid input pixel for */
2935 : /* output pixels that are outside the footprint of the source */
2936 : /* dataset, so limit the area we update in the target dataset from
2937 : */
2938 : /* the suggested warp output (only in cases where
2939 : * SKIP_NOSOURCE=YES) */
2940 : /* --------------------------------------------------------------------
2941 : */
2942 1154 : int nWarpDstXOff = 0;
2943 1154 : int nWarpDstYOff = 0;
2944 1154 : int nWarpDstXSize = GDALGetRasterXSize(hDstDS);
2945 1154 : int nWarpDstYSize = GDALGetRasterYSize(hDstDS);
2946 :
2947 1154 : if (!AdjustOutputExtentForRPC(hSrcDS, hDstDS, pfnTransformer,
2948 : hTransformArg.get(), psWO.get(),
2949 : psOptions, nWarpDstXOff, nWarpDstYOff,
2950 : nWarpDstXSize, nWarpDstYSize))
2951 : {
2952 1 : GDALReleaseDataset(hWrkSrcDS);
2953 1 : continue;
2954 : }
2955 :
2956 : /* We need to recreate the transform when operating on an overview */
2957 1153 : if (poSrcOvrDS != nullptr)
2958 : {
2959 15 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
2960 15 : hWrkSrcDS, hDstDS, psOptions->aosTransformerOptions.List()));
2961 : }
2962 :
2963 1153 : bool bUseApproxTransformer = psOptions->dfErrorThreshold != 0.0;
2964 :
2965 1153 : if (!psOptions->bNoVShift)
2966 : {
2967 : // Can modify psWO->papszWarpOptions
2968 1152 : if (ApplyVerticalShift(hWrkSrcDS, psOptions, psWO.get()))
2969 : {
2970 20 : bUseApproxTransformer = false;
2971 : }
2972 : }
2973 :
2974 : /* --------------------------------------------------------------------
2975 : */
2976 : /* Warp the transformer with a linear approximator unless the */
2977 : /* acceptable error is zero. */
2978 : /* --------------------------------------------------------------------
2979 : */
2980 1153 : if (bUseApproxTransformer)
2981 : {
2982 1119 : hTransformArg.reset(GDALCreateApproxTransformer(
2983 : GDALGenImgProjTransform, hTransformArg.release(),
2984 : psOptions->dfErrorThreshold));
2985 1119 : pfnTransformer = GDALApproxTransform;
2986 1119 : GDALApproxTransformerOwnsSubtransformer(hTransformArg.get(), TRUE);
2987 : }
2988 :
2989 : /* --------------------------------------------------------------------
2990 : */
2991 : /* If we have a cutline, transform it into the source */
2992 : /* pixel/line coordinate system and insert into warp options. */
2993 : /* --------------------------------------------------------------------
2994 : */
2995 1153 : if (poCutline)
2996 : {
2997 : CPLErr eError;
2998 46 : eError = TransformCutlineToSource(
2999 : GDALDataset::FromHandle(hWrkSrcDS), poCutline.get(),
3000 46 : &(psWO->papszWarpOptions),
3001 46 : psOptions->aosTransformerOptions.List());
3002 46 : if (eError == CE_Failure)
3003 : {
3004 1 : GDALReleaseDataset(hWrkSrcDS);
3005 1 : GDALReleaseDataset(hDstDS);
3006 1 : return nullptr;
3007 : }
3008 : }
3009 :
3010 : /* --------------------------------------------------------------------
3011 : */
3012 : /* If we are producing VRT output, then just initialize it with */
3013 : /* the warp options and write out now rather than proceeding */
3014 : /* with the operations. */
3015 : /* --------------------------------------------------------------------
3016 : */
3017 1152 : if (bVRT)
3018 : {
3019 177 : GDALSetMetadataItem(hDstDS, "SrcOvrLevel",
3020 : CPLSPrintf("%d", psOptions->nOvLevel), nullptr);
3021 :
3022 : // In case of success, hDstDS has become the owner of hTransformArg
3023 : // so we need to release it
3024 177 : psWO->pfnTransformer = pfnTransformer;
3025 177 : psWO->pTransformerArg = hTransformArg.release();
3026 177 : CPLErr eErr = GDALInitializeWarpedVRT(hDstDS, psWO.get());
3027 177 : if (eErr != CE_None)
3028 : {
3029 : // In case of error, reacquire psWO->pTransformerArg
3030 1 : hTransformArg.reset(psWO->pTransformerArg);
3031 : }
3032 177 : GDALReleaseDataset(hWrkSrcDS);
3033 177 : if (eErr != CE_None)
3034 : {
3035 1 : GDALReleaseDataset(hDstDS);
3036 1 : return nullptr;
3037 : }
3038 :
3039 176 : if (!EQUAL(pszDest, ""))
3040 : {
3041 : const bool bWasFailureBefore =
3042 20 : (CPLGetLastErrorType() == CE_Failure);
3043 20 : GDALFlushCache(hDstDS);
3044 20 : if (!bWasFailureBefore && CPLGetLastErrorType() == CE_Failure)
3045 : {
3046 1 : GDALReleaseDataset(hDstDS);
3047 1 : hDstDS = nullptr;
3048 : }
3049 : }
3050 :
3051 176 : if (hDstDS)
3052 175 : oProgress.Do(1);
3053 :
3054 176 : return hDstDS;
3055 : }
3056 :
3057 : /* --------------------------------------------------------------------
3058 : */
3059 : /* Initialize and execute the warp. */
3060 : /* --------------------------------------------------------------------
3061 : */
3062 1950 : GDALWarpOperation oWO;
3063 :
3064 1950 : if (oWO.Initialize(psWO.get(), pfnTransformer,
3065 1950 : std::move(hTransformArg)) == CE_None)
3066 : {
3067 : CPLErr eErr;
3068 972 : if (psOptions->bMulti)
3069 6 : eErr = oWO.ChunkAndWarpMulti(nWarpDstXOff, nWarpDstYOff,
3070 : nWarpDstXSize, nWarpDstYSize);
3071 : else
3072 966 : eErr = oWO.ChunkAndWarpImage(nWarpDstXOff, nWarpDstYOff,
3073 : nWarpDstXSize, nWarpDstYSize);
3074 972 : if (eErr != CE_None)
3075 5 : bHasGotErr = true;
3076 : }
3077 : else
3078 : {
3079 3 : bHasGotErr = true;
3080 : }
3081 :
3082 : /* --------------------------------------------------------------------
3083 : */
3084 : /* Cleanup */
3085 : /* --------------------------------------------------------------------
3086 : */
3087 975 : GDALReleaseDataset(hWrkSrcDS);
3088 : }
3089 :
3090 : /* -------------------------------------------------------------------- */
3091 : /* Final Cleanup. */
3092 : /* -------------------------------------------------------------------- */
3093 954 : const bool bWasFailureBefore = (CPLGetLastErrorType() == CE_Failure);
3094 954 : GDALFlushCache(hDstDS);
3095 954 : if (!bWasFailureBefore && CPLGetLastErrorType() == CE_Failure)
3096 : {
3097 1 : bHasGotErr = true;
3098 : }
3099 :
3100 954 : if (bHasGotErr || bDropDstDSRef)
3101 100 : GDALReleaseDataset(hDstDS);
3102 :
3103 : #ifdef DEBUG
3104 954 : if (!bHasGotErr || bDropDstDSRef)
3105 : {
3106 945 : CPLAssert(poDstDS->GetRefCount() == nExpectedRefCountAtEnd);
3107 : }
3108 : #endif
3109 :
3110 954 : return bHasGotErr ? nullptr : hDstDS;
3111 : }
3112 :
3113 : /************************************************************************/
3114 : /* ValidateCutline() */
3115 : /* Same as OGR_G_IsValid() except that it processes polygon per polygon*/
3116 : /* without paying attention to MultiPolygon specific validity rules. */
3117 : /************************************************************************/
3118 :
3119 234 : static bool ValidateCutline(const OGRGeometry *poGeom, bool bVerbose)
3120 : {
3121 234 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3122 234 : if (eType == wkbMultiPolygon)
3123 : {
3124 151 : for (const auto *poSubGeom : *(poGeom->toMultiPolygon()))
3125 : {
3126 81 : if (!ValidateCutline(poSubGeom, bVerbose))
3127 5 : return false;
3128 : }
3129 : }
3130 159 : else if (eType == wkbPolygon)
3131 : {
3132 158 : std::string osReason;
3133 158 : if (OGRGeometryFactory::haveGEOS() && !poGeom->IsValid(&osReason))
3134 : {
3135 6 : if (!bVerbose)
3136 6 : return false;
3137 :
3138 2 : char *pszWKT = nullptr;
3139 2 : poGeom->exportToWkt(&pszWKT);
3140 2 : CPLDebug("GDALWARP", "WKT = \"%s\"", pszWKT ? pszWKT : "(null)");
3141 : const char *pszFile =
3142 2 : CPLGetConfigOption("GDALWARP_DUMP_WKT_TO_FILE", nullptr);
3143 2 : if (pszFile && pszWKT)
3144 : {
3145 : FILE *f =
3146 0 : EQUAL(pszFile, "stderr") ? stderr : fopen(pszFile, "wb");
3147 0 : if (f)
3148 : {
3149 0 : fprintf(f, "id,WKT\n");
3150 0 : fprintf(f, "1,\"%s\"\n", pszWKT);
3151 0 : if (!EQUAL(pszFile, "stderr"))
3152 0 : fclose(f);
3153 : }
3154 : }
3155 2 : CPLFree(pszWKT);
3156 :
3157 2 : if (CPLTestBool(
3158 : CPLGetConfigOption("GDALWARP_IGNORE_BAD_CUTLINE", "NO")))
3159 : {
3160 0 : CPLError(CE_Warning, CPLE_AppDefined,
3161 : "Cutline polygon is invalid: %s.", osReason.c_str());
3162 : }
3163 : else
3164 : {
3165 2 : CPLError(CE_Failure, CPLE_AppDefined,
3166 : "Cutline polygon is invalid: %s.", osReason.c_str());
3167 2 : return false;
3168 : }
3169 : }
3170 : }
3171 : else
3172 : {
3173 1 : if (bVerbose)
3174 : {
3175 1 : CPLError(CE_Failure, CPLE_AppDefined,
3176 : "Cutline not of polygon type.");
3177 : }
3178 1 : return false;
3179 : }
3180 :
3181 222 : return true;
3182 : }
3183 :
3184 : /************************************************************************/
3185 : /* LoadCutline() */
3186 : /* */
3187 : /* Load blend cutline from OGR datasource. */
3188 : /************************************************************************/
3189 :
3190 58 : static CPLErr LoadCutline(const std::string &osCutlineDSNameOrWKT,
3191 : const std::string &osSRS, const std::string &osCLayer,
3192 : const std::string &osCWHERE,
3193 : const std::string &osCSQL, OGRGeometryH *phCutlineRet)
3194 :
3195 : {
3196 58 : if (STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "POLYGON(") ||
3197 58 : STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "POLYGON (") ||
3198 160 : STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "MULTIPOLYGON(") ||
3199 44 : STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "MULTIPOLYGON ("))
3200 : {
3201 28 : OGRSpatialReferenceRefCountedPtr poSRS;
3202 14 : if (!osSRS.empty())
3203 : {
3204 8 : poSRS = OGRSpatialReferenceRefCountedPtr::makeInstance();
3205 8 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3206 8 : poSRS->SetFromUserInput(osSRS.c_str());
3207 : }
3208 :
3209 14 : auto [poGeom, _] = OGRGeometryFactory::createFromWkt(
3210 14 : osCutlineDSNameOrWKT.c_str(), poSRS.get());
3211 14 : *phCutlineRet = OGRGeometry::ToHandle(poGeom.release());
3212 14 : return *phCutlineRet ? CE_None : CE_Failure;
3213 : }
3214 :
3215 : /* -------------------------------------------------------------------- */
3216 : /* Open source vector dataset. */
3217 : /* -------------------------------------------------------------------- */
3218 : auto poDS = std::unique_ptr<GDALDataset>(
3219 88 : GDALDataset::Open(osCutlineDSNameOrWKT.c_str(), GDAL_OF_VECTOR));
3220 44 : if (poDS == nullptr)
3221 : {
3222 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot open %s.",
3223 : osCutlineDSNameOrWKT.c_str());
3224 1 : return CE_Failure;
3225 : }
3226 :
3227 : /* -------------------------------------------------------------------- */
3228 : /* Get the source layer */
3229 : /* -------------------------------------------------------------------- */
3230 43 : OGRLayer *poLayer = nullptr;
3231 :
3232 43 : if (!osCSQL.empty())
3233 2 : poLayer = poDS->ExecuteSQL(osCSQL.c_str(), nullptr, nullptr);
3234 41 : else if (!osCLayer.empty())
3235 15 : poLayer = poDS->GetLayerByName(osCLayer.c_str());
3236 : else
3237 26 : poLayer = poDS->GetLayer(0);
3238 :
3239 43 : if (poLayer == nullptr)
3240 : {
3241 1 : CPLError(CE_Failure, CPLE_AppDefined,
3242 : "Failed to identify source layer from datasource.");
3243 1 : return CE_Failure;
3244 : }
3245 :
3246 : /* -------------------------------------------------------------------- */
3247 : /* Apply WHERE clause if there is one. */
3248 : /* -------------------------------------------------------------------- */
3249 42 : if (!osCWHERE.empty())
3250 1 : poLayer->SetAttributeFilter(osCWHERE.c_str());
3251 :
3252 : /* -------------------------------------------------------------------- */
3253 : /* Collect the geometries from this layer, and build list of */
3254 : /* burn values. */
3255 : /* -------------------------------------------------------------------- */
3256 84 : auto poMultiPolygon = std::make_unique<OGRMultiPolygon>();
3257 :
3258 81 : for (auto &&poFeature : poLayer)
3259 : {
3260 42 : auto poGeom = std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
3261 42 : if (poGeom == nullptr)
3262 : {
3263 1 : CPLError(CE_Failure, CPLE_AppDefined,
3264 : "Cutline feature without a geometry.");
3265 1 : goto error;
3266 : }
3267 :
3268 41 : if (!ValidateCutline(poGeom.get(), true))
3269 : {
3270 2 : goto error;
3271 : }
3272 :
3273 39 : OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3274 :
3275 39 : if (eType == wkbPolygon)
3276 36 : poMultiPolygon->addGeometry(std::move(poGeom));
3277 3 : else if (eType == wkbMultiPolygon)
3278 : {
3279 7 : for (const auto *poSubGeom : poGeom->toMultiPolygon())
3280 : {
3281 4 : poMultiPolygon->addGeometry(poSubGeom);
3282 : }
3283 : }
3284 : }
3285 :
3286 39 : if (poMultiPolygon->IsEmpty())
3287 : {
3288 1 : CPLError(CE_Failure, CPLE_AppDefined,
3289 : "Did not get any cutline features.");
3290 1 : goto error;
3291 : }
3292 :
3293 : /* -------------------------------------------------------------------- */
3294 : /* Ensure the coordinate system gets set on the geometry. */
3295 : /* -------------------------------------------------------------------- */
3296 38 : if (!osSRS.empty())
3297 : {
3298 2 : auto poSRS = OGRSpatialReferenceRefCountedPtr::makeInstance();
3299 1 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3300 1 : poSRS->SetFromUserInput(osSRS.c_str());
3301 1 : poMultiPolygon->assignSpatialReference(poSRS.get());
3302 : }
3303 : else
3304 : {
3305 37 : poMultiPolygon->assignSpatialReference(poLayer->GetSpatialRef());
3306 : }
3307 :
3308 38 : *phCutlineRet = OGRGeometry::ToHandle(poMultiPolygon.release());
3309 :
3310 : /* -------------------------------------------------------------------- */
3311 : /* Cleanup */
3312 : /* -------------------------------------------------------------------- */
3313 38 : if (!osCSQL.empty())
3314 1 : poDS->ReleaseResultSet(poLayer);
3315 :
3316 38 : return CE_None;
3317 :
3318 4 : error:
3319 4 : if (!osCSQL.empty())
3320 1 : poDS->ReleaseResultSet(poLayer);
3321 :
3322 4 : return CE_Failure;
3323 : }
3324 :
3325 : /************************************************************************/
3326 : /* GDALWarpCreateOutput() */
3327 : /* */
3328 : /* Create the output file based on various command line options, */
3329 : /* and the input file. */
3330 : /* If there's just one source file, then hUniqueTransformArg will */
3331 : /* be set in order them to be reused by main function. This saves */
3332 : /* transform recomputation, which can be expensive in the -tps case*/
3333 : /************************************************************************/
3334 :
3335 1068 : static GDALDatasetH GDALWarpCreateOutput(
3336 : int nSrcCount, GDALDatasetH *pahSrcDS, const char *pszFilename,
3337 : const char *pszFormat, char **papszTO, CSLConstList papszCreateOptions,
3338 : GDALDataType eDT, GDALTransformerArgUniquePtr &hUniqueTransformArg,
3339 : bool bSetColorInterpretation, GDALWarpAppOptions *psOptions,
3340 : bool bUpdateTransformerWithDestGT)
3341 :
3342 : {
3343 : GDALDriverH hDriver;
3344 : GDALDatasetH hDstDS;
3345 1068 : GDALRasterAttributeTableH hRAT = nullptr;
3346 1068 : double dfWrkMinX = 0, dfWrkMaxX = 0, dfWrkMinY = 0, dfWrkMaxY = 0;
3347 1068 : double dfWrkResX = 0, dfWrkResY = 0;
3348 1068 : int nDstBandCount = 0;
3349 2136 : std::vector<GDALColorInterp> apeColorInterpretations;
3350 1068 : bool bVRT = false;
3351 :
3352 1068 : if (EQUAL(pszFormat, "VRT"))
3353 181 : bVRT = true;
3354 :
3355 : // Special case for geographic to Mercator (typically EPSG:4326 to EPSG:3857)
3356 : // where latitudes close to 90 go to infinity
3357 : // We clamp latitudes between ~ -85 and ~ 85 degrees.
3358 1068 : const char *pszDstSRS = CSLFetchNameValue(papszTO, "DST_SRS");
3359 1068 : if (nSrcCount == 1 && pszDstSRS && psOptions->dfMinX == 0.0 &&
3360 172 : psOptions->dfMinY == 0.0 && psOptions->dfMaxX == 0.0 &&
3361 172 : psOptions->dfMaxY == 0.0)
3362 : {
3363 172 : auto hSrcDS = pahSrcDS[0];
3364 344 : const auto osSrcSRS = GetSrcDSProjection(pahSrcDS[0], papszTO);
3365 344 : OGRSpatialReference oSrcSRS;
3366 172 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3367 172 : oSrcSRS.SetFromUserInput(osSrcSRS.c_str());
3368 344 : OGRSpatialReference oDstSRS;
3369 172 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3370 172 : oDstSRS.SetFromUserInput(pszDstSRS);
3371 172 : const char *pszProjection = oDstSRS.GetAttrValue("PROJECTION");
3372 172 : const char *pszMethod = FetchSrcMethod(papszTO);
3373 172 : double adfSrcGT[6] = {0};
3374 : // This MAX_LAT values is equivalent to the semi_major_axis * PI
3375 : // easting/northing value only for EPSG:3857, but it is also quite
3376 : // reasonable for other Mercator projections
3377 172 : constexpr double MAX_LAT = 85.0511287798066;
3378 172 : constexpr double EPS = 1e-3;
3379 5 : const auto GetMinLon = [&adfSrcGT]() { return adfSrcGT[0]; };
3380 5 : const auto GetMaxLon = [&adfSrcGT, hSrcDS]()
3381 5 : { return adfSrcGT[0] + adfSrcGT[1] * GDALGetRasterXSize(hSrcDS); };
3382 5 : const auto GetMinLat = [&adfSrcGT, hSrcDS]()
3383 5 : { return adfSrcGT[3] + adfSrcGT[5] * GDALGetRasterYSize(hSrcDS); };
3384 6 : const auto GetMaxLat = [&adfSrcGT]() { return adfSrcGT[3]; };
3385 238 : if (oSrcSRS.IsGeographic() && !oSrcSRS.IsDerivedGeographic() &&
3386 39 : pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) &&
3387 3 : oDstSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0) == 0 &&
3388 0 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
3389 3 : CSLFetchNameValue(papszTO, "COORDINATE_OPERATION") == nullptr &&
3390 3 : CSLFetchNameValue(papszTO, "SRC_METHOD") == nullptr &&
3391 3 : CSLFetchNameValue(papszTO, "DST_METHOD") == nullptr &&
3392 3 : GDALGetGeoTransform(hSrcDS, adfSrcGT) == CE_None &&
3393 6 : adfSrcGT[2] == 0 && adfSrcGT[4] == 0 && adfSrcGT[5] < 0 &&
3394 9 : GetMinLon() >= -180 - EPS && GetMaxLon() <= 180 + EPS &&
3395 6 : ((GetMaxLat() > MAX_LAT && GetMinLat() < MAX_LAT) ||
3396 2 : (GetMaxLat() > -MAX_LAT && GetMinLat() < -MAX_LAT)) &&
3397 240 : GDALGetMetadata(hSrcDS, "GEOLOC_ARRAY") == nullptr &&
3398 2 : GDALGetMetadata(hSrcDS, "RPC") == nullptr)
3399 : {
3400 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
3401 4 : OGRCreateCoordinateTransformation(&oSrcSRS, &oDstSRS));
3402 2 : if (poCT)
3403 : {
3404 2 : double xLL = std::max(GetMinLon(), -180.0);
3405 2 : double yLL = std::max(GetMinLat(), -MAX_LAT);
3406 2 : double xUR = std::min(GetMaxLon(), 180.0);
3407 2 : double yUR = std::min(GetMaxLat(), MAX_LAT);
3408 4 : if (poCT->Transform(1, &xLL, &yLL) &&
3409 2 : poCT->Transform(1, &xUR, &yUR))
3410 : {
3411 2 : psOptions->dfMinX = xLL;
3412 2 : psOptions->dfMinY = yLL;
3413 2 : psOptions->dfMaxX = xUR;
3414 2 : psOptions->dfMaxY = yUR;
3415 2 : CPLError(CE_Warning, CPLE_AppDefined,
3416 : "Clamping output bounds to (%f,%f) -> (%f, %f)",
3417 : psOptions->dfMinX, psOptions->dfMinY,
3418 : psOptions->dfMaxX, psOptions->dfMaxY);
3419 : }
3420 : }
3421 : }
3422 : }
3423 :
3424 : /* If (-ts and -te) or (-tr and -te) are specified, we don't need to compute
3425 : * the suggested output extent */
3426 1068 : const bool bNeedsSuggestedWarpOutput =
3427 2136 : !(((psOptions->nForcePixels != 0 && psOptions->nForceLines != 0) ||
3428 900 : (psOptions->dfXRes != 0 && psOptions->dfYRes != 0)) &&
3429 261 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
3430 133 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0));
3431 :
3432 : // If -te is specified, not not -tr and -ts
3433 1068 : const bool bKnownTargetExtentButNotResolution =
3434 868 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
3435 863 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0) &&
3436 207 : psOptions->nForcePixels == 0 && psOptions->nForceLines == 0 &&
3437 2136 : psOptions->dfXRes == 0 && psOptions->dfYRes == 0;
3438 :
3439 : /* -------------------------------------------------------------------- */
3440 : /* Find the output driver. */
3441 : /* -------------------------------------------------------------------- */
3442 1068 : hDriver = GDALGetDriverByName(pszFormat);
3443 2136 : if (hDriver == nullptr ||
3444 1068 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) == nullptr &&
3445 0 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) ==
3446 : nullptr))
3447 : {
3448 : auto poMissingDriver =
3449 0 : GetGDALDriverManager()->GetHiddenDriverByName(pszFormat);
3450 0 : if (poMissingDriver)
3451 : {
3452 : const std::string msg =
3453 0 : GDALGetMessageAboutMissingPluginDriver(poMissingDriver);
3454 0 : printf("Output driver `%s' not found but is known. However plugin "
3455 : "%s\n",
3456 : pszFormat, msg.c_str());
3457 0 : return nullptr;
3458 : }
3459 :
3460 0 : printf("Output driver `%s' not recognised or does not support\n",
3461 : pszFormat);
3462 0 : printf("direct output file creation or CreateCopy. "
3463 : "The following format drivers are eligible for warp output:\n");
3464 :
3465 0 : for (int iDr = 0; iDr < GDALGetDriverCount(); iDr++)
3466 : {
3467 0 : hDriver = GDALGetDriver(iDr);
3468 :
3469 0 : if (GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, nullptr) !=
3470 0 : nullptr &&
3471 0 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) !=
3472 0 : nullptr ||
3473 0 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) !=
3474 : nullptr))
3475 : {
3476 0 : printf(" %s: %s\n", GDALGetDriverShortName(hDriver),
3477 : GDALGetDriverLongName(hDriver));
3478 : }
3479 : }
3480 0 : printf("\n");
3481 0 : return nullptr;
3482 : }
3483 :
3484 : /* -------------------------------------------------------------------- */
3485 : /* For virtual output files, we have to set a special subclass */
3486 : /* of dataset to create. */
3487 : /* -------------------------------------------------------------------- */
3488 2136 : CPLStringList aosCreateOptions(CSLDuplicate(papszCreateOptions));
3489 1068 : if (bVRT)
3490 181 : aosCreateOptions.SetNameValue("SUBCLASS", "VRTWarpedDataset");
3491 :
3492 : /* -------------------------------------------------------------------- */
3493 : /* Loop over all input files to collect extents. */
3494 : /* -------------------------------------------------------------------- */
3495 2136 : CPLString osThisTargetSRS;
3496 : {
3497 1068 : const char *pszThisTargetSRS = CSLFetchNameValue(papszTO, "DST_SRS");
3498 1068 : if (pszThisTargetSRS != nullptr)
3499 286 : osThisTargetSRS = pszThisTargetSRS;
3500 : }
3501 :
3502 2136 : CPLStringList aoTOList(papszTO, FALSE);
3503 :
3504 1068 : double dfResFromSourceAndTargetExtent =
3505 : std::numeric_limits<double>::infinity();
3506 :
3507 : /* -------------------------------------------------------------------- */
3508 : /* Establish list of files of output dataset if it already exists. */
3509 : /* -------------------------------------------------------------------- */
3510 2136 : std::set<std::string> oSetExistingDestFiles;
3511 : {
3512 1068 : CPLPushErrorHandler(CPLQuietErrorHandler);
3513 1068 : const char *const apszAllowedDrivers[] = {pszFormat, nullptr};
3514 : auto poExistingOutputDS = std::unique_ptr<GDALDataset>(
3515 2136 : GDALDataset::Open(pszFilename, GDAL_OF_RASTER, apszAllowedDrivers));
3516 1068 : if (poExistingOutputDS)
3517 : {
3518 74 : for (const char *pszFilenameInList :
3519 70 : CPLStringList(poExistingOutputDS->GetFileList()))
3520 : {
3521 : oSetExistingDestFiles.insert(
3522 37 : CPLString(pszFilenameInList).replaceAll('\\', '/'));
3523 : }
3524 : }
3525 1068 : CPLPopErrorHandler();
3526 : }
3527 2136 : std::set<std::string> oSetExistingDestFilesFoundInSource;
3528 1068 : std::unique_ptr<GDALColorTable> poCT;
3529 :
3530 2150 : for (int iSrc = 0; iSrc < nSrcCount; iSrc++)
3531 : {
3532 : /* --------------------------------------------------------------------
3533 : */
3534 : /* Check that there's at least one raster band */
3535 : /* --------------------------------------------------------------------
3536 : */
3537 1092 : GDALDatasetH hSrcDS = pahSrcDS[iSrc];
3538 1092 : if (GDALGetRasterCount(hSrcDS) == 0)
3539 : {
3540 1 : CPLError(CE_Failure, CPLE_AppDefined,
3541 : "Input file %s has no raster bands.",
3542 : GDALGetDescription(hSrcDS));
3543 10 : return nullptr;
3544 : }
3545 :
3546 : // Examine desired overview level and retrieve the corresponding dataset
3547 : // if it exists.
3548 0 : std::unique_ptr<GDALDataset> oDstDSOverview;
3549 1091 : if (psOptions->nOvLevel >= 0)
3550 : {
3551 5 : oDstDSOverview.reset(GDALCreateOverviewDataset(
3552 : GDALDataset::FromHandle(hSrcDS), psOptions->nOvLevel,
3553 : /* bThisLevelOnly = */ true));
3554 5 : if (oDstDSOverview)
3555 4 : hSrcDS = oDstDSOverview.get();
3556 : }
3557 :
3558 : /* --------------------------------------------------------------------
3559 : */
3560 : /* Check if the source dataset shares some files with the dest
3561 : * one.*/
3562 : /* --------------------------------------------------------------------
3563 : */
3564 1091 : if (!oSetExistingDestFiles.empty())
3565 : {
3566 : // We need to reopen in a temporary dataset for the particular
3567 : // case of overwritten a .tif.ovr file from a .tif
3568 : // If we probe the file list of the .tif, it will then open the
3569 : // .tif.ovr !
3570 36 : auto poSrcDS = GDALDataset::FromHandle(hSrcDS);
3571 36 : const char *const apszAllowedDrivers[] = {
3572 36 : poSrcDS->GetDriver() ? poSrcDS->GetDriver()->GetDescription()
3573 : : nullptr,
3574 36 : nullptr};
3575 : auto poSrcDSTmp = std::unique_ptr<GDALDataset>(GDALDataset::Open(
3576 72 : poSrcDS->GetDescription(), GDAL_OF_RASTER, apszAllowedDrivers));
3577 36 : if (poSrcDSTmp)
3578 : {
3579 38 : for (const char *pszFilenameInList :
3580 74 : CPLStringList(poSrcDSTmp->GetFileList()))
3581 : {
3582 : std::string osFilename =
3583 76 : CPLString(pszFilenameInList).replaceAll('\\', '/');
3584 38 : if (oSetExistingDestFiles.find(osFilename) !=
3585 76 : oSetExistingDestFiles.end())
3586 : {
3587 : oSetExistingDestFilesFoundInSource.insert(
3588 5 : std::move(osFilename));
3589 : }
3590 : }
3591 : }
3592 : }
3593 :
3594 1091 : if (eDT == GDT_Unknown)
3595 1030 : eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));
3596 :
3597 : /* --------------------------------------------------------------------
3598 : */
3599 : /* If we are processing the first file, and it has a raster */
3600 : /* attribute table, then we will copy it to the destination file.
3601 : */
3602 : /* --------------------------------------------------------------------
3603 : */
3604 1091 : if (iSrc == 0)
3605 : {
3606 1065 : hRAT = GDALGetDefaultRAT(GDALGetRasterBand(hSrcDS, 1));
3607 1065 : if (hRAT != nullptr)
3608 : {
3609 0 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
3610 0 : psOptions->eResampleAlg != GRA_Mode &&
3611 0 : GDALRATGetTableType(hRAT) == GRTT_THEMATIC)
3612 : {
3613 0 : if (!psOptions->bQuiet)
3614 : {
3615 0 : CPLError(CE_Warning, CPLE_AppDefined,
3616 : "Warning: Input file %s has a thematic RAT, "
3617 : "which will likely lead "
3618 : "to bad results when using a resampling "
3619 : "method other than nearest neighbour "
3620 : "or mode so we are discarding it.\n",
3621 : GDALGetDescription(hSrcDS));
3622 : }
3623 0 : hRAT = nullptr;
3624 : }
3625 : else
3626 : {
3627 0 : if (!psOptions->bQuiet)
3628 0 : printf("Copying raster attribute table from %s to new "
3629 : "file.\n",
3630 : GDALGetDescription(hSrcDS));
3631 : }
3632 : }
3633 : }
3634 :
3635 : /* --------------------------------------------------------------------
3636 : */
3637 : /* If we are processing the first file, and it has a color */
3638 : /* table, then we will copy it to the destination file. */
3639 : /* --------------------------------------------------------------------
3640 : */
3641 1091 : if (iSrc == 0)
3642 : {
3643 1065 : auto hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1));
3644 1065 : if (hCT != nullptr)
3645 : {
3646 9 : poCT.reset(
3647 : GDALColorTable::FromHandle(GDALCloneColorTable(hCT)));
3648 9 : if (!psOptions->bQuiet)
3649 0 : printf("Copying color table from %s to new file.\n",
3650 : GDALGetDescription(hSrcDS));
3651 : }
3652 :
3653 1065 : if (psOptions->anDstBands.empty())
3654 : {
3655 1039 : nDstBandCount = GDALGetRasterCount(hSrcDS);
3656 2474 : for (int iBand = 0; iBand < nDstBandCount; iBand++)
3657 : {
3658 1435 : GDALColorInterp eInterp = GDALGetRasterColorInterpretation(
3659 1435 : GDALGetRasterBand(hSrcDS, iBand + 1));
3660 1435 : apeColorInterpretations.push_back(eInterp);
3661 : }
3662 :
3663 : // Do we want to generate an alpha band in the output file?
3664 1039 : if (psOptions->bEnableSrcAlpha)
3665 31 : nDstBandCount--;
3666 :
3667 1039 : if (psOptions->bEnableDstAlpha)
3668 90 : nDstBandCount++;
3669 : }
3670 : else
3671 : {
3672 68 : for (int nSrcBand : psOptions->anSrcBands)
3673 : {
3674 42 : auto hBand = GDALGetRasterBand(hSrcDS, nSrcBand);
3675 : GDALColorInterp eInterp =
3676 42 : hBand ? GDALGetRasterColorInterpretation(hBand)
3677 42 : : GCI_Undefined;
3678 42 : apeColorInterpretations.push_back(eInterp);
3679 : }
3680 26 : nDstBandCount = static_cast<int>(psOptions->anDstBands.size());
3681 26 : if (psOptions->bEnableDstAlpha)
3682 : {
3683 9 : nDstBandCount++;
3684 9 : apeColorInterpretations.push_back(GCI_AlphaBand);
3685 : }
3686 17 : else if (GDALGetRasterCount(hSrcDS) &&
3687 17 : GDALGetRasterColorInterpretation(GDALGetRasterBand(
3688 : hSrcDS, GDALGetRasterCount(hSrcDS))) ==
3689 34 : GCI_AlphaBand &&
3690 1 : !psOptions->bDisableSrcAlpha)
3691 : {
3692 0 : nDstBandCount++;
3693 0 : apeColorInterpretations.push_back(GCI_AlphaBand);
3694 : }
3695 : }
3696 : }
3697 :
3698 : /* --------------------------------------------------------------------
3699 : */
3700 : /* If we are processing the first file, get the source srs from */
3701 : /* dataset, if not set already. */
3702 : /* --------------------------------------------------------------------
3703 : */
3704 1091 : const auto osThisSourceSRS = GetSrcDSProjection(hSrcDS, papszTO);
3705 1091 : if (iSrc == 0 && osThisTargetSRS.empty())
3706 : {
3707 779 : if (!osThisSourceSRS.empty())
3708 : {
3709 577 : osThisTargetSRS = osThisSourceSRS;
3710 577 : aoTOList.SetNameValue("DST_SRS", osThisSourceSRS);
3711 : }
3712 : }
3713 :
3714 : /* --------------------------------------------------------------------
3715 : */
3716 : /* Create a transformation object from the source to */
3717 : /* destination coordinate system. */
3718 : /* --------------------------------------------------------------------
3719 : */
3720 0 : GDALTransformerArgUniquePtr hTransformArg;
3721 1091 : if (hUniqueTransformArg)
3722 2 : hTransformArg = std::move(hUniqueTransformArg);
3723 : else
3724 : {
3725 1089 : hTransformArg.reset(GDALCreateGenImgProjTransformer2(
3726 1089 : hSrcDS, nullptr, aoTOList.List()));
3727 1089 : if (hTransformArg == nullptr)
3728 : {
3729 9 : return nullptr;
3730 : }
3731 : }
3732 :
3733 : GDALTransformerInfo *psInfo =
3734 1082 : static_cast<GDALTransformerInfo *>(hTransformArg.get());
3735 :
3736 : /* --------------------------------------------------------------------
3737 : */
3738 : /* Get approximate output resolution */
3739 : /* --------------------------------------------------------------------
3740 : */
3741 :
3742 1082 : if (bKnownTargetExtentButNotResolution)
3743 : {
3744 : // Sample points along a grid in target CRS
3745 84 : constexpr int nPointsX = 10;
3746 84 : constexpr int nPointsY = 10;
3747 84 : constexpr int nPoints = 3 * nPointsX * nPointsY;
3748 168 : std::vector<double> padfX;
3749 168 : std::vector<double> padfY;
3750 168 : std::vector<double> padfZ(nPoints);
3751 168 : std::vector<int> pabSuccess(nPoints);
3752 : const double dfEps =
3753 168 : std::min(psOptions->dfMaxX - psOptions->dfMinX,
3754 84 : std::abs(psOptions->dfMaxY - psOptions->dfMinY)) /
3755 84 : 1000;
3756 924 : for (int iY = 0; iY < nPointsY; iY++)
3757 : {
3758 9240 : for (int iX = 0; iX < nPointsX; iX++)
3759 : {
3760 8400 : const double dfX =
3761 8400 : psOptions->dfMinX +
3762 8400 : static_cast<double>(iX) *
3763 8400 : (psOptions->dfMaxX - psOptions->dfMinX) /
3764 : (nPointsX - 1);
3765 8400 : const double dfY =
3766 8400 : psOptions->dfMinY +
3767 8400 : static_cast<double>(iY) *
3768 8400 : (psOptions->dfMaxY - psOptions->dfMinY) /
3769 : (nPointsY - 1);
3770 :
3771 : // Reproject each destination sample point and its
3772 : // neighbours at (x+1,y) and (x,y+1), so as to get the local
3773 : // scale.
3774 8400 : padfX.push_back(dfX);
3775 8400 : padfY.push_back(dfY);
3776 :
3777 15960 : padfX.push_back((iX == nPointsX - 1) ? dfX - dfEps
3778 7560 : : dfX + dfEps);
3779 8400 : padfY.push_back(dfY);
3780 :
3781 8400 : padfX.push_back(dfX);
3782 15960 : padfY.push_back((iY == nPointsY - 1) ? dfY - dfEps
3783 7560 : : dfY + dfEps);
3784 : }
3785 : }
3786 :
3787 84 : bool transformedToSrcCRS{false};
3788 :
3789 : GDALGenImgProjTransformInfo *psTransformInfo{
3790 : static_cast<GDALGenImgProjTransformInfo *>(
3791 84 : hTransformArg.get())};
3792 :
3793 : // If a transformer is available, use an extent that covers the
3794 : // target extent instead of the real source image extent, but also
3795 : // check for target extent compatibility with source CRS extent
3796 84 : if (psTransformInfo && psTransformInfo->pReprojectArg &&
3797 64 : psTransformInfo->sSrcParams.pTransformer == nullptr)
3798 : {
3799 64 : const GDALReprojectionTransformInfo *psRTI =
3800 : static_cast<const GDALReprojectionTransformInfo *>(
3801 : psTransformInfo->pReprojectArg);
3802 64 : if (psRTI->poReverseTransform)
3803 : {
3804 :
3805 : // Compute new geotransform from transformed target extent
3806 : double adfGeoTransform[6];
3807 64 : if (GDALGetGeoTransform(hSrcDS, adfGeoTransform) ==
3808 64 : CE_None &&
3809 64 : adfGeoTransform[2] == 0 && adfGeoTransform[4] == 0)
3810 : {
3811 :
3812 : // Transform target extent to source CRS
3813 64 : double dfMinX = psOptions->dfMinX;
3814 64 : double dfMinY = psOptions->dfMinY;
3815 :
3816 : // Need this to check if the target extent is compatible with the source extent
3817 64 : double dfMaxX = psOptions->dfMaxX;
3818 64 : double dfMaxY = psOptions->dfMaxY;
3819 :
3820 : // Clone of psRTI->poReverseTransform with CHECK_WITH_INVERT_PROJ set to TRUE
3821 : // to detect out of source CRS bounds destination extent and fall back to original
3822 : // algorithm if needed
3823 : CPLConfigOptionSetter oSetter("CHECK_WITH_INVERT_PROJ",
3824 128 : "TRUE", false);
3825 128 : OGRCoordinateTransformationOptions options;
3826 : auto poReverseTransform =
3827 : std::unique_ptr<OGRCoordinateTransformation>(
3828 : OGRCreateCoordinateTransformation(
3829 64 : psRTI->poReverseTransform->GetSourceCS(),
3830 64 : psRTI->poReverseTransform->GetTargetCS(),
3831 128 : options));
3832 :
3833 64 : if (poReverseTransform)
3834 : {
3835 :
3836 128 : poReverseTransform->Transform(
3837 64 : 1, &dfMinX, &dfMinY, nullptr, &pabSuccess[0]);
3838 :
3839 64 : if (pabSuccess[0])
3840 : {
3841 60 : adfGeoTransform[0] = dfMinX;
3842 60 : adfGeoTransform[3] = dfMinY;
3843 :
3844 120 : poReverseTransform->Transform(1, &dfMaxX,
3845 : &dfMaxY, nullptr,
3846 60 : &pabSuccess[0]);
3847 :
3848 60 : if (pabSuccess[0])
3849 : {
3850 :
3851 : // Reproject to source image CRS
3852 59 : psRTI->poReverseTransform->Transform(
3853 59 : nPoints, &padfX[0], &padfY[0],
3854 59 : &padfZ[0], &pabSuccess[0]);
3855 :
3856 : // Transform back to source image coordinate space using geotransform
3857 17759 : for (size_t i = 0; i < padfX.size(); i++)
3858 : {
3859 35400 : padfX[i] =
3860 17700 : (padfX[i] - adfGeoTransform[0]) /
3861 17700 : adfGeoTransform[1];
3862 17700 : padfY[i] = std::abs(
3863 17700 : (padfY[i] - adfGeoTransform[3]) /
3864 17700 : adfGeoTransform[5]);
3865 : }
3866 :
3867 59 : transformedToSrcCRS = true;
3868 : }
3869 : }
3870 : }
3871 : }
3872 : }
3873 : }
3874 :
3875 84 : if (!transformedToSrcCRS)
3876 : {
3877 : // Transform to source image coordinate space
3878 25 : psInfo->pfnTransform(hTransformArg.get(), TRUE, nPoints,
3879 25 : &padfX[0], &padfY[0], &padfZ[0],
3880 25 : &pabSuccess[0]);
3881 : }
3882 :
3883 : // Compute the resolution at sampling points
3884 168 : std::vector<std::pair<double, double>> aoResPairs;
3885 :
3886 15062 : const auto Distance = [](double x, double y)
3887 15062 : { return sqrt(x * x + y * y); };
3888 :
3889 84 : const int nSrcXSize = GDALGetRasterXSize(hSrcDS);
3890 84 : const int nSrcYSize = GDALGetRasterYSize(hSrcDS);
3891 :
3892 8484 : for (int i = 0; i < nPoints; i += 3)
3893 : {
3894 16800 : if (pabSuccess[i] && pabSuccess[i + 1] && pabSuccess[i + 2] &&
3895 18925 : padfX[i] >= 0 && padfY[i] >= 0 &&
3896 2125 : (transformedToSrcCRS ||
3897 2125 : (padfX[i] <= nSrcXSize && padfY[i] <= nSrcYSize)))
3898 : {
3899 : const double dfRes1 =
3900 15062 : std::abs(dfEps) / Distance(padfX[i + 1] - padfX[i],
3901 7531 : padfY[i + 1] - padfY[i]);
3902 : const double dfRes2 =
3903 15062 : std::abs(dfEps) / Distance(padfX[i + 2] - padfX[i],
3904 7531 : padfY[i + 2] - padfY[i]);
3905 7531 : if (std::isfinite(dfRes1) && std::isfinite(dfRes2))
3906 : {
3907 7511 : aoResPairs.push_back(
3908 15022 : std::pair<double, double>(dfRes1, dfRes2));
3909 : }
3910 : }
3911 : }
3912 :
3913 : // Find the minimum resolution that is at least 10 times greater
3914 : // than the median, to remove outliers.
3915 : // Start first by doing that on dfRes1, then dfRes2 and then their
3916 : // average.
3917 84 : std::sort(aoResPairs.begin(), aoResPairs.end(),
3918 38467 : [](const std::pair<double, double> &oPair1,
3919 : const std::pair<double, double> &oPair2)
3920 38467 : { return oPair1.first < oPair2.first; });
3921 :
3922 84 : if (!aoResPairs.empty())
3923 : {
3924 168 : std::vector<std::pair<double, double>> aoResPairsNew;
3925 : const double dfMedian1 =
3926 84 : aoResPairs[aoResPairs.size() / 2].first;
3927 7595 : for (const auto &oPair : aoResPairs)
3928 : {
3929 7511 : if (oPair.first > dfMedian1 / 10)
3930 : {
3931 7481 : aoResPairsNew.push_back(oPair);
3932 : }
3933 : }
3934 :
3935 84 : aoResPairs = std::move(aoResPairsNew);
3936 84 : std::sort(aoResPairs.begin(), aoResPairs.end(),
3937 40580 : [](const std::pair<double, double> &oPair1,
3938 : const std::pair<double, double> &oPair2)
3939 40580 : { return oPair1.second < oPair2.second; });
3940 84 : if (!aoResPairs.empty())
3941 : {
3942 168 : std::vector<double> adfRes;
3943 : const double dfMedian2 =
3944 84 : aoResPairs[aoResPairs.size() / 2].second;
3945 7565 : for (const auto &oPair : aoResPairs)
3946 : {
3947 7481 : if (oPair.second > dfMedian2 / 10)
3948 : {
3949 7471 : adfRes.push_back((oPair.first + oPair.second) / 2);
3950 : }
3951 : }
3952 :
3953 84 : std::sort(adfRes.begin(), adfRes.end());
3954 84 : if (!adfRes.empty())
3955 : {
3956 84 : const double dfMedian = adfRes[adfRes.size() / 2];
3957 84 : for (const double dfRes : adfRes)
3958 : {
3959 84 : if (dfRes > dfMedian / 10)
3960 : {
3961 84 : dfResFromSourceAndTargetExtent = std::min(
3962 84 : dfResFromSourceAndTargetExtent, dfRes);
3963 84 : break;
3964 : }
3965 : }
3966 : }
3967 : }
3968 : }
3969 : }
3970 :
3971 : /* --------------------------------------------------------------------
3972 : */
3973 : /* Get approximate output definition. */
3974 : /* --------------------------------------------------------------------
3975 : */
3976 : double adfThisGeoTransform[6];
3977 : double adfExtent[4];
3978 1082 : if (bNeedsSuggestedWarpOutput)
3979 : {
3980 : int nThisPixels, nThisLines;
3981 :
3982 : // For sum, round-up dimension, to be sure that the output extent
3983 : // includes all source pixels, to have the sum preserving property.
3984 1904 : int nOptions = (psOptions->eResampleAlg == GRA_Sum)
3985 952 : ? GDAL_SWO_ROUND_UP_SIZE
3986 : : 0;
3987 952 : if (psOptions->bSquarePixels)
3988 : {
3989 1 : nOptions |= GDAL_SWO_FORCE_SQUARE_PIXEL;
3990 : }
3991 :
3992 952 : if (GDALSuggestedWarpOutput2(
3993 : hSrcDS, psInfo->pfnTransform, hTransformArg.get(),
3994 : adfThisGeoTransform, &nThisPixels, &nThisLines, adfExtent,
3995 952 : nOptions) != CE_None)
3996 : {
3997 0 : return nullptr;
3998 : }
3999 :
4000 952 : if (CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", nullptr) ==
4001 : nullptr)
4002 : {
4003 952 : double MinX = adfExtent[0];
4004 952 : double MaxX = adfExtent[2];
4005 952 : double MaxY = adfExtent[3];
4006 952 : double MinY = adfExtent[1];
4007 952 : int bSuccess = TRUE;
4008 :
4009 : // +/-180 deg in longitude do not roundtrip sometimes
4010 952 : if (MinX == -180)
4011 51 : MinX += 1e-6;
4012 952 : if (MaxX == 180)
4013 50 : MaxX -= 1e-6;
4014 :
4015 : // +/-90 deg in latitude do not roundtrip sometimes
4016 952 : if (MinY == -90)
4017 52 : MinY += 1e-6;
4018 952 : if (MaxY == 90)
4019 54 : MaxY -= 1e-6;
4020 :
4021 : /* Check that the edges of the target image are in the validity
4022 : * area */
4023 : /* of the target projection */
4024 952 : const int N_STEPS = 20;
4025 20339 : for (int i = 0; i <= N_STEPS && bSuccess; i++)
4026 : {
4027 425998 : for (int j = 0; j <= N_STEPS && bSuccess; j++)
4028 : {
4029 406611 : const double dfRatioI = i * 1.0 / N_STEPS;
4030 406611 : const double dfRatioJ = j * 1.0 / N_STEPS;
4031 406611 : const double expected_x =
4032 406611 : (1 - dfRatioI) * MinX + dfRatioI * MaxX;
4033 406611 : const double expected_y =
4034 406611 : (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
4035 406611 : double x = expected_x;
4036 406611 : double y = expected_y;
4037 406611 : double z = 0;
4038 : /* Target SRS coordinates to source image pixel
4039 : * coordinates */
4040 406611 : if (!psInfo->pfnTransform(hTransformArg.get(), TRUE, 1,
4041 813202 : &x, &y, &z, &bSuccess) ||
4042 406591 : !bSuccess)
4043 20 : bSuccess = FALSE;
4044 : /* Source image pixel coordinates to target SRS
4045 : * coordinates */
4046 406611 : if (!psInfo->pfnTransform(hTransformArg.get(), FALSE, 1,
4047 813200 : &x, &y, &z, &bSuccess) ||
4048 406589 : !bSuccess)
4049 22 : bSuccess = FALSE;
4050 406611 : if (fabs(x - expected_x) >
4051 406611 : (MaxX - MinX) / nThisPixels ||
4052 406580 : fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
4053 31 : bSuccess = FALSE;
4054 : }
4055 : }
4056 :
4057 : /* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces
4058 : * ogrct.cpp */
4059 : /* to check the consistency of each requested projection result
4060 : * with the */
4061 : /* invert projection */
4062 952 : if (!bSuccess)
4063 : {
4064 31 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ",
4065 : "TRUE");
4066 31 : CPLDebug("WARP", "Recompute out extent with "
4067 : "CHECK_WITH_INVERT_PROJ=TRUE");
4068 :
4069 31 : const CPLErr eErr = GDALSuggestedWarpOutput2(
4070 : hSrcDS, psInfo->pfnTransform, hTransformArg.get(),
4071 : adfThisGeoTransform, &nThisPixels, &nThisLines,
4072 : adfExtent, 0);
4073 31 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ",
4074 : nullptr);
4075 31 : if (eErr != CE_None)
4076 : {
4077 0 : return nullptr;
4078 : }
4079 : }
4080 : }
4081 : }
4082 :
4083 : // If no reprojection or geometry change is involved, and that the
4084 : // source image is north-up, preserve source resolution instead of
4085 : // forcing square pixels.
4086 1082 : const char *pszMethod = FetchSrcMethod(papszTO);
4087 : double adfThisGeoTransformTmp[6];
4088 1081 : if (!psOptions->bSquarePixels && bNeedsSuggestedWarpOutput &&
4089 951 : psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
4090 923 : psOptions->nForcePixels == 0 && psOptions->nForceLines == 0 &&
4091 32 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
4092 780 : CSLFetchNameValue(papszTO, "COORDINATE_OPERATION") == nullptr &&
4093 776 : CSLFetchNameValue(papszTO, "SRC_METHOD") == nullptr &&
4094 776 : CSLFetchNameValue(papszTO, "DST_METHOD") == nullptr &&
4095 773 : GDALGetGeoTransform(hSrcDS, adfThisGeoTransformTmp) == CE_None &&
4096 753 : adfThisGeoTransformTmp[2] == 0 && adfThisGeoTransformTmp[4] == 0 &&
4097 751 : adfThisGeoTransformTmp[5] < 0 &&
4098 2803 : GDALGetMetadata(hSrcDS, "GEOLOC_ARRAY") == nullptr &&
4099 640 : GDALGetMetadata(hSrcDS, "RPC") == nullptr)
4100 : {
4101 640 : bool bIsSameHorizontal = osThisSourceSRS == osThisTargetSRS;
4102 640 : if (!bIsSameHorizontal)
4103 : {
4104 328 : OGRSpatialReference oSrcSRS;
4105 328 : OGRSpatialReference oDstSRS;
4106 328 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
4107 164 : if (oSrcSRS.SetFromUserInput(osThisSourceSRS.c_str()) ==
4108 164 : OGRERR_NONE &&
4109 164 : oDstSRS.SetFromUserInput(osThisTargetSRS.c_str()) ==
4110 164 : OGRERR_NONE &&
4111 164 : (oSrcSRS.GetAxesCount() == 3 ||
4112 146 : oDstSRS.GetAxesCount() == 3) &&
4113 350 : oSrcSRS.DemoteTo2D(nullptr) == OGRERR_NONE &&
4114 22 : oDstSRS.DemoteTo2D(nullptr) == OGRERR_NONE)
4115 : {
4116 22 : bIsSameHorizontal = CPL_TO_BOOL(oSrcSRS.IsSame(&oDstSRS));
4117 : }
4118 : }
4119 640 : if (bIsSameHorizontal)
4120 : {
4121 482 : memcpy(adfThisGeoTransform, adfThisGeoTransformTmp,
4122 : 6 * sizeof(double));
4123 482 : adfExtent[0] = adfThisGeoTransform[0];
4124 482 : adfExtent[1] =
4125 964 : adfThisGeoTransform[3] +
4126 482 : GDALGetRasterYSize(hSrcDS) * adfThisGeoTransform[5];
4127 482 : adfExtent[2] =
4128 964 : adfThisGeoTransform[0] +
4129 482 : GDALGetRasterXSize(hSrcDS) * adfThisGeoTransform[1];
4130 482 : adfExtent[3] = adfThisGeoTransform[3];
4131 482 : dfResFromSourceAndTargetExtent =
4132 482 : std::numeric_limits<double>::infinity();
4133 : }
4134 : }
4135 :
4136 1082 : if (bNeedsSuggestedWarpOutput)
4137 : {
4138 : /* --------------------------------------------------------------------
4139 : */
4140 : /* Expand the working bounds to include this region, ensure the
4141 : */
4142 : /* working resolution is no more than this resolution. */
4143 : /* --------------------------------------------------------------------
4144 : */
4145 952 : if (dfWrkMaxX == 0.0 && dfWrkMinX == 0.0)
4146 : {
4147 928 : dfWrkMinX = adfExtent[0];
4148 928 : dfWrkMaxX = adfExtent[2];
4149 928 : dfWrkMaxY = adfExtent[3];
4150 928 : dfWrkMinY = adfExtent[1];
4151 928 : dfWrkResX = adfThisGeoTransform[1];
4152 928 : dfWrkResY = std::abs(adfThisGeoTransform[5]);
4153 : }
4154 : else
4155 : {
4156 24 : dfWrkMinX = std::min(dfWrkMinX, adfExtent[0]);
4157 24 : dfWrkMaxX = std::max(dfWrkMaxX, adfExtent[2]);
4158 24 : dfWrkMaxY = std::max(dfWrkMaxY, adfExtent[3]);
4159 24 : dfWrkMinY = std::min(dfWrkMinY, adfExtent[1]);
4160 24 : dfWrkResX = std::min(dfWrkResX, adfThisGeoTransform[1]);
4161 24 : dfWrkResY =
4162 24 : std::min(dfWrkResY, std::abs(adfThisGeoTransform[5]));
4163 : }
4164 : }
4165 :
4166 1082 : if (nSrcCount == 1)
4167 : {
4168 1030 : hUniqueTransformArg = std::move(hTransformArg);
4169 : }
4170 : }
4171 :
4172 : // If the source file(s) and the dest one share some files in common,
4173 : // only remove the files that are *not* in common
4174 1058 : if (!oSetExistingDestFilesFoundInSource.empty())
4175 : {
4176 14 : for (const std::string &osFilename : oSetExistingDestFiles)
4177 : {
4178 9 : if (oSetExistingDestFilesFoundInSource.find(osFilename) ==
4179 18 : oSetExistingDestFilesFoundInSource.end())
4180 : {
4181 4 : VSIUnlink(osFilename.c_str());
4182 : }
4183 : }
4184 : }
4185 :
4186 1058 : if (std::isfinite(dfResFromSourceAndTargetExtent))
4187 : {
4188 26 : dfWrkResX = dfResFromSourceAndTargetExtent;
4189 26 : dfWrkResY = dfResFromSourceAndTargetExtent;
4190 : }
4191 :
4192 : /* -------------------------------------------------------------------- */
4193 : /* Did we have any usable sources? */
4194 : /* -------------------------------------------------------------------- */
4195 1058 : if (nDstBandCount == 0)
4196 : {
4197 3 : CPLError(CE_Failure, CPLE_AppDefined, "No usable source images.");
4198 3 : return nullptr;
4199 : }
4200 :
4201 : /* -------------------------------------------------------------------- */
4202 : /* Turn the suggested region into a geotransform and suggested */
4203 : /* number of pixels and lines. */
4204 : /* -------------------------------------------------------------------- */
4205 1055 : double adfDstGeoTransform[6] = {0, 0, 0, 0, 0, 0};
4206 1055 : int nPixels = 0;
4207 1055 : int nLines = 0;
4208 :
4209 200 : const auto ComputePixelsFromResAndExtent = [psOptions]()
4210 : {
4211 400 : return std::max(1.0,
4212 400 : std::round((psOptions->dfMaxX - psOptions->dfMinX) /
4213 200 : psOptions->dfXRes));
4214 1055 : };
4215 :
4216 202 : const auto ComputeLinesFromResAndExtent = [psOptions]()
4217 : {
4218 : return std::max(
4219 404 : 1.0, std::round(std::fabs(psOptions->dfMaxY - psOptions->dfMinY) /
4220 202 : psOptions->dfYRes));
4221 1055 : };
4222 :
4223 1055 : if (bNeedsSuggestedWarpOutput)
4224 : {
4225 926 : adfDstGeoTransform[0] = dfWrkMinX;
4226 926 : adfDstGeoTransform[1] = dfWrkResX;
4227 926 : adfDstGeoTransform[2] = 0.0;
4228 926 : adfDstGeoTransform[3] = dfWrkMaxY;
4229 926 : adfDstGeoTransform[4] = 0.0;
4230 926 : adfDstGeoTransform[5] = -1 * dfWrkResY;
4231 :
4232 926 : const double dfPixels = (dfWrkMaxX - dfWrkMinX) / dfWrkResX;
4233 926 : const double dfLines = (dfWrkMaxY - dfWrkMinY) / dfWrkResY;
4234 : // guaranteed by GDALSuggestedWarpOutput2() behavior
4235 926 : CPLAssert(std::round(dfPixels) <= INT_MAX);
4236 926 : CPLAssert(std::round(dfLines) <= INT_MAX);
4237 926 : nPixels =
4238 926 : static_cast<int>(std::min<double>(std::round(dfPixels), INT_MAX));
4239 926 : nLines =
4240 926 : static_cast<int>(std::min<double>(std::round(dfLines), INT_MAX));
4241 : }
4242 :
4243 : /* -------------------------------------------------------------------- */
4244 : /* Did the user override some parameters? */
4245 : /* -------------------------------------------------------------------- */
4246 1055 : if (UseTEAndTSAndTRConsistently(psOptions))
4247 : {
4248 26 : adfDstGeoTransform[0] = psOptions->dfMinX;
4249 26 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4250 26 : adfDstGeoTransform[1] = psOptions->dfXRes;
4251 26 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4252 :
4253 26 : nPixels = psOptions->nForcePixels;
4254 26 : nLines = psOptions->nForceLines;
4255 : }
4256 1029 : else if (psOptions->dfXRes != 0.0 && psOptions->dfYRes != 0.0)
4257 : {
4258 93 : bool bDetectBlankBorders = false;
4259 :
4260 93 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4261 29 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4262 : {
4263 28 : bDetectBlankBorders = bNeedsSuggestedWarpOutput;
4264 :
4265 28 : psOptions->dfMinX = adfDstGeoTransform[0];
4266 28 : psOptions->dfMaxX =
4267 28 : adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
4268 28 : psOptions->dfMaxY = adfDstGeoTransform[3];
4269 28 : psOptions->dfMinY =
4270 28 : adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
4271 : }
4272 :
4273 182 : if (psOptions->bTargetAlignedPixels ||
4274 89 : (psOptions->bCropToCutline &&
4275 1 : psOptions->aosWarpOptions.FetchBool("CUTLINE_ALL_TOUCHED", false)))
4276 : {
4277 4 : if ((psOptions->bTargetAlignedPixels &&
4278 10 : bNeedsSuggestedWarpOutput) ||
4279 2 : (psOptions->bCropToCutline &&
4280 1 : psOptions->aosWarpOptions.FetchBool("CUTLINE_ALL_TOUCHED",
4281 : false)))
4282 : {
4283 4 : bDetectBlankBorders = true;
4284 : }
4285 5 : constexpr double EPS = 1e-8;
4286 5 : psOptions->dfMinX =
4287 5 : floor(psOptions->dfMinX / psOptions->dfXRes + EPS) *
4288 5 : psOptions->dfXRes;
4289 5 : psOptions->dfMaxX =
4290 5 : ceil(psOptions->dfMaxX / psOptions->dfXRes - EPS) *
4291 5 : psOptions->dfXRes;
4292 5 : psOptions->dfMinY =
4293 5 : floor(psOptions->dfMinY / psOptions->dfYRes + EPS) *
4294 5 : psOptions->dfYRes;
4295 5 : psOptions->dfMaxY =
4296 5 : ceil(psOptions->dfMaxY / psOptions->dfYRes - EPS) *
4297 5 : psOptions->dfYRes;
4298 : }
4299 :
4300 120 : const auto UpdateGeoTransformandAndPixelLines = [&]()
4301 : {
4302 120 : const double dfPixels = ComputePixelsFromResAndExtent();
4303 120 : const double dfLines = ComputeLinesFromResAndExtent();
4304 120 : if (dfPixels > INT_MAX || dfLines > INT_MAX)
4305 : {
4306 2 : CPLError(CE_Failure, CPLE_AppDefined,
4307 : "Too large output raster size: %f x %f", dfPixels,
4308 : dfLines);
4309 2 : return false;
4310 : }
4311 :
4312 118 : nPixels = static_cast<int>(dfPixels);
4313 118 : nLines = static_cast<int>(dfLines);
4314 354 : adfDstGeoTransform[0] = psOptions->dfMinX;
4315 118 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4316 118 : adfDstGeoTransform[1] = psOptions->dfXRes;
4317 236 : adfDstGeoTransform[5] = (psOptions->dfMaxY > psOptions->dfMinY)
4318 118 : ? -psOptions->dfYRes
4319 2 : : psOptions->dfYRes;
4320 118 : return true;
4321 93 : };
4322 :
4323 122 : if (bDetectBlankBorders && nSrcCount == 1 && hUniqueTransformArg &&
4324 : // to avoid too large memory allocations
4325 29 : std::max(nPixels, nLines) < 100 * 1000 * 1000)
4326 : {
4327 : // Try to detect if the edge of the raster would be blank
4328 : // Cf https://github.com/OSGeo/gdal/issues/7905
4329 30 : while (nPixels > 1 || nLines > 1)
4330 : {
4331 27 : if (!UpdateGeoTransformandAndPixelLines())
4332 0 : return nullptr;
4333 :
4334 27 : GDALSetGenImgProjTransformerDstGeoTransform(
4335 : hUniqueTransformArg.get(), adfDstGeoTransform);
4336 :
4337 27 : std::vector<double> adfX(std::max(nPixels, nLines));
4338 27 : std::vector<double> adfY(adfX.size());
4339 27 : std::vector<double> adfZ(adfX.size());
4340 27 : std::vector<int> abSuccess(adfX.size());
4341 :
4342 : const auto DetectBlankBorder =
4343 108 : [&](int nValues,
4344 : std::function<bool(double, double)> funcIsOK)
4345 : {
4346 108 : if (nValues > 3)
4347 : {
4348 : // First try with just a subsample of 3 points
4349 52 : double adf3X[3] = {adfX[0], adfX[nValues / 2],
4350 52 : adfX[nValues - 1]};
4351 52 : double adf3Y[3] = {adfY[0], adfY[nValues / 2],
4352 52 : adfY[nValues - 1]};
4353 52 : double adf3Z[3] = {0};
4354 52 : GDALGenImgProjTransform(hUniqueTransformArg.get(), TRUE,
4355 : 3, &adf3X[0], &adf3Y[0],
4356 52 : &adf3Z[0], &abSuccess[0]);
4357 66 : for (int i = 0; i < 3; ++i)
4358 : {
4359 63 : if (abSuccess[i] && funcIsOK(adf3X[i], adf3Y[i]))
4360 : {
4361 49 : return false;
4362 : }
4363 : }
4364 : }
4365 :
4366 : // Do on full border to confirm
4367 59 : GDALGenImgProjTransform(hUniqueTransformArg.get(), TRUE,
4368 59 : nValues, &adfX[0], &adfY[0],
4369 59 : &adfZ[0], &abSuccess[0]);
4370 125 : for (int i = 0; i < nValues; ++i)
4371 : {
4372 122 : if (abSuccess[i] && funcIsOK(adfX[i], adfY[i]))
4373 : {
4374 56 : return false;
4375 : }
4376 : }
4377 :
4378 3 : return true;
4379 27 : };
4380 :
4381 9912 : for (int i = 0; i < nPixels; ++i)
4382 : {
4383 9885 : adfX[i] = i + 0.5;
4384 9885 : adfY[i] = 0.5;
4385 9885 : adfZ[i] = 0;
4386 : }
4387 27 : const bool bTopBlankLine = DetectBlankBorder(
4388 42 : nPixels, [](double, double y) { return y >= 0; });
4389 :
4390 9912 : for (int i = 0; i < nPixels; ++i)
4391 : {
4392 9885 : adfX[i] = i + 0.5;
4393 9885 : adfY[i] = nLines - 0.5;
4394 9885 : adfZ[i] = 0;
4395 : }
4396 27 : const int nSrcLines = GDALGetRasterYSize(pahSrcDS[0]);
4397 : const bool bBottomBlankLine =
4398 27 : DetectBlankBorder(nPixels, [nSrcLines](double, double y)
4399 27 : { return y <= nSrcLines; });
4400 :
4401 5139 : for (int i = 0; i < nLines; ++i)
4402 : {
4403 5112 : adfX[i] = 0.5;
4404 5112 : adfY[i] = i + 0.5;
4405 5112 : adfZ[i] = 0;
4406 : }
4407 27 : const bool bLeftBlankCol = DetectBlankBorder(
4408 58 : nLines, [](double x, double) { return x >= 0; });
4409 :
4410 5139 : for (int i = 0; i < nLines; ++i)
4411 : {
4412 5112 : adfX[i] = nPixels - 0.5;
4413 5112 : adfY[i] = i + 0.5;
4414 5112 : adfZ[i] = 0;
4415 : }
4416 27 : const int nSrcCols = GDALGetRasterXSize(pahSrcDS[0]);
4417 : const bool bRightBlankCol =
4418 27 : DetectBlankBorder(nLines, [nSrcCols](double x, double)
4419 54 : { return x <= nSrcCols; });
4420 :
4421 27 : if (!bTopBlankLine && !bBottomBlankLine && !bLeftBlankCol &&
4422 25 : !bRightBlankCol)
4423 25 : break;
4424 :
4425 2 : if (bTopBlankLine)
4426 : {
4427 1 : if (psOptions->dfMaxY - psOptions->dfMinY <=
4428 1 : 2 * psOptions->dfYRes)
4429 0 : break;
4430 1 : psOptions->dfMaxY -= psOptions->dfYRes;
4431 : }
4432 2 : if (bBottomBlankLine)
4433 : {
4434 0 : if (psOptions->dfMaxY - psOptions->dfMinY <=
4435 0 : 2 * psOptions->dfYRes)
4436 0 : break;
4437 0 : psOptions->dfMinY += psOptions->dfYRes;
4438 : }
4439 2 : if (bLeftBlankCol)
4440 : {
4441 1 : if (psOptions->dfMaxX - psOptions->dfMinX <=
4442 1 : 2 * psOptions->dfXRes)
4443 0 : break;
4444 1 : psOptions->dfMinX += psOptions->dfXRes;
4445 : }
4446 2 : if (bRightBlankCol)
4447 : {
4448 1 : if (psOptions->dfMaxX - psOptions->dfMinX <=
4449 1 : 2 * psOptions->dfXRes)
4450 0 : break;
4451 1 : psOptions->dfMaxX -= psOptions->dfXRes;
4452 : }
4453 : }
4454 : }
4455 :
4456 93 : if (!UpdateGeoTransformandAndPixelLines())
4457 93 : return nullptr;
4458 : }
4459 :
4460 936 : else if (psOptions->nForcePixels != 0 && psOptions->nForceLines != 0)
4461 : {
4462 142 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4463 104 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4464 : {
4465 104 : psOptions->dfMinX = dfWrkMinX;
4466 104 : psOptions->dfMaxX = dfWrkMaxX;
4467 104 : psOptions->dfMaxY = dfWrkMaxY;
4468 104 : psOptions->dfMinY = dfWrkMinY;
4469 : }
4470 :
4471 142 : psOptions->dfXRes =
4472 142 : (psOptions->dfMaxX - psOptions->dfMinX) / psOptions->nForcePixels;
4473 142 : psOptions->dfYRes =
4474 142 : (psOptions->dfMaxY - psOptions->dfMinY) / psOptions->nForceLines;
4475 :
4476 142 : adfDstGeoTransform[0] = psOptions->dfMinX;
4477 142 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4478 142 : adfDstGeoTransform[1] = psOptions->dfXRes;
4479 142 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4480 :
4481 142 : nPixels = psOptions->nForcePixels;
4482 142 : nLines = psOptions->nForceLines;
4483 : }
4484 :
4485 794 : else if (psOptions->nForcePixels != 0)
4486 : {
4487 5 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4488 4 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4489 : {
4490 4 : psOptions->dfMinX = dfWrkMinX;
4491 4 : psOptions->dfMaxX = dfWrkMaxX;
4492 4 : psOptions->dfMaxY = dfWrkMaxY;
4493 4 : psOptions->dfMinY = dfWrkMinY;
4494 : }
4495 :
4496 5 : psOptions->dfXRes =
4497 5 : (psOptions->dfMaxX - psOptions->dfMinX) / psOptions->nForcePixels;
4498 5 : psOptions->dfYRes = psOptions->dfXRes;
4499 :
4500 5 : adfDstGeoTransform[0] = psOptions->dfMinX;
4501 5 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4502 5 : adfDstGeoTransform[1] = psOptions->dfXRes;
4503 10 : adfDstGeoTransform[5] = (psOptions->dfMaxY > psOptions->dfMinY)
4504 5 : ? -psOptions->dfYRes
4505 0 : : psOptions->dfYRes;
4506 :
4507 5 : nPixels = psOptions->nForcePixels;
4508 5 : const double dfLines = ComputeLinesFromResAndExtent();
4509 5 : if (dfLines > INT_MAX)
4510 : {
4511 1 : CPLError(CE_Failure, CPLE_AppDefined,
4512 : "Too large output raster size: %d x %f", nPixels, dfLines);
4513 1 : return nullptr;
4514 : }
4515 4 : nLines = static_cast<int>(dfLines);
4516 : }
4517 :
4518 789 : else if (psOptions->nForceLines != 0)
4519 : {
4520 3 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4521 3 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4522 : {
4523 3 : psOptions->dfMinX = dfWrkMinX;
4524 3 : psOptions->dfMaxX = dfWrkMaxX;
4525 3 : psOptions->dfMaxY = dfWrkMaxY;
4526 3 : psOptions->dfMinY = dfWrkMinY;
4527 : }
4528 :
4529 3 : psOptions->dfYRes =
4530 3 : (psOptions->dfMaxY - psOptions->dfMinY) / psOptions->nForceLines;
4531 3 : psOptions->dfXRes = std::fabs(psOptions->dfYRes);
4532 :
4533 3 : adfDstGeoTransform[0] = psOptions->dfMinX;
4534 3 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4535 3 : adfDstGeoTransform[1] = psOptions->dfXRes;
4536 3 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4537 :
4538 3 : const double dfPixels = ComputePixelsFromResAndExtent();
4539 3 : nLines = psOptions->nForceLines;
4540 3 : if (dfPixels > INT_MAX)
4541 : {
4542 1 : CPLError(CE_Failure, CPLE_AppDefined,
4543 : "Too large output raster size: %f x %d", dfPixels, nLines);
4544 1 : return nullptr;
4545 : }
4546 2 : nPixels = static_cast<int>(dfPixels);
4547 : }
4548 :
4549 786 : else if (psOptions->dfMinX != 0.0 || psOptions->dfMinY != 0.0 ||
4550 710 : psOptions->dfMaxX != 0.0 || psOptions->dfMaxY != 0.0)
4551 : {
4552 77 : psOptions->dfXRes = adfDstGeoTransform[1];
4553 77 : psOptions->dfYRes = fabs(adfDstGeoTransform[5]);
4554 :
4555 77 : const double dfPixels = ComputePixelsFromResAndExtent();
4556 77 : const double dfLines = ComputeLinesFromResAndExtent();
4557 77 : if (dfPixels > INT_MAX || dfLines > INT_MAX)
4558 : {
4559 1 : CPLError(CE_Failure, CPLE_AppDefined,
4560 : "Too large output raster size: %f x %f", dfPixels,
4561 : dfLines);
4562 1 : return nullptr;
4563 : }
4564 :
4565 76 : nPixels = static_cast<int>(dfPixels);
4566 76 : nLines = static_cast<int>(dfLines);
4567 :
4568 76 : psOptions->dfXRes = (psOptions->dfMaxX - psOptions->dfMinX) / nPixels;
4569 76 : psOptions->dfYRes = (psOptions->dfMaxY - psOptions->dfMinY) / nLines;
4570 :
4571 76 : adfDstGeoTransform[0] = psOptions->dfMinX;
4572 76 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4573 76 : adfDstGeoTransform[1] = psOptions->dfXRes;
4574 76 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4575 : }
4576 :
4577 1050 : if (EQUAL(pszFormat, "GTiff"))
4578 : {
4579 :
4580 : /* --------------------------------------------------------------------
4581 : */
4582 : /* Automatically set PHOTOMETRIC=RGB for GTiff when appropriate */
4583 : /* --------------------------------------------------------------------
4584 : */
4585 455 : if (apeColorInterpretations.size() >= 3 &&
4586 35 : apeColorInterpretations[0] == GCI_RedBand &&
4587 35 : apeColorInterpretations[1] == GCI_GreenBand &&
4588 525 : apeColorInterpretations[2] == GCI_BlueBand &&
4589 35 : aosCreateOptions.FetchNameValue("PHOTOMETRIC") == nullptr)
4590 : {
4591 35 : aosCreateOptions.SetNameValue("PHOTOMETRIC", "RGB");
4592 :
4593 : // Preserve potential ALPHA=PREMULTIPLIED from source alpha band
4594 70 : if (aosCreateOptions.FetchNameValue("ALPHA") == nullptr &&
4595 35 : apeColorInterpretations.size() == 4 &&
4596 75 : apeColorInterpretations[3] == GCI_AlphaBand &&
4597 5 : GDALGetRasterCount(pahSrcDS[0]) == 4)
4598 : {
4599 : const char *pszAlpha =
4600 5 : GDALGetMetadataItem(GDALGetRasterBand(pahSrcDS[0], 4),
4601 : "ALPHA", "IMAGE_STRUCTURE");
4602 5 : if (pszAlpha)
4603 : {
4604 1 : aosCreateOptions.SetNameValue("ALPHA", pszAlpha);
4605 : }
4606 : }
4607 : }
4608 :
4609 : /* The GTiff driver now supports writing band color interpretation */
4610 : /* in the TIFF_GDAL_METADATA tag */
4611 455 : bSetColorInterpretation = true;
4612 : }
4613 :
4614 : /* -------------------------------------------------------------------- */
4615 : /* Create the output file. */
4616 : /* -------------------------------------------------------------------- */
4617 1050 : if (!psOptions->bQuiet)
4618 78 : printf("Creating output file that is %dP x %dL.\n", nPixels, nLines);
4619 :
4620 1050 : hDstDS = GDALCreate(hDriver, pszFilename, nPixels, nLines, nDstBandCount,
4621 1050 : eDT, aosCreateOptions.List());
4622 :
4623 1050 : if (hDstDS == nullptr)
4624 : {
4625 1 : return nullptr;
4626 : }
4627 :
4628 1049 : if (psOptions->bDeleteOutputFileOnceCreated)
4629 : {
4630 12 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
4631 6 : GDALDeleteDataset(hDriver, pszFilename);
4632 : }
4633 :
4634 : /* -------------------------------------------------------------------- */
4635 : /* Write out the projection definition. */
4636 : /* -------------------------------------------------------------------- */
4637 1049 : const char *pszDstMethod = CSLFetchNameValue(papszTO, "DST_METHOD");
4638 1049 : if (pszDstMethod == nullptr || !EQUAL(pszDstMethod, "NO_GEOTRANSFORM"))
4639 : {
4640 1047 : OGRSpatialReference oTargetSRS;
4641 1047 : oTargetSRS.SetFromUserInput(osThisTargetSRS);
4642 1047 : oTargetSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4643 :
4644 1047 : if (oTargetSRS.IsDynamic())
4645 : {
4646 620 : double dfCoordEpoch = CPLAtof(CSLFetchNameValueDef(
4647 : papszTO, "DST_COORDINATE_EPOCH",
4648 : CSLFetchNameValueDef(papszTO, "COORDINATE_EPOCH", "0")));
4649 620 : if (dfCoordEpoch == 0)
4650 : {
4651 : const OGRSpatialReferenceH hSrcSRS =
4652 620 : GDALGetSpatialRef(pahSrcDS[0]);
4653 620 : const char *pszMethod = FetchSrcMethod(papszTO);
4654 620 : if (hSrcSRS &&
4655 7 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")))
4656 : {
4657 577 : dfCoordEpoch = OSRGetCoordinateEpoch(hSrcSRS);
4658 : }
4659 : }
4660 620 : if (dfCoordEpoch > 0)
4661 1 : oTargetSRS.SetCoordinateEpoch(dfCoordEpoch);
4662 : }
4663 :
4664 1047 : if (GDALSetSpatialRef(hDstDS, OGRSpatialReference::ToHandle(
4665 2094 : &oTargetSRS)) == CE_Failure ||
4666 1047 : GDALSetGeoTransform(hDstDS, adfDstGeoTransform) == CE_Failure)
4667 : {
4668 0 : GDALClose(hDstDS);
4669 0 : return nullptr;
4670 1047 : }
4671 : }
4672 : else
4673 : {
4674 2 : adfDstGeoTransform[3] += adfDstGeoTransform[5] * nLines;
4675 2 : adfDstGeoTransform[5] = fabs(adfDstGeoTransform[5]);
4676 : }
4677 :
4678 1049 : if (hUniqueTransformArg && bUpdateTransformerWithDestGT)
4679 : {
4680 1023 : GDALSetGenImgProjTransformerDstGeoTransform(hUniqueTransformArg.get(),
4681 : adfDstGeoTransform);
4682 :
4683 1023 : void *pTransformerArg = hUniqueTransformArg.get();
4684 1023 : if (GDALIsTransformer(pTransformerArg,
4685 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4686 : {
4687 : // Detect if there is a change of coordinate operation in the area of
4688 : // interest. The underlying proj_trans_get_last_used_operation() is
4689 : // quite costly due to using proj_clone() internally, so only do that
4690 : // on a restricted set of points.
4691 1023 : GDALGenImgProjTransformInfo *psTransformInfo{
4692 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg)};
4693 1023 : GDALTransformerInfo *psInfo = &psTransformInfo->sTI;
4694 :
4695 1023 : void *pReprojectArg = psTransformInfo->pReprojectArg;
4696 1023 : if (GDALIsTransformer(pReprojectArg,
4697 : GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4698 : {
4699 1 : const auto *pApproxInfo =
4700 : static_cast<const GDALApproxTransformInfo *>(pReprojectArg);
4701 1 : pReprojectArg = pApproxInfo->pBaseCBData;
4702 : }
4703 :
4704 1023 : if (GDALIsTransformer(pReprojectArg,
4705 : GDAL_REPROJECTION_TRANSFORMER_CLASS_NAME))
4706 : {
4707 618 : const GDALReprojectionTransformInfo *psRTI =
4708 : static_cast<const GDALReprojectionTransformInfo *>(
4709 : pReprojectArg);
4710 618 : if (psRTI->poReverseTransform)
4711 : {
4712 1236 : std::vector<double> adfX, adfY, adfZ;
4713 1236 : std::vector<int> abSuccess;
4714 :
4715 618 : GDALDatasetH hSrcDS = pahSrcDS[0];
4716 :
4717 : // Sample points on a N x N grid in the source raster
4718 618 : constexpr int N = 10;
4719 618 : const int nSrcXSize = GDALGetRasterXSize(hSrcDS);
4720 618 : const int nSrcYSize = GDALGetRasterYSize(hSrcDS);
4721 7416 : for (int j = 0; j <= N; ++j)
4722 : {
4723 81576 : for (int i = 0; i <= N; ++i)
4724 : {
4725 74778 : adfX.push_back(static_cast<double>(i) / N *
4726 : nSrcXSize);
4727 74778 : adfY.push_back(static_cast<double>(j) / N *
4728 : nSrcYSize);
4729 74778 : adfZ.push_back(0);
4730 74778 : abSuccess.push_back(0);
4731 : }
4732 : }
4733 :
4734 : {
4735 1236 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
4736 :
4737 : // Transform from source raster coordinates to target raster
4738 : // coordinates
4739 618 : psInfo->pfnTransform(hUniqueTransformArg.get(), FALSE,
4740 618 : static_cast<int>(adfX.size()),
4741 618 : &adfX[0], &adfY[0], &adfZ[0],
4742 618 : &abSuccess[0]);
4743 :
4744 618 : const int nDstXSize = nPixels;
4745 618 : const int nDstYSize = nLines;
4746 :
4747 : // Clamp target raster coordinates
4748 75396 : for (size_t i = 0; i < adfX.size(); ++i)
4749 : {
4750 74778 : if (adfX[i] < 0)
4751 1100 : adfX[i] = 0;
4752 74778 : if (adfX[i] > nDstXSize)
4753 2981 : adfX[i] = nDstXSize;
4754 74778 : if (adfY[i] < 0)
4755 1414 : adfY[i] = 0;
4756 74778 : if (adfY[i] > nDstYSize)
4757 3568 : adfY[i] = nDstYSize;
4758 : }
4759 :
4760 : // Start recording if different coordinate operations are
4761 : // going to be used
4762 618 : OGRProjCTDifferentOperationsStart(
4763 618 : psRTI->poReverseTransform);
4764 :
4765 : // Transform back to source raster coordinates.
4766 618 : psInfo->pfnTransform(hUniqueTransformArg.get(), TRUE,
4767 618 : static_cast<int>(adfX.size()),
4768 618 : &adfX[0], &adfY[0], &adfZ[0],
4769 618 : &abSuccess[0]);
4770 : }
4771 :
4772 618 : if (OGRProjCTDifferentOperationsUsed(
4773 618 : psRTI->poReverseTransform))
4774 : {
4775 0 : CPLError(
4776 : CE_Warning, CPLE_AppDefined,
4777 : "Several coordinate operations are going to be "
4778 : "used. Artifacts may appear. You may consider "
4779 : "using the -to ALLOW_BALLPARK=NO and/or "
4780 : "-to ONLY_BEST=YES transform options, or specify "
4781 : "a particular coordinate operation with -ct");
4782 : }
4783 :
4784 : // Stop recording
4785 618 : OGRProjCTDifferentOperationsStop(psRTI->poReverseTransform);
4786 : }
4787 : }
4788 : }
4789 : }
4790 :
4791 : /* -------------------------------------------------------------------- */
4792 : /* Try to set color interpretation of source bands to target */
4793 : /* dataset. */
4794 : /* FIXME? We should likely do that for other drivers than VRT & */
4795 : /* GTiff but it might create spurious .aux.xml files (at least */
4796 : /* with HFA, and netCDF) */
4797 : /* -------------------------------------------------------------------- */
4798 1049 : if (bVRT || bSetColorInterpretation)
4799 : {
4800 632 : int nBandsToCopy = static_cast<int>(apeColorInterpretations.size());
4801 632 : if (psOptions->bEnableSrcAlpha)
4802 30 : nBandsToCopy--;
4803 1447 : for (int iBand = 0; iBand < nBandsToCopy; iBand++)
4804 : {
4805 815 : GDALSetRasterColorInterpretation(
4806 : GDALGetRasterBand(hDstDS, iBand + 1),
4807 815 : apeColorInterpretations[iBand]);
4808 : }
4809 : }
4810 :
4811 : /* -------------------------------------------------------------------- */
4812 : /* Try to set color interpretation of output file alpha band. */
4813 : /* -------------------------------------------------------------------- */
4814 1049 : if (psOptions->bEnableDstAlpha)
4815 : {
4816 99 : GDALSetRasterColorInterpretation(
4817 : GDALGetRasterBand(hDstDS, nDstBandCount), GCI_AlphaBand);
4818 : }
4819 :
4820 : /* -------------------------------------------------------------------- */
4821 : /* Copy the raster attribute table, if required. */
4822 : /* -------------------------------------------------------------------- */
4823 1049 : if (hRAT != nullptr)
4824 : {
4825 0 : GDALSetDefaultRAT(GDALGetRasterBand(hDstDS, 1), hRAT);
4826 : }
4827 :
4828 : /* -------------------------------------------------------------------- */
4829 : /* Copy the color table, if required. */
4830 : /* -------------------------------------------------------------------- */
4831 1049 : if (poCT)
4832 : {
4833 6 : GDALSetRasterColorTable(GDALGetRasterBand(hDstDS, 1),
4834 : GDALColorTable::ToHandle(poCT.get()));
4835 : }
4836 :
4837 : /* -------------------------------------------------------------------- */
4838 : /* Copy scale/offset if found on source */
4839 : /* -------------------------------------------------------------------- */
4840 1049 : if (nSrcCount == 1)
4841 : {
4842 1025 : GDALDataset *poSrcDS = GDALDataset::FromHandle(pahSrcDS[0]);
4843 1025 : GDALDataset *poDstDS = GDALDataset::FromHandle(hDstDS);
4844 :
4845 1025 : int nBandsToCopy = nDstBandCount;
4846 1025 : if (psOptions->bEnableDstAlpha)
4847 96 : nBandsToCopy--;
4848 1025 : nBandsToCopy = std::min(nBandsToCopy, poSrcDS->GetRasterCount());
4849 :
4850 2419 : for (int i = 0; i < nBandsToCopy; i++)
4851 : {
4852 1434 : auto poSrcBand = poSrcDS->GetRasterBand(
4853 1394 : psOptions->anSrcBands.empty() ? i + 1
4854 40 : : psOptions->anSrcBands[i]);
4855 1434 : auto poDstBand = poDstDS->GetRasterBand(
4856 1394 : psOptions->anDstBands.empty() ? i + 1
4857 40 : : psOptions->anDstBands[i]);
4858 1394 : if (poSrcBand && poDstBand)
4859 : {
4860 1392 : int bHasScale = FALSE;
4861 1392 : const double dfScale = poSrcBand->GetScale(&bHasScale);
4862 1392 : if (bHasScale)
4863 42 : poDstBand->SetScale(dfScale);
4864 :
4865 1392 : int bHasOffset = FALSE;
4866 1392 : const double dfOffset = poSrcBand->GetOffset(&bHasOffset);
4867 1392 : if (bHasOffset)
4868 42 : poDstBand->SetOffset(dfOffset);
4869 : }
4870 : }
4871 : }
4872 :
4873 1049 : return hDstDS;
4874 : }
4875 :
4876 : /************************************************************************/
4877 : /* GeoTransform_Transformer() */
4878 : /* */
4879 : /* Convert points from georef coordinates to pixel/line based */
4880 : /* on a geotransform. */
4881 : /************************************************************************/
4882 : namespace
4883 : {
4884 : class CutlineTransformer final : public OGRCoordinateTransformation
4885 : {
4886 : CPL_DISALLOW_COPY_ASSIGN(CutlineTransformer)
4887 :
4888 : public:
4889 : void *hSrcImageTransformer = nullptr;
4890 :
4891 46 : explicit CutlineTransformer(void *hTransformArg)
4892 46 : : hSrcImageTransformer(hTransformArg)
4893 : {
4894 46 : }
4895 :
4896 0 : const OGRSpatialReference *GetSourceCS() const override
4897 : {
4898 0 : return nullptr;
4899 : }
4900 :
4901 185 : const OGRSpatialReference *GetTargetCS() const override
4902 : {
4903 185 : return nullptr;
4904 : }
4905 :
4906 : ~CutlineTransformer() override;
4907 :
4908 72 : virtual int Transform(size_t nCount, double *x, double *y, double *z,
4909 : double * /* t */, int *pabSuccess) override
4910 : {
4911 72 : CPLAssert(nCount <=
4912 : static_cast<size_t>(std::numeric_limits<int>::max()));
4913 72 : return GDALGenImgProjTransform(hSrcImageTransformer, TRUE,
4914 : static_cast<int>(nCount), x, y, z,
4915 72 : pabSuccess);
4916 : }
4917 :
4918 0 : OGRCoordinateTransformation *Clone() const override
4919 : {
4920 : return new CutlineTransformer(
4921 0 : GDALCloneTransformer(hSrcImageTransformer));
4922 : }
4923 :
4924 0 : OGRCoordinateTransformation *GetInverse() const override
4925 : {
4926 0 : return nullptr;
4927 : }
4928 : };
4929 :
4930 46 : CutlineTransformer::~CutlineTransformer()
4931 : {
4932 46 : GDALDestroyTransformer(hSrcImageTransformer);
4933 46 : }
4934 : } // namespace
4935 :
4936 312 : static double GetMaximumSegmentLength(const OGRGeometry *poGeom)
4937 : {
4938 312 : switch (wkbFlatten(poGeom->getGeometryType()))
4939 : {
4940 120 : case wkbLineString:
4941 : {
4942 120 : const OGRLineString *poLS = poGeom->toLineString();
4943 120 : double dfMaxSquaredLength = 0.0;
4944 16620 : for (int i = 0; i < poLS->getNumPoints() - 1; i++)
4945 : {
4946 16500 : double dfDeltaX = poLS->getX(i + 1) - poLS->getX(i);
4947 16500 : double dfDeltaY = poLS->getY(i + 1) - poLS->getY(i);
4948 16500 : double dfSquaredLength =
4949 16500 : dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY;
4950 16500 : dfMaxSquaredLength =
4951 16500 : std::max(dfMaxSquaredLength, dfSquaredLength);
4952 : }
4953 120 : return sqrt(dfMaxSquaredLength);
4954 : }
4955 :
4956 118 : case wkbPolygon:
4957 : {
4958 118 : const OGRPolygon *poPoly = poGeom->toPolygon();
4959 118 : double dfMaxLength = 0;
4960 238 : for (const auto *poRing : *poPoly)
4961 : {
4962 120 : dfMaxLength =
4963 120 : std::max(dfMaxLength, GetMaximumSegmentLength(poRing));
4964 : }
4965 118 : return dfMaxLength;
4966 : }
4967 :
4968 74 : case wkbMultiPolygon:
4969 : {
4970 74 : const OGRMultiPolygon *poMP = poGeom->toMultiPolygon();
4971 74 : double dfMaxLength = 0.0;
4972 152 : for (const auto *poPoly : *poMP)
4973 : {
4974 78 : dfMaxLength =
4975 78 : std::max(dfMaxLength, GetMaximumSegmentLength(poPoly));
4976 : }
4977 74 : return dfMaxLength;
4978 : }
4979 :
4980 0 : default:
4981 0 : CPLAssert(false);
4982 : return 0.0;
4983 : }
4984 : }
4985 :
4986 : /************************************************************************/
4987 : /* RemoveZeroWidthSlivers() */
4988 : /* */
4989 : /* Such slivers can cause issues after reprojection. */
4990 : /************************************************************************/
4991 :
4992 130 : static void RemoveZeroWidthSlivers(OGRGeometry *poGeom)
4993 : {
4994 130 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
4995 130 : if (eType == wkbMultiPolygon)
4996 : {
4997 32 : auto poMP = poGeom->toMultiPolygon();
4998 32 : int nNumGeometries = poMP->getNumGeometries();
4999 66 : for (int i = 0; i < nNumGeometries; /* incremented in loop */)
5000 : {
5001 34 : auto poPoly = poMP->getGeometryRef(i);
5002 34 : RemoveZeroWidthSlivers(poPoly);
5003 34 : if (poPoly->IsEmpty())
5004 : {
5005 1 : CPLDebug("WARP",
5006 : "RemoveZeroWidthSlivers: removing empty polygon");
5007 1 : poMP->removeGeometry(i, /* bDelete = */ true);
5008 1 : --nNumGeometries;
5009 : }
5010 : else
5011 : {
5012 33 : ++i;
5013 : }
5014 : }
5015 : }
5016 98 : else if (eType == wkbPolygon)
5017 : {
5018 48 : auto poPoly = poGeom->toPolygon();
5019 48 : if (auto poExteriorRing = poPoly->getExteriorRing())
5020 : {
5021 48 : RemoveZeroWidthSlivers(poExteriorRing);
5022 48 : if (poExteriorRing->getNumPoints() < 4)
5023 : {
5024 1 : poPoly->empty();
5025 1 : return;
5026 : }
5027 : }
5028 47 : int nNumInteriorRings = poPoly->getNumInteriorRings();
5029 49 : for (int i = 0; i < nNumInteriorRings; /* incremented in loop */)
5030 : {
5031 2 : auto poRing = poPoly->getInteriorRing(i);
5032 2 : RemoveZeroWidthSlivers(poRing);
5033 2 : if (poRing->getNumPoints() < 4)
5034 : {
5035 1 : CPLDebug(
5036 : "WARP",
5037 : "RemoveZeroWidthSlivers: removing empty interior ring");
5038 1 : constexpr int OFFSET_EXTERIOR_RING = 1;
5039 1 : poPoly->removeRing(i + OFFSET_EXTERIOR_RING,
5040 : /* bDelete = */ true);
5041 1 : --nNumInteriorRings;
5042 : }
5043 : else
5044 : {
5045 1 : ++i;
5046 : }
5047 : }
5048 : }
5049 50 : else if (eType == wkbLineString)
5050 : {
5051 50 : OGRLineString *poLS = poGeom->toLineString();
5052 50 : int numPoints = poLS->getNumPoints();
5053 691 : for (int i = 1; i < numPoints - 1;)
5054 : {
5055 641 : const double x1 = poLS->getX(i - 1);
5056 641 : const double y1 = poLS->getY(i - 1);
5057 641 : const double x2 = poLS->getX(i);
5058 641 : const double y2 = poLS->getY(i);
5059 641 : const double x3 = poLS->getX(i + 1);
5060 641 : const double y3 = poLS->getY(i + 1);
5061 641 : const double dx1 = x2 - x1;
5062 641 : const double dy1 = y2 - y1;
5063 641 : const double dx2 = x3 - x2;
5064 641 : const double dy2 = y3 - y2;
5065 641 : const double scalar_product = dx1 * dx2 + dy1 * dy2;
5066 641 : const double square_scalar_product =
5067 : scalar_product * scalar_product;
5068 641 : const double square_norm1 = dx1 * dx1 + dy1 * dy1;
5069 641 : const double square_norm2 = dx2 * dx2 + dy2 * dy2;
5070 641 : const double square_norm1_mult_norm2 = square_norm1 * square_norm2;
5071 641 : if (scalar_product < 0 &&
5072 75 : fabs(square_scalar_product - square_norm1_mult_norm2) <=
5073 75 : 1e-15 * square_norm1_mult_norm2)
5074 : {
5075 5 : CPLDebug("WARP",
5076 : "RemoveZeroWidthSlivers: removing point %.10g %.10g",
5077 : x2, y2);
5078 5 : poLS->removePoint(i);
5079 5 : numPoints--;
5080 : }
5081 : else
5082 : {
5083 636 : ++i;
5084 : }
5085 : }
5086 : }
5087 : }
5088 :
5089 : /************************************************************************/
5090 : /* TransformCutlineToSource() */
5091 : /* */
5092 : /* Transform cutline from its SRS to source pixel/line coordinates.*/
5093 : /************************************************************************/
5094 46 : static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,
5095 : OGRGeometry *poCutline,
5096 : char ***ppapszWarpOptions,
5097 : CSLConstList papszTO_In)
5098 :
5099 : {
5100 46 : RemoveZeroWidthSlivers(poCutline);
5101 :
5102 92 : auto poMultiPolygon = std::unique_ptr<OGRGeometry>(poCutline->clone());
5103 :
5104 : /* -------------------------------------------------------------------- */
5105 : /* Checkout that if there's a cutline SRS, there's also a raster */
5106 : /* one. */
5107 : /* -------------------------------------------------------------------- */
5108 46 : std::unique_ptr<OGRSpatialReference> poRasterSRS;
5109 : const CPLString osProjection =
5110 92 : GetSrcDSProjection(GDALDataset::ToHandle(poSrcDS), papszTO_In);
5111 46 : if (!osProjection.empty())
5112 : {
5113 42 : poRasterSRS = std::make_unique<OGRSpatialReference>();
5114 42 : poRasterSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
5115 42 : if (poRasterSRS->SetFromUserInput(osProjection) != OGRERR_NONE)
5116 : {
5117 0 : poRasterSRS.reset();
5118 : }
5119 : }
5120 :
5121 46 : std::unique_ptr<OGRSpatialReference> poDstSRS;
5122 46 : const char *pszThisTargetSRS = CSLFetchNameValue(papszTO_In, "DST_SRS");
5123 46 : if (pszThisTargetSRS)
5124 : {
5125 8 : poDstSRS = std::make_unique<OGRSpatialReference>();
5126 8 : poDstSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
5127 8 : if (poDstSRS->SetFromUserInput(pszThisTargetSRS) != OGRERR_NONE)
5128 : {
5129 0 : return CE_Failure;
5130 : }
5131 : }
5132 38 : else if (poRasterSRS)
5133 : {
5134 34 : poDstSRS.reset(poRasterSRS->Clone());
5135 : }
5136 :
5137 : /* -------------------------------------------------------------------- */
5138 : /* Extract the cutline SRS. */
5139 : /* -------------------------------------------------------------------- */
5140 : const OGRSpatialReference *poCutlineSRS =
5141 46 : poMultiPolygon->getSpatialReference();
5142 :
5143 : /* -------------------------------------------------------------------- */
5144 : /* Detect if there's no transform at all involved, in which case */
5145 : /* we can avoid densification. */
5146 : /* -------------------------------------------------------------------- */
5147 46 : bool bMayNeedDensify = true;
5148 88 : if (poRasterSRS && poCutlineSRS && poRasterSRS->IsSame(poCutlineSRS) &&
5149 32 : poSrcDS->GetGCPCount() == 0 && !poSrcDS->GetMetadata("RPC") &&
5150 28 : !poSrcDS->GetMetadata("GEOLOCATION") &&
5151 116 : !CSLFetchNameValue(papszTO_In, "GEOLOC_ARRAY") &&
5152 28 : !CSLFetchNameValue(papszTO_In, "SRC_GEOLOC_ARRAY"))
5153 : {
5154 56 : CPLStringList aosTOTmp(papszTO_In);
5155 28 : aosTOTmp.SetNameValue("SRC_SRS", nullptr);
5156 28 : aosTOTmp.SetNameValue("DST_SRS", nullptr);
5157 28 : if (aosTOTmp.size() == 0)
5158 : {
5159 22 : bMayNeedDensify = false;
5160 : }
5161 : }
5162 :
5163 : /* -------------------------------------------------------------------- */
5164 : /* Compare source raster SRS and cutline SRS */
5165 : /* -------------------------------------------------------------------- */
5166 46 : if (poRasterSRS && poCutlineSRS)
5167 : {
5168 : /* OK, we will reproject */
5169 : }
5170 8 : else if (poRasterSRS && !poCutlineSRS)
5171 : {
5172 4 : CPLError(
5173 : CE_Warning, CPLE_AppDefined,
5174 : "the source raster dataset has a SRS, but the cutline features\n"
5175 : "not. We assume that the cutline coordinates are expressed in the "
5176 : "destination SRS.\n"
5177 : "If not, cutline results may be incorrect.");
5178 : }
5179 4 : else if (!poRasterSRS && poCutlineSRS)
5180 : {
5181 0 : CPLError(CE_Warning, CPLE_AppDefined,
5182 : "the input vector layer has a SRS, but the source raster "
5183 : "dataset does not.\n"
5184 : "Cutline results may be incorrect.");
5185 : }
5186 :
5187 : auto poCTCutlineToSrc = CreateCTCutlineToSrc(
5188 92 : poRasterSRS.get(), poDstSRS.get(), poCutlineSRS, papszTO_In);
5189 :
5190 92 : CPLStringList aosTO(papszTO_In);
5191 :
5192 46 : if (pszThisTargetSRS && !osProjection.empty())
5193 : {
5194 : // Avoid any reprojection when using the GenImgProjTransformer
5195 8 : aosTO.SetNameValue("DST_SRS", osProjection.c_str());
5196 : }
5197 46 : aosTO.SetNameValue("SRC_COORDINATE_EPOCH", nullptr);
5198 46 : aosTO.SetNameValue("DST_COORDINATE_EPOCH", nullptr);
5199 46 : aosTO.SetNameValue("COORDINATE_OPERATION", nullptr);
5200 :
5201 : /* -------------------------------------------------------------------- */
5202 : /* It may be unwise to let the mask geometry be re-wrapped by */
5203 : /* the CENTER_LONG machinery as this can easily screw up world */
5204 : /* spanning masks and invert the mask topology. */
5205 : /* -------------------------------------------------------------------- */
5206 46 : aosTO.SetNameValue("INSERT_CENTER_LONG", "FALSE");
5207 :
5208 : /* -------------------------------------------------------------------- */
5209 : /* Transform the geometry to pixel/line coordinates. */
5210 : /* -------------------------------------------------------------------- */
5211 : /* The cutline transformer will *invert* the hSrcImageTransformer */
5212 : /* so it will convert from the source SRS to the source pixel/line */
5213 : /* coordinates */
5214 : CutlineTransformer oTransformer(GDALCreateGenImgProjTransformer2(
5215 92 : GDALDataset::ToHandle(poSrcDS), nullptr, aosTO.List()));
5216 :
5217 46 : if (oTransformer.hSrcImageTransformer == nullptr)
5218 : {
5219 0 : return CE_Failure;
5220 : }
5221 :
5222 : // Some transforms like RPC can transform a valid geometry into an invalid
5223 : // one if the node density of the input geometry isn't sufficient before
5224 : // reprojection. So after an initial reprojection, we check that the
5225 : // maximum length of a segment is no longer than 1 pixel, and if not,
5226 : // we densify the input geometry before doing a new reprojection
5227 : const double dfMaxLengthInSpatUnits =
5228 46 : GetMaximumSegmentLength(poMultiPolygon.get());
5229 46 : OGRErr eErr = OGRERR_NONE;
5230 46 : if (poCTCutlineToSrc)
5231 : {
5232 12 : poMultiPolygon.reset(OGRGeometryFactory::transformWithOptions(
5233 6 : poMultiPolygon.get(), poCTCutlineToSrc.get(), nullptr));
5234 6 : if (!poMultiPolygon)
5235 : {
5236 0 : eErr = OGRERR_FAILURE;
5237 0 : poMultiPolygon.reset(poCutline->clone());
5238 0 : poMultiPolygon->transform(poCTCutlineToSrc.get());
5239 : }
5240 : }
5241 46 : if (poMultiPolygon->transform(&oTransformer) != OGRERR_NONE)
5242 : {
5243 0 : CPLError(CE_Failure, CPLE_AppDefined,
5244 : "poMultiPolygon->transform(&oTransformer) failed at line %d",
5245 : __LINE__);
5246 0 : eErr = OGRERR_FAILURE;
5247 : }
5248 : const double dfInitialMaxLengthInPixels =
5249 46 : GetMaximumSegmentLength(poMultiPolygon.get());
5250 :
5251 46 : CPLPushErrorHandler(CPLQuietErrorHandler);
5252 : const bool bWasValidInitially =
5253 46 : ValidateCutline(poMultiPolygon.get(), false);
5254 46 : CPLPopErrorHandler();
5255 46 : if (!bWasValidInitially)
5256 : {
5257 3 : CPLDebug("WARP", "Cutline is not valid after initial reprojection");
5258 3 : char *pszWKT = nullptr;
5259 3 : poMultiPolygon->exportToWkt(&pszWKT);
5260 3 : CPLDebug("GDALWARP", "WKT = \"%s\"", pszWKT ? pszWKT : "(null)");
5261 3 : CPLFree(pszWKT);
5262 : }
5263 :
5264 46 : bool bDensify = false;
5265 46 : if (bMayNeedDensify && eErr == OGRERR_NONE &&
5266 : dfInitialMaxLengthInPixels > 1.0)
5267 : {
5268 : const char *pszDensifyCutline =
5269 22 : CPLGetConfigOption("GDALWARP_DENSIFY_CUTLINE", "YES");
5270 22 : if (EQUAL(pszDensifyCutline, "ONLY_IF_INVALID"))
5271 : {
5272 1 : bDensify = (OGRGeometryFactory::haveGEOS() && !bWasValidInitially);
5273 : }
5274 21 : else if (CSLFetchNameValue(*ppapszWarpOptions, "CUTLINE_BLEND_DIST") !=
5275 21 : nullptr &&
5276 0 : CPLGetConfigOption("GDALWARP_DENSIFY_CUTLINE", nullptr) ==
5277 : nullptr)
5278 : {
5279 : // TODO: we should only emit this message if a
5280 : // transform/reprojection will be actually done
5281 0 : CPLDebug("WARP",
5282 : "Densification of cutline could perhaps be useful but as "
5283 : "CUTLINE_BLEND_DIST is used, this could be very slow. So "
5284 : "disabled "
5285 : "unless GDALWARP_DENSIFY_CUTLINE=YES is explicitly "
5286 : "specified as configuration option");
5287 : }
5288 : else
5289 : {
5290 21 : bDensify = CPLTestBool(pszDensifyCutline);
5291 : }
5292 : }
5293 46 : if (bDensify)
5294 : {
5295 21 : CPLDebug("WARP",
5296 : "Cutline maximum segment size was %.0f pixel after "
5297 : "reprojection to source coordinates.",
5298 : dfInitialMaxLengthInPixels);
5299 :
5300 : // Densify and reproject with the aim of having a 1 pixel density
5301 21 : double dfSegmentSize =
5302 : dfMaxLengthInSpatUnits / dfInitialMaxLengthInPixels;
5303 21 : const int MAX_ITERATIONS = 10;
5304 22 : for (int i = 0; i < MAX_ITERATIONS; i++)
5305 : {
5306 22 : poMultiPolygon.reset(poCutline->clone());
5307 22 : poMultiPolygon->segmentize(dfSegmentSize);
5308 22 : if (i == MAX_ITERATIONS - 1)
5309 : {
5310 0 : char *pszWKT = nullptr;
5311 0 : poMultiPolygon->exportToWkt(&pszWKT);
5312 0 : CPLDebug("WARP",
5313 : "WKT of polygon after densification with segment size "
5314 : "= %f: %s",
5315 : dfSegmentSize, pszWKT);
5316 0 : CPLFree(pszWKT);
5317 : }
5318 22 : eErr = OGRERR_NONE;
5319 22 : if (poCTCutlineToSrc)
5320 : {
5321 12 : poMultiPolygon.reset(OGRGeometryFactory::transformWithOptions(
5322 6 : poMultiPolygon.get(), poCTCutlineToSrc.get(), nullptr));
5323 6 : if (!poMultiPolygon)
5324 : {
5325 0 : eErr = OGRERR_FAILURE;
5326 0 : break;
5327 : }
5328 : }
5329 22 : if (poMultiPolygon->transform(&oTransformer) != OGRERR_NONE)
5330 0 : eErr = OGRERR_FAILURE;
5331 22 : if (eErr == OGRERR_NONE)
5332 : {
5333 : const double dfMaxLengthInPixels =
5334 22 : GetMaximumSegmentLength(poMultiPolygon.get());
5335 22 : if (bWasValidInitially)
5336 : {
5337 : // In some cases, the densification itself results in a
5338 : // reprojected invalid polygon due to the non-linearity of
5339 : // RPC DEM transformation, so in those cases, try a less
5340 : // dense cutline
5341 20 : CPLPushErrorHandler(CPLQuietErrorHandler);
5342 : const bool bIsValid =
5343 20 : ValidateCutline(poMultiPolygon.get(), false);
5344 20 : CPLPopErrorHandler();
5345 20 : if (!bIsValid)
5346 : {
5347 1 : if (i == MAX_ITERATIONS - 1)
5348 : {
5349 0 : char *pszWKT = nullptr;
5350 0 : poMultiPolygon->exportToWkt(&pszWKT);
5351 0 : CPLDebug("WARP",
5352 : "After densification, cutline maximum "
5353 : "segment size is now %.0f pixel, "
5354 : "but cutline is invalid. %s",
5355 : dfMaxLengthInPixels, pszWKT);
5356 0 : CPLFree(pszWKT);
5357 0 : break;
5358 : }
5359 1 : CPLDebug("WARP",
5360 : "After densification, cutline maximum segment "
5361 : "size is now %.0f pixel, "
5362 : "but cutline is invalid. So trying a less "
5363 : "dense cutline.",
5364 : dfMaxLengthInPixels);
5365 1 : dfSegmentSize *= 2;
5366 1 : continue;
5367 : }
5368 : }
5369 21 : CPLDebug("WARP",
5370 : "After densification, cutline maximum segment size is "
5371 : "now %.0f pixel.",
5372 : dfMaxLengthInPixels);
5373 : }
5374 21 : break;
5375 : }
5376 : }
5377 :
5378 46 : if (eErr == OGRERR_FAILURE)
5379 : {
5380 0 : if (CPLTestBool(
5381 : CPLGetConfigOption("GDALWARP_IGNORE_BAD_CUTLINE", "NO")))
5382 0 : CPLError(CE_Warning, CPLE_AppDefined,
5383 : "Cutline transformation failed");
5384 : else
5385 : {
5386 0 : CPLError(CE_Failure, CPLE_AppDefined,
5387 : "Cutline transformation failed");
5388 0 : return CE_Failure;
5389 : }
5390 : }
5391 46 : else if (!ValidateCutline(poMultiPolygon.get(), true))
5392 : {
5393 1 : return CE_Failure;
5394 : }
5395 :
5396 : // Optimization: if the cutline contains the footprint of the source
5397 : // dataset, no need to use the cutline.
5398 45 : if (OGRGeometryFactory::haveGEOS()
5399 : #ifdef DEBUG
5400 : // Env var just for debugging purposes
5401 45 : && !CPLTestBool(CPLGetConfigOption(
5402 : "GDALWARP_SKIP_CUTLINE_CONTAINMENT_TEST", "NO"))
5403 : #endif
5404 : )
5405 : {
5406 44 : const double dfCutlineBlendDist = CPLAtof(CSLFetchNameValueDef(
5407 : *ppapszWarpOptions, "CUTLINE_BLEND_DIST", "0"));
5408 44 : auto poRing = std::make_unique<OGRLinearRing>();
5409 44 : poRing->addPoint(-dfCutlineBlendDist, -dfCutlineBlendDist);
5410 44 : poRing->addPoint(-dfCutlineBlendDist,
5411 44 : dfCutlineBlendDist + poSrcDS->GetRasterYSize());
5412 88 : poRing->addPoint(dfCutlineBlendDist + poSrcDS->GetRasterXSize(),
5413 44 : dfCutlineBlendDist + poSrcDS->GetRasterYSize());
5414 44 : poRing->addPoint(dfCutlineBlendDist + poSrcDS->GetRasterXSize(),
5415 : -dfCutlineBlendDist);
5416 44 : poRing->addPoint(-dfCutlineBlendDist, -dfCutlineBlendDist);
5417 44 : OGRPolygon oSrcDSFootprint;
5418 44 : oSrcDSFootprint.addRing(std::move(poRing));
5419 44 : OGREnvelope sSrcDSEnvelope;
5420 44 : oSrcDSFootprint.getEnvelope(&sSrcDSEnvelope);
5421 44 : OGREnvelope sCutlineEnvelope;
5422 44 : poMultiPolygon->getEnvelope(&sCutlineEnvelope);
5423 55 : if (sCutlineEnvelope.Contains(sSrcDSEnvelope) &&
5424 11 : poMultiPolygon->Contains(&oSrcDSFootprint))
5425 : {
5426 3 : CPLDebug("WARP", "Source dataset fully contained within cutline.");
5427 3 : return CE_None;
5428 : }
5429 : }
5430 :
5431 : /* -------------------------------------------------------------------- */
5432 : /* Convert aggregate geometry into WKT. */
5433 : /* -------------------------------------------------------------------- */
5434 42 : char *pszWKT = nullptr;
5435 42 : poMultiPolygon->exportToWkt(&pszWKT);
5436 : // fprintf(stderr, "WKT = \"%s\"\n", pszWKT ? pszWKT : "(null)");
5437 :
5438 42 : *ppapszWarpOptions = CSLSetNameValue(*ppapszWarpOptions, "CUTLINE", pszWKT);
5439 42 : CPLFree(pszWKT);
5440 42 : return CE_None;
5441 : }
5442 :
5443 47 : static void RemoveConflictingMetadata(GDALMajorObjectH hObj,
5444 : CSLConstList papszSrcMetadata,
5445 : const char *pszValueConflict)
5446 : {
5447 47 : if (hObj == nullptr)
5448 0 : return;
5449 :
5450 38 : for (const auto &[pszKey, pszValue] :
5451 85 : cpl::IterateNameValue(papszSrcMetadata))
5452 : {
5453 19 : const char *pszValueComp = GDALGetMetadataItem(hObj, pszKey, nullptr);
5454 19 : if (pszValueComp == nullptr || (!EQUAL(pszValue, pszValueComp) &&
5455 2 : !EQUAL(pszValueComp, pszValueConflict)))
5456 : {
5457 3 : if (STARTS_WITH(pszKey, "STATISTICS_"))
5458 1 : GDALSetMetadataItem(hObj, pszKey, nullptr, nullptr);
5459 : else
5460 2 : GDALSetMetadataItem(hObj, pszKey, pszValueConflict, nullptr);
5461 : }
5462 : }
5463 : }
5464 :
5465 : /************************************************************************/
5466 : /* IsValidSRS */
5467 : /************************************************************************/
5468 :
5469 331 : static bool IsValidSRS(const char *pszUserInput)
5470 :
5471 : {
5472 : OGRSpatialReferenceH hSRS;
5473 331 : bool bRes = true;
5474 :
5475 331 : hSRS = OSRNewSpatialReference(nullptr);
5476 331 : if (OSRSetFromUserInput(hSRS, pszUserInput) != OGRERR_NONE)
5477 : {
5478 0 : bRes = false;
5479 : }
5480 :
5481 331 : OSRDestroySpatialReference(hSRS);
5482 :
5483 331 : return bRes;
5484 : }
5485 :
5486 : /************************************************************************/
5487 : /* GDALWarpAppOptionsGetParser() */
5488 : /************************************************************************/
5489 :
5490 : static std::unique_ptr<GDALArgumentParser>
5491 1182 : GDALWarpAppOptionsGetParser(GDALWarpAppOptions *psOptions,
5492 : GDALWarpAppOptionsForBinary *psOptionsForBinary)
5493 : {
5494 : auto argParser = std::make_unique<GDALArgumentParser>(
5495 1182 : "gdalwarp", /* bForBinary=*/psOptionsForBinary != nullptr);
5496 :
5497 1182 : argParser->add_description(_("Image reprojection and warping utility."));
5498 :
5499 1182 : argParser->add_epilog(
5500 1182 : _("For more details, consult https://gdal.org/programs/gdalwarp.html"));
5501 :
5502 : argParser->add_quiet_argument(
5503 1182 : psOptionsForBinary ? &psOptionsForBinary->bQuiet : nullptr);
5504 :
5505 1182 : argParser->add_argument("-overwrite")
5506 1182 : .flag()
5507 : .action(
5508 90 : [psOptionsForBinary](const std::string &)
5509 : {
5510 49 : if (psOptionsForBinary)
5511 41 : psOptionsForBinary->bOverwrite = true;
5512 1182 : })
5513 1182 : .help(_("Overwrite the target dataset if it already exists."));
5514 :
5515 1182 : argParser->add_output_format_argument(psOptions->osFormat);
5516 :
5517 1182 : argParser->add_argument("-co")
5518 2364 : .metavar("<NAME>=<VALUE>")
5519 1182 : .append()
5520 : .action(
5521 305 : [psOptions, psOptionsForBinary](const std::string &s)
5522 : {
5523 149 : psOptions->aosCreateOptions.AddString(s.c_str());
5524 149 : psOptions->bCreateOutput = true;
5525 :
5526 149 : if (psOptionsForBinary)
5527 7 : psOptionsForBinary->aosCreateOptions.AddString(s.c_str());
5528 1182 : })
5529 1182 : .help(_("Creation option(s)."));
5530 :
5531 1182 : argParser->add_argument("-s_srs")
5532 2364 : .metavar("<srs_def>")
5533 : .action(
5534 64 : [psOptions](const std::string &s)
5535 : {
5536 32 : if (!IsValidSRS(s.c_str()))
5537 : {
5538 0 : throw std::invalid_argument("Invalid SRS for -s_srs");
5539 : }
5540 : psOptions->aosTransformerOptions.SetNameValue("SRC_SRS",
5541 32 : s.c_str());
5542 1214 : })
5543 1182 : .help(_("Set source spatial reference."));
5544 :
5545 1182 : argParser->add_argument("-t_srs")
5546 2364 : .metavar("<srs_def>")
5547 : .action(
5548 568 : [psOptions](const std::string &s)
5549 : {
5550 284 : if (!IsValidSRS(s.c_str()))
5551 : {
5552 0 : throw std::invalid_argument("Invalid SRS for -t_srs");
5553 : }
5554 : psOptions->aosTransformerOptions.SetNameValue("DST_SRS",
5555 284 : s.c_str());
5556 1466 : })
5557 1182 : .help(_("Set target spatial reference."));
5558 :
5559 : {
5560 1182 : auto &group = argParser->add_mutually_exclusive_group();
5561 1182 : group.add_argument("-srcalpha")
5562 1182 : .flag()
5563 1182 : .store_into(psOptions->bEnableSrcAlpha)
5564 : .help(_("Force the last band of a source image to be considered as "
5565 1182 : "a source alpha band."));
5566 1182 : group.add_argument("-nosrcalpha")
5567 1182 : .flag()
5568 1182 : .store_into(psOptions->bDisableSrcAlpha)
5569 : .help(_("Prevent the alpha band of a source image to be considered "
5570 1182 : "as such."));
5571 : }
5572 :
5573 1182 : argParser->add_argument("-dstalpha")
5574 1182 : .flag()
5575 1182 : .store_into(psOptions->bEnableDstAlpha)
5576 : .help(_("Create an output alpha band to identify nodata "
5577 1182 : "(unset/transparent) pixels."));
5578 :
5579 : // Parsing of that option is done in a preprocessing stage
5580 1182 : argParser->add_argument("-tr")
5581 2364 : .metavar("<xres> <yres>|square")
5582 1182 : .help(_("Target resolution."));
5583 :
5584 1182 : argParser->add_argument("-ts")
5585 2364 : .metavar("<width> <height>")
5586 1182 : .nargs(2)
5587 1182 : .scan<'i', int>()
5588 1182 : .help(_("Set output file size in pixels and lines."));
5589 :
5590 1182 : argParser->add_argument("-te")
5591 2364 : .metavar("<xmin> <ymin> <xmax> <ymax>")
5592 1182 : .nargs(4)
5593 1182 : .scan<'g', double>()
5594 1182 : .help(_("Set georeferenced extents of output file to be created."));
5595 :
5596 1182 : argParser->add_argument("-te_srs")
5597 2364 : .metavar("<srs_def>")
5598 : .action(
5599 18 : [psOptions](const std::string &s)
5600 : {
5601 6 : if (!IsValidSRS(s.c_str()))
5602 : {
5603 0 : throw std::invalid_argument("Invalid SRS for -te_srs");
5604 : }
5605 6 : psOptions->osTE_SRS = s;
5606 6 : psOptions->bCreateOutput = true;
5607 1188 : })
5608 1182 : .help(_("Set source spatial reference."));
5609 :
5610 1182 : argParser->add_argument("-r")
5611 : .metavar("near|bilinear|cubic|cubicspline|lanczos|average|rms|mode|min|"
5612 2364 : "max|med|q1|q3|sum")
5613 : .action(
5614 1399 : [psOptions](const std::string &s)
5615 : {
5616 700 : GDALGetWarpResampleAlg(s.c_str(), psOptions->eResampleAlg,
5617 : /*bThrow=*/true);
5618 699 : psOptions->bResampleAlgSpecifiedByUser = true;
5619 1182 : })
5620 1182 : .help(_("Resampling method to use."));
5621 :
5622 1182 : argParser->add_output_type_argument(psOptions->eOutputType);
5623 :
5624 : ///////////////////////////////////////////////////////////////////////
5625 1182 : argParser->add_group("Advanced options");
5626 :
5627 32 : const auto CheckSingleMethod = [psOptions]()
5628 : {
5629 : const char *pszMethod =
5630 16 : FetchSrcMethod(psOptions->aosTransformerOptions);
5631 16 : if (pszMethod)
5632 0 : CPLError(CE_Warning, CPLE_IllegalArg,
5633 : "Warning: only one METHOD can be used. Method %s is "
5634 : "already defined.",
5635 : pszMethod);
5636 : const char *pszMAX_GCP_ORDER =
5637 16 : psOptions->aosTransformerOptions.FetchNameValue("MAX_GCP_ORDER");
5638 16 : if (pszMAX_GCP_ORDER)
5639 0 : CPLError(CE_Warning, CPLE_IllegalArg,
5640 : "Warning: only one METHOD can be used. -order %s "
5641 : "option was specified, so it is likely that "
5642 : "GCP_POLYNOMIAL was implied.",
5643 : pszMAX_GCP_ORDER);
5644 1198 : };
5645 :
5646 1182 : argParser->add_argument("-wo")
5647 2364 : .metavar("<NAME>=<VALUE>")
5648 1182 : .append()
5649 229 : .action([psOptions](const std::string &s)
5650 1411 : { psOptions->aosWarpOptions.AddString(s.c_str()); })
5651 1182 : .help(_("Warping option(s)."));
5652 :
5653 1182 : argParser->add_argument("-multi")
5654 1182 : .flag()
5655 1182 : .store_into(psOptions->bMulti)
5656 1182 : .help(_("Multithreaded input/output."));
5657 :
5658 1182 : argParser->add_argument("-s_coord_epoch")
5659 2364 : .metavar("<epoch>")
5660 : .action(
5661 0 : [psOptions](const std::string &s)
5662 : {
5663 : psOptions->aosTransformerOptions.SetNameValue(
5664 0 : "SRC_COORDINATE_EPOCH", s.c_str());
5665 1182 : })
5666 : .help(_("Assign a coordinate epoch, linked with the source SRS when "
5667 1182 : "-s_srs is used."));
5668 :
5669 1182 : argParser->add_argument("-t_coord_epoch")
5670 2364 : .metavar("<epoch>")
5671 : .action(
5672 0 : [psOptions](const std::string &s)
5673 : {
5674 : psOptions->aosTransformerOptions.SetNameValue(
5675 0 : "DST_COORDINATE_EPOCH", s.c_str());
5676 1182 : })
5677 : .help(_("Assign a coordinate epoch, linked with the output SRS when "
5678 1182 : "-t_srs is used."));
5679 :
5680 1182 : argParser->add_argument("-ct")
5681 2364 : .metavar("<string>")
5682 : .action(
5683 4 : [psOptions](const std::string &s)
5684 : {
5685 : psOptions->aosTransformerOptions.SetNameValue(
5686 4 : "COORDINATE_OPERATION", s.c_str());
5687 1182 : })
5688 1182 : .help(_("Set a coordinate transformation."));
5689 :
5690 : {
5691 1182 : auto &group = argParser->add_mutually_exclusive_group();
5692 1182 : group.add_argument("-tps")
5693 1182 : .flag()
5694 : .action(
5695 10 : [psOptions, CheckSingleMethod](const std::string &)
5696 : {
5697 5 : CheckSingleMethod();
5698 : psOptions->aosTransformerOptions.SetNameValue("SRC_METHOD",
5699 5 : "GCP_TPS");
5700 1182 : })
5701 : .help(_("Force use of thin plate spline transformer based on "
5702 1182 : "available GCPs."));
5703 :
5704 1182 : group.add_argument("-rpc")
5705 1182 : .flag()
5706 : .action(
5707 4 : [psOptions, CheckSingleMethod](const std::string &)
5708 : {
5709 2 : CheckSingleMethod();
5710 : psOptions->aosTransformerOptions.SetNameValue("SRC_METHOD",
5711 2 : "RPC");
5712 1182 : })
5713 1182 : .help(_("Force use of RPCs."));
5714 :
5715 1182 : group.add_argument("-geoloc")
5716 1182 : .flag()
5717 : .action(
5718 18 : [psOptions, CheckSingleMethod](const std::string &)
5719 : {
5720 9 : CheckSingleMethod();
5721 : psOptions->aosTransformerOptions.SetNameValue(
5722 9 : "SRC_METHOD", "GEOLOC_ARRAY");
5723 1182 : })
5724 1182 : .help(_("Force use of Geolocation Arrays."));
5725 : }
5726 :
5727 1182 : argParser->add_argument("-order")
5728 2364 : .metavar("<1|2|3>")
5729 1182 : .choices("1", "2", "3")
5730 : .action(
5731 0 : [psOptions](const std::string &s)
5732 : {
5733 : const char *pszMethod =
5734 0 : FetchSrcMethod(psOptions->aosTransformerOptions);
5735 0 : if (pszMethod)
5736 0 : CPLError(
5737 : CE_Warning, CPLE_IllegalArg,
5738 : "Warning: only one METHOD can be used. Method %s is "
5739 : "already defined",
5740 : pszMethod);
5741 : psOptions->aosTransformerOptions.SetNameValue("MAX_GCP_ORDER",
5742 0 : s.c_str());
5743 1182 : })
5744 1182 : .help(_("Order of polynomial used for GCP warping."));
5745 :
5746 : // Parsing of that option is done in a preprocessing stage
5747 1182 : argParser->add_argument("-refine_gcps")
5748 2364 : .metavar("<tolerance> [<minimum_gcps>]")
5749 1182 : .help(_("Refines the GCPs by automatically eliminating outliers."));
5750 :
5751 1182 : argParser->add_argument("-to")
5752 2364 : .metavar("<NAME>=<VALUE>")
5753 1182 : .append()
5754 64 : .action([psOptions](const std::string &s)
5755 1246 : { psOptions->aosTransformerOptions.AddString(s.c_str()); })
5756 1182 : .help(_("Transform option(s)."));
5757 :
5758 1182 : argParser->add_argument("-et")
5759 2364 : .metavar("<err_threshold>")
5760 1182 : .store_into(psOptions->dfErrorThreshold)
5761 : .action(
5762 26 : [psOptions](const std::string &)
5763 : {
5764 13 : if (psOptions->dfErrorThreshold < 0)
5765 : {
5766 : throw std::invalid_argument(
5767 0 : "Invalid value for error threshold");
5768 : }
5769 : psOptions->aosWarpOptions.AddString(CPLSPrintf(
5770 13 : "ERROR_THRESHOLD=%.16g", psOptions->dfErrorThreshold));
5771 1195 : })
5772 1182 : .help(_("Error threshold."));
5773 :
5774 1182 : argParser->add_argument("-wm")
5775 2364 : .metavar("<memory_in_mb>")
5776 : .action(
5777 74 : [psOptions](const std::string &s)
5778 : {
5779 38 : bool bUnitSpecified = false;
5780 : GIntBig nBytes;
5781 38 : if (CPLParseMemorySize(s.c_str(), &nBytes, &bUnitSpecified) ==
5782 : CE_None)
5783 : {
5784 36 : if (!bUnitSpecified && nBytes < 10000)
5785 : {
5786 7 : nBytes *= (1024 * 1024);
5787 : }
5788 36 : psOptions->dfWarpMemoryLimit = static_cast<double>(nBytes);
5789 : }
5790 : else
5791 : {
5792 2 : throw std::invalid_argument("Failed to parse value of -wm");
5793 : }
5794 1218 : })
5795 1182 : .help(_("Set max warp memory."));
5796 :
5797 1182 : argParser->add_argument("-srcnodata")
5798 2364 : .metavar("\"<value>[ <value>]...\"")
5799 1182 : .store_into(psOptions->osSrcNodata)
5800 1182 : .help(_("Nodata masking values for input bands."));
5801 :
5802 1182 : argParser->add_argument("-dstnodata")
5803 2364 : .metavar("\"<value>[ <value>]...\"")
5804 1182 : .store_into(psOptions->osDstNodata)
5805 1182 : .help(_("Nodata masking values for output bands."));
5806 :
5807 1182 : argParser->add_argument("-tap")
5808 1182 : .flag()
5809 1182 : .store_into(psOptions->bTargetAlignedPixels)
5810 1182 : .help(_("Force target aligned pixels."));
5811 :
5812 1182 : argParser->add_argument("-wt")
5813 2364 : .metavar("Byte|Int8|[U]Int{16|32|64}|CInt{16|32}|[C]Float{32|64}")
5814 : .action(
5815 0 : [psOptions](const std::string &s)
5816 : {
5817 0 : psOptions->eWorkingType = GDALGetDataTypeByName(s.c_str());
5818 0 : if (psOptions->eWorkingType == GDT_Unknown)
5819 : {
5820 : throw std::invalid_argument(
5821 0 : std::string("Unknown output pixel type: ").append(s));
5822 : }
5823 1182 : })
5824 1182 : .help(_("Working data type."));
5825 :
5826 : // Non-documented alias of -r nearest
5827 1182 : argParser->add_argument("-rn")
5828 1182 : .flag()
5829 1182 : .hidden()
5830 1 : .action([psOptions](const std::string &)
5831 1182 : { psOptions->eResampleAlg = GRA_NearestNeighbour; })
5832 1182 : .help(_("Nearest neighbour resampling."));
5833 :
5834 : // Non-documented alias of -r bilinear
5835 1182 : argParser->add_argument("-rb")
5836 1182 : .flag()
5837 1182 : .hidden()
5838 2 : .action([psOptions](const std::string &)
5839 1182 : { psOptions->eResampleAlg = GRA_Bilinear; })
5840 1182 : .help(_("Bilinear resampling."));
5841 :
5842 : // Non-documented alias of -r cubic
5843 1182 : argParser->add_argument("-rc")
5844 1182 : .flag()
5845 1182 : .hidden()
5846 1 : .action([psOptions](const std::string &)
5847 1182 : { psOptions->eResampleAlg = GRA_Cubic; })
5848 1182 : .help(_("Cubic resampling."));
5849 :
5850 : // Non-documented alias of -r cubicspline
5851 1182 : argParser->add_argument("-rcs")
5852 1182 : .flag()
5853 1182 : .hidden()
5854 1 : .action([psOptions](const std::string &)
5855 1182 : { psOptions->eResampleAlg = GRA_CubicSpline; })
5856 1182 : .help(_("Cubic spline resampling."));
5857 :
5858 : // Non-documented alias of -r lanczos
5859 1182 : argParser->add_argument("-rl")
5860 1182 : .flag()
5861 1182 : .hidden()
5862 0 : .action([psOptions](const std::string &)
5863 1182 : { psOptions->eResampleAlg = GRA_Lanczos; })
5864 1182 : .help(_("Lanczos resampling."));
5865 :
5866 : // Non-documented alias of -r average
5867 1182 : argParser->add_argument("-ra")
5868 1182 : .flag()
5869 1182 : .hidden()
5870 0 : .action([psOptions](const std::string &)
5871 1182 : { psOptions->eResampleAlg = GRA_Average; })
5872 1182 : .help(_("Average resampling."));
5873 :
5874 : // Non-documented alias of -r rms
5875 1182 : argParser->add_argument("-rrms")
5876 1182 : .flag()
5877 1182 : .hidden()
5878 0 : .action([psOptions](const std::string &)
5879 1182 : { psOptions->eResampleAlg = GRA_RMS; })
5880 1182 : .help(_("RMS resampling."));
5881 :
5882 : // Non-documented alias of -r mode
5883 1182 : argParser->add_argument("-rm")
5884 1182 : .flag()
5885 1182 : .hidden()
5886 0 : .action([psOptions](const std::string &)
5887 1182 : { psOptions->eResampleAlg = GRA_Mode; })
5888 1182 : .help(_("Mode resampling."));
5889 :
5890 1182 : argParser->add_argument("-cutline")
5891 2364 : .metavar("<datasource>|<WKT>")
5892 1182 : .store_into(psOptions->osCutlineDSNameOrWKT)
5893 : .help(_("Enable use of a blend cutline from the name of a vector "
5894 1182 : "dataset or a WKT geometry."));
5895 :
5896 1182 : argParser->add_argument("-cutline_srs")
5897 2364 : .metavar("<srs_def>")
5898 : .action(
5899 18 : [psOptions](const std::string &s)
5900 : {
5901 9 : if (!IsValidSRS(s.c_str()))
5902 : {
5903 0 : throw std::invalid_argument("Invalid SRS for -cutline_srs");
5904 : }
5905 9 : psOptions->osCutlineSRS = s;
5906 1191 : })
5907 1182 : .help(_("Sets/overrides cutline SRS."));
5908 :
5909 1182 : argParser->add_argument("-cwhere")
5910 2364 : .metavar("<expression>")
5911 1182 : .store_into(psOptions->osCWHERE)
5912 1182 : .help(_("Restrict desired cutline features based on attribute query."));
5913 :
5914 : {
5915 1182 : auto &group = argParser->add_mutually_exclusive_group();
5916 1182 : group.add_argument("-cl")
5917 2364 : .metavar("<layername>")
5918 1182 : .store_into(psOptions->osCLayer)
5919 1182 : .help(_("Select the named layer from the cutline datasource."));
5920 :
5921 1182 : group.add_argument("-csql")
5922 2364 : .metavar("<query>")
5923 1182 : .store_into(psOptions->osCSQL)
5924 1182 : .help(_("Select cutline features using an SQL query."));
5925 : }
5926 :
5927 1182 : argParser->add_argument("-cblend")
5928 2364 : .metavar("<distance>")
5929 : .action(
5930 0 : [psOptions](const std::string &s)
5931 : {
5932 : psOptions->aosWarpOptions.SetNameValue("CUTLINE_BLEND_DIST",
5933 0 : s.c_str());
5934 1182 : })
5935 : .help(_(
5936 1182 : "Set a blend distance to use to blend over cutlines (in pixels)."));
5937 :
5938 1182 : argParser->add_argument("-crop_to_cutline")
5939 1182 : .flag()
5940 : .action(
5941 18 : [psOptions](const std::string &)
5942 : {
5943 18 : psOptions->bCropToCutline = true;
5944 18 : psOptions->bCreateOutput = true;
5945 1182 : })
5946 : .help(_("Crop the extent of the target dataset to the extent of the "
5947 1182 : "cutline."));
5948 :
5949 1182 : argParser->add_argument("-nomd")
5950 1182 : .flag()
5951 : .action(
5952 0 : [psOptions](const std::string &)
5953 : {
5954 0 : psOptions->bCopyMetadata = false;
5955 0 : psOptions->bCopyBandInfo = false;
5956 1182 : })
5957 1182 : .help(_("Do not copy metadata."));
5958 :
5959 1182 : argParser->add_argument("-cvmd")
5960 2364 : .metavar("<meta_conflict_value>")
5961 1182 : .store_into(psOptions->osMDConflictValue)
5962 : .help(_("Value to set metadata items that conflict between source "
5963 1182 : "datasets."));
5964 :
5965 1182 : argParser->add_argument("-setci")
5966 1182 : .flag()
5967 1182 : .store_into(psOptions->bSetColorInterpretation)
5968 : .help(_("Set the color interpretation of the bands of the target "
5969 1182 : "dataset from the source dataset."));
5970 :
5971 : argParser->add_open_options_argument(
5972 1182 : psOptionsForBinary ? &(psOptionsForBinary->aosOpenOptions) : nullptr);
5973 :
5974 1182 : argParser->add_argument("-doo")
5975 2364 : .metavar("<NAME>=<VALUE>")
5976 1182 : .append()
5977 : .action(
5978 0 : [psOptionsForBinary](const std::string &s)
5979 : {
5980 0 : if (psOptionsForBinary)
5981 0 : psOptionsForBinary->aosDestOpenOptions.AddString(s.c_str());
5982 1182 : })
5983 1182 : .help(_("Open option(s) for output dataset."));
5984 :
5985 1182 : argParser->add_argument("-ovr")
5986 2364 : .metavar("<level>|AUTO|AUTO-<n>|NONE")
5987 : .action(
5988 24 : [psOptions](const std::string &s)
5989 : {
5990 12 : const char *pszOvLevel = s.c_str();
5991 12 : if (EQUAL(pszOvLevel, "AUTO"))
5992 1 : psOptions->nOvLevel = OVR_LEVEL_AUTO;
5993 11 : else if (STARTS_WITH_CI(pszOvLevel, "AUTO-"))
5994 1 : psOptions->nOvLevel =
5995 1 : OVR_LEVEL_AUTO - atoi(pszOvLevel + strlen("AUTO-"));
5996 10 : else if (EQUAL(pszOvLevel, "NONE"))
5997 5 : psOptions->nOvLevel = OVR_LEVEL_NONE;
5998 5 : else if (CPLGetValueType(pszOvLevel) == CPL_VALUE_INTEGER)
5999 5 : psOptions->nOvLevel = atoi(pszOvLevel);
6000 : else
6001 : {
6002 : throw std::invalid_argument(CPLSPrintf(
6003 0 : "Invalid value '%s' for -ov option", pszOvLevel));
6004 : }
6005 1194 : })
6006 1182 : .help(_("Specify which overview level of source files must be used."));
6007 :
6008 : {
6009 1182 : auto &group = argParser->add_mutually_exclusive_group();
6010 1182 : group.add_argument("-vshift")
6011 1182 : .flag()
6012 1182 : .store_into(psOptions->bVShift)
6013 1182 : .help(_("Force the use of vertical shift."));
6014 1182 : group.add_argument("-novshift", "-novshiftgrid")
6015 1182 : .flag()
6016 1182 : .store_into(psOptions->bNoVShift)
6017 1182 : .help(_("Disable the use of vertical shift."));
6018 : }
6019 :
6020 : argParser->add_input_format_argument(
6021 : psOptionsForBinary ? &psOptionsForBinary->aosAllowedInputDrivers
6022 1182 : : nullptr);
6023 :
6024 1182 : argParser->add_argument("-b", "-srcband")
6025 2364 : .metavar("<band>")
6026 1182 : .append()
6027 1182 : .store_into(psOptions->anSrcBands)
6028 1182 : .help(_("Specify input band(s) number to warp."));
6029 :
6030 1182 : argParser->add_argument("-dstband")
6031 2364 : .metavar("<band>")
6032 1182 : .append()
6033 1182 : .store_into(psOptions->anDstBands)
6034 1182 : .help(_("Specify the output band number in which to warp."));
6035 :
6036 1182 : if (psOptionsForBinary)
6037 : {
6038 92 : argParser->add_argument("src_dataset_name")
6039 184 : .metavar("<src_dataset_name>")
6040 92 : .nargs(argparse::nargs_pattern::at_least_one)
6041 93 : .action([psOptionsForBinary](const std::string &s)
6042 185 : { psOptionsForBinary->aosSrcFiles.AddString(s.c_str()); })
6043 92 : .help(_("Input dataset(s)."));
6044 :
6045 92 : argParser->add_argument("dst_dataset_name")
6046 184 : .metavar("<dst_dataset_name>")
6047 92 : .store_into(psOptionsForBinary->osDstFilename)
6048 92 : .help(_("Output dataset."));
6049 : }
6050 :
6051 2364 : return argParser;
6052 : }
6053 :
6054 : /************************************************************************/
6055 : /* GDALWarpAppGetParserUsage() */
6056 : /************************************************************************/
6057 :
6058 2 : std::string GDALWarpAppGetParserUsage()
6059 : {
6060 : try
6061 : {
6062 4 : GDALWarpAppOptions sOptions;
6063 4 : GDALWarpAppOptionsForBinary sOptionsForBinary;
6064 : auto argParser =
6065 4 : GDALWarpAppOptionsGetParser(&sOptions, &sOptionsForBinary);
6066 2 : return argParser->usage();
6067 : }
6068 0 : catch (const std::exception &err)
6069 : {
6070 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
6071 0 : err.what());
6072 0 : return std::string();
6073 : }
6074 : }
6075 :
6076 : /************************************************************************/
6077 : /* GDALWarpAppOptionsNew() */
6078 : /************************************************************************/
6079 :
6080 : #ifndef CheckHasEnoughAdditionalArgs_defined
6081 : #define CheckHasEnoughAdditionalArgs_defined
6082 :
6083 120 : static bool CheckHasEnoughAdditionalArgs(CSLConstList papszArgv, int i,
6084 : int nExtraArg, int nArgc)
6085 : {
6086 120 : if (i + nExtraArg >= nArgc)
6087 : {
6088 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6089 0 : "%s option requires %d argument%s", papszArgv[i], nExtraArg,
6090 : nExtraArg == 1 ? "" : "s");
6091 0 : return false;
6092 : }
6093 120 : return true;
6094 : }
6095 : #endif
6096 :
6097 : #define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
6098 : if (!CheckHasEnoughAdditionalArgs(papszArgv, i, nExtraArg, nArgc)) \
6099 : { \
6100 : return nullptr; \
6101 : }
6102 :
6103 : /**
6104 : * Allocates a GDALWarpAppOptions struct.
6105 : *
6106 : * @param papszArgv NULL terminated list of options (potentially including
6107 : * filename and open options too), or NULL. The accepted options are the ones of
6108 : * the <a href="/programs/gdalwarp.html">gdalwarp</a> utility.
6109 : * @param psOptionsForBinary (output) may be NULL (and should generally be
6110 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
6111 : * GDALWarpAppOptionsForBinaryNew() prior to this
6112 : * function. Will be filled with potentially present filename, open options,...
6113 : * @return pointer to the allocated GDALWarpAppOptions struct. Must be freed
6114 : * with GDALWarpAppOptionsFree().
6115 : *
6116 : * @since GDAL 2.1
6117 : */
6118 :
6119 : GDALWarpAppOptions *
6120 1180 : GDALWarpAppOptionsNew(char **papszArgv,
6121 : GDALWarpAppOptionsForBinary *psOptionsForBinary)
6122 : {
6123 2360 : auto psOptions = std::make_unique<GDALWarpAppOptions>();
6124 :
6125 1180 : psOptions->aosArgs.Assign(CSLDuplicate(papszArgv), true);
6126 :
6127 : /* -------------------------------------------------------------------- */
6128 : /* Pre-processing for custom syntax that ArgumentParser does not */
6129 : /* support. */
6130 : /* -------------------------------------------------------------------- */
6131 :
6132 2360 : CPLStringList aosArgv;
6133 1180 : const int nArgc = CSLCount(papszArgv);
6134 8095 : for (int i = 0;
6135 8095 : i < nArgc && papszArgv != nullptr && papszArgv[i] != nullptr; i++)
6136 : {
6137 6915 : if (EQUAL(papszArgv[i], "-refine_gcps"))
6138 : {
6139 0 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
6140 0 : psOptions->aosTransformerOptions.SetNameValue("REFINE_TOLERANCE",
6141 0 : papszArgv[++i]);
6142 0 : if (CPLAtof(papszArgv[i]) < 0)
6143 : {
6144 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6145 : "The tolerance for -refine_gcps may not be negative.");
6146 0 : return nullptr;
6147 : }
6148 0 : if (i < nArgc - 1 && atoi(papszArgv[i + 1]) >= 0 &&
6149 0 : isdigit(static_cast<unsigned char>(papszArgv[i + 1][0])))
6150 : {
6151 0 : psOptions->aosTransformerOptions.SetNameValue(
6152 0 : "REFINE_MINIMUM_GCPS", papszArgv[++i]);
6153 : }
6154 : else
6155 : {
6156 0 : psOptions->aosTransformerOptions.SetNameValue(
6157 0 : "REFINE_MINIMUM_GCPS", "-1");
6158 : }
6159 : }
6160 6915 : else if (EQUAL(papszArgv[i], "-tr") && i + 1 < nArgc &&
6161 121 : EQUAL(papszArgv[i + 1], "square"))
6162 : {
6163 1 : ++i;
6164 1 : psOptions->bSquarePixels = true;
6165 1 : psOptions->bCreateOutput = true;
6166 : }
6167 6914 : else if (EQUAL(papszArgv[i], "-tr"))
6168 : {
6169 120 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
6170 120 : psOptions->dfXRes = CPLAtofM(papszArgv[++i]);
6171 120 : psOptions->dfYRes = fabs(CPLAtofM(papszArgv[++i]));
6172 120 : if (psOptions->dfXRes == 0 || psOptions->dfYRes == 0)
6173 : {
6174 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6175 : "Wrong value for -tr parameters.");
6176 0 : return nullptr;
6177 : }
6178 120 : psOptions->bCreateOutput = true;
6179 : }
6180 : // argparser will be confused if the value of a string argument
6181 : // starts with a negative sign.
6182 6794 : else if (EQUAL(papszArgv[i], "-srcnodata") && i + 1 < nArgc)
6183 : {
6184 48 : ++i;
6185 48 : psOptions->osSrcNodata = papszArgv[i];
6186 : }
6187 : // argparser will be confused if the value of a string argument
6188 : // starts with a negative sign.
6189 6746 : else if (EQUAL(papszArgv[i], "-dstnodata") && i + 1 < nArgc)
6190 : {
6191 140 : ++i;
6192 140 : psOptions->osDstNodata = papszArgv[i];
6193 : }
6194 : else
6195 : {
6196 6606 : aosArgv.AddString(papszArgv[i]);
6197 : }
6198 : }
6199 :
6200 : try
6201 : {
6202 : auto argParser =
6203 2360 : GDALWarpAppOptionsGetParser(psOptions.get(), psOptionsForBinary);
6204 :
6205 1180 : argParser->parse_args_without_binary_name(aosArgv.List());
6206 :
6207 1350 : if (auto oTS = argParser->present<std::vector<int>>("-ts"))
6208 : {
6209 174 : psOptions->nForcePixels = (*oTS)[0];
6210 174 : psOptions->nForceLines = (*oTS)[1];
6211 174 : psOptions->bCreateOutput = true;
6212 : }
6213 :
6214 1367 : if (auto oTE = argParser->present<std::vector<double>>("-te"))
6215 : {
6216 191 : psOptions->dfMinX = (*oTE)[0];
6217 191 : psOptions->dfMinY = (*oTE)[1];
6218 191 : psOptions->dfMaxX = (*oTE)[2];
6219 191 : psOptions->dfMaxY = (*oTE)[3];
6220 191 : psOptions->aosTransformerOptions.SetNameValue(
6221 : "TARGET_EXTENT",
6222 191 : CPLSPrintf("%.17g,%.17g,%.17g,%.17g", psOptions->dfMinX,
6223 382 : psOptions->dfMinY, psOptions->dfMaxX,
6224 191 : psOptions->dfMaxY));
6225 191 : psOptions->bCreateOutput = true;
6226 : }
6227 :
6228 1181 : if (!psOptions->anDstBands.empty() &&
6229 5 : psOptions->anSrcBands.size() != psOptions->anDstBands.size())
6230 : {
6231 1 : CPLError(
6232 : CE_Failure, CPLE_IllegalArg,
6233 : "-srcband should be specified as many times as -dstband is");
6234 1 : return nullptr;
6235 : }
6236 1202 : else if (!psOptions->anSrcBands.empty() &&
6237 27 : psOptions->anDstBands.empty())
6238 : {
6239 60 : for (int i = 0; i < static_cast<int>(psOptions->anSrcBands.size());
6240 : ++i)
6241 : {
6242 37 : psOptions->anDstBands.push_back(i + 1);
6243 : }
6244 : }
6245 :
6246 1686 : if (!psOptions->osFormat.empty() ||
6247 511 : psOptions->eOutputType != GDT_Unknown)
6248 : {
6249 670 : psOptions->bCreateOutput = true;
6250 : }
6251 :
6252 1175 : if (psOptionsForBinary)
6253 88 : psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
6254 :
6255 1175 : return psOptions.release();
6256 : }
6257 4 : catch (const std::exception &err)
6258 : {
6259 4 : CPLError(CE_Failure, CPLE_AppDefined, "%s", err.what());
6260 4 : return nullptr;
6261 : }
6262 : }
6263 :
6264 : /************************************************************************/
6265 : /* GDALWarpAppOptionsFree() */
6266 : /************************************************************************/
6267 :
6268 : /**
6269 : * Frees the GDALWarpAppOptions struct.
6270 : *
6271 : * @param psOptions the options struct for GDALWarp().
6272 : *
6273 : * @since GDAL 2.1
6274 : */
6275 :
6276 1175 : void GDALWarpAppOptionsFree(GDALWarpAppOptions *psOptions)
6277 : {
6278 1175 : delete psOptions;
6279 1175 : }
6280 :
6281 : /************************************************************************/
6282 : /* GDALWarpAppOptionsSetProgress() */
6283 : /************************************************************************/
6284 :
6285 : /**
6286 : * Set a progress function.
6287 : *
6288 : * @param psOptions the options struct for GDALWarp().
6289 : * @param pfnProgress the progress callback.
6290 : * @param pProgressData the user data for the progress callback.
6291 : *
6292 : * @since GDAL 2.1
6293 : */
6294 :
6295 136 : void GDALWarpAppOptionsSetProgress(GDALWarpAppOptions *psOptions,
6296 : GDALProgressFunc pfnProgress,
6297 : void *pProgressData)
6298 : {
6299 136 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
6300 136 : psOptions->pProgressData = pProgressData;
6301 136 : if (pfnProgress == GDALTermProgress)
6302 0 : psOptions->bQuiet = false;
6303 136 : }
6304 :
6305 : /************************************************************************/
6306 : /* GDALWarpAppOptionsSetQuiet() */
6307 : /************************************************************************/
6308 :
6309 : /**
6310 : * Set a progress function.
6311 : *
6312 : * @param psOptions the options struct for GDALWarp().
6313 : * @param bQuiet whether GDALWarp() should emit messages on stdout.
6314 : *
6315 : * @since GDAL 2.3
6316 : */
6317 :
6318 81 : void GDALWarpAppOptionsSetQuiet(GDALWarpAppOptions *psOptions, int bQuiet)
6319 : {
6320 81 : psOptions->bQuiet = CPL_TO_BOOL(bQuiet);
6321 81 : }
6322 :
6323 : /************************************************************************/
6324 : /* GDALWarpAppOptionsSetWarpOption() */
6325 : /************************************************************************/
6326 :
6327 : /**
6328 : * Set a warp option
6329 : *
6330 : * @param psOptions the options struct for GDALWarp().
6331 : * @param pszKey key.
6332 : * @param pszValue value.
6333 : *
6334 : * @since GDAL 2.1
6335 : */
6336 :
6337 0 : void GDALWarpAppOptionsSetWarpOption(GDALWarpAppOptions *psOptions,
6338 : const char *pszKey, const char *pszValue)
6339 : {
6340 0 : psOptions->aosWarpOptions.SetNameValue(pszKey, pszValue);
6341 0 : }
6342 :
6343 : #undef CHECK_HAS_ENOUGH_ADDITIONAL_ARGS
|
