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