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