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 963 : static CPLString GetSrcDSProjection(GDALDatasetH hDS, CSLConstList papszTO)
329 : {
330 963 : const char *pszProjection = CSLFetchNameValue(papszTO, "SRC_SRS");
331 963 : if (pszProjection != nullptr || hDS == nullptr)
332 : {
333 59 : return pszProjection ? pszProjection : "";
334 : }
335 :
336 904 : const char *pszMethod = CSLFetchNameValue(papszTO, "METHOD");
337 904 : char **papszMD = nullptr;
338 904 : const OGRSpatialReferenceH hSRS = GDALGetSpatialRef(hDS);
339 : const char *pszGeolocationDataset =
340 904 : CSLFetchNameValueDef(papszTO, "SRC_GEOLOC_ARRAY",
341 : CSLFetchNameValue(papszTO, "GEOLOC_ARRAY"));
342 904 : 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 903 : 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 146 : else if (GDALGetGCPProjection(hDS) != nullptr &&
370 146 : strlen(GDALGetGCPProjection(hDS)) > 0 &&
371 296 : GDALGetGCPCount(hDS) > 1 &&
372 4 : (pszMethod == nullptr || STARTS_WITH_CI(pszMethod, "GCP_")))
373 : {
374 26 : pszProjection = GDALGetGCPProjection(hDS);
375 : }
376 123 : else if (GDALGetMetadata(hDS, "RPC") != nullptr &&
377 3 : (pszMethod == nullptr || EQUAL(pszMethod, "RPC")))
378 : {
379 14 : pszProjection = SRS_WKT_WGS84_LAT_LONG;
380 : }
381 115 : else if ((papszMD = GDALGetMetadata(hDS, "GEOLOCATION")) != nullptr &&
382 9 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")))
383 : {
384 16 : pszProjection = CSLFetchNameValue(papszMD, "SRS");
385 : }
386 146 : 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 1771 : static bool MustApplyVerticalShift(GDALDatasetH hWrkSrcDS,
622 : const GDALWarpAppOptions *psOptions,
623 : OGRSpatialReference &oSRSSrc,
624 : OGRSpatialReference &oSRSDst,
625 : bool &bSrcHasVertAxis, bool &bDstHasVertAxis)
626 : {
627 1771 : bool bApplyVShift = psOptions->bVShift;
628 :
629 : // Check if we must do vertical shift grid transform
630 : const char *pszSrcWKT =
631 1771 : psOptions->aosTransformerOptions.FetchNameValue("SRC_SRS");
632 1771 : if (pszSrcWKT)
633 76 : oSRSSrc.SetFromUserInput(pszSrcWKT);
634 : else
635 : {
636 1695 : auto hSRS = GDALGetSpatialRef(hWrkSrcDS);
637 1695 : if (hSRS)
638 1329 : oSRSSrc = *(OGRSpatialReference::FromHandle(hSRS));
639 : else
640 366 : 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 885 : static bool ApplyVerticalShift(GDALDatasetH hWrkSrcDS,
674 : const GDALWarpAppOptions *psOptions,
675 : GDALWarpOptions *psWO)
676 : {
677 885 : if (psOptions->bVShift)
678 : {
679 0 : psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,
680 : "APPLY_VERTICAL_SHIFT", "YES");
681 : }
682 :
683 1770 : OGRSpatialReference oSRSSrc;
684 885 : OGRSpatialReference oSRSDst;
685 885 : bool bSrcHasVertAxis = false;
686 885 : bool bDstHasVertAxis = false;
687 : bool bApplyVShift =
688 885 : MustApplyVerticalShift(hWrkSrcDS, psOptions, oSRSSrc, oSRSDst,
689 : bSrcHasVertAxis, bDstHasVertAxis);
690 :
691 1590 : if ((GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&
692 705 : (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 1770 : 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 897 : GDALDatasetH GDALWarp(const char *pszDest, GDALDatasetH hDstDS, int nSrcCount,
1387 : GDALDatasetH *pahSrcDS,
1388 : const GDALWarpAppOptions *psOptionsIn, int *pbUsageError)
1389 : {
1390 897 : CPLErrorReset();
1391 :
1392 1815 : for (int i = 0; i < nSrcCount; i++)
1393 : {
1394 918 : if (!pahSrcDS[i])
1395 0 : return nullptr;
1396 : }
1397 :
1398 1794 : GDALWarpAppOptions oOptionsTmp;
1399 897 : if (psOptionsIn)
1400 895 : oOptionsTmp = *psOptionsIn;
1401 897 : GDALWarpAppOptions *psOptions = &oOptionsTmp;
1402 :
1403 897 : if (hDstDS == nullptr)
1404 : {
1405 808 : 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 808 : auto hDriver = GDALGetDriverByName(psOptions->osFormat.c_str());
1416 1616 : if (hDriver != nullptr &&
1417 808 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) ==
1418 1616 : 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 888 : auto ret = GDALWarpDirect(pszDest, hDstDS, nSrcCount, pahSrcDS, psOptions,
1429 : pbUsageError);
1430 :
1431 888 : return ret;
1432 : }
1433 :
1434 : /************************************************************************/
1435 : /* UseTEAndTSAndTRConsistently() */
1436 : /************************************************************************/
1437 :
1438 810 : 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 810 : 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 969 : RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP &&
1453 50 : fabs((psOptions->dfMaxY - psOptions->dfMinY) / psOptions->dfYRes -
1454 50 : psOptions->nForceLines) <=
1455 810 : RELATIVE_ERROR_RES_SHARED_BY_COG_AND_GDALWARP;
1456 : }
1457 :
1458 : /************************************************************************/
1459 : /* CheckOptions() */
1460 : /************************************************************************/
1461 :
1462 892 : static bool CheckOptions(const char *pszDest, GDALDatasetH hDstDS,
1463 : int nSrcCount, GDALDatasetH *pahSrcDS,
1464 : GDALWarpAppOptions *psOptions, bool &bVRT,
1465 : int *pbUsageError)
1466 : {
1467 :
1468 892 : 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 1962 : if ((psOptions->osFormat.empty() &&
1479 1784 : EQUAL(CPLGetExtensionSafe(pszDest).c_str(), "VRT")) ||
1480 892 : (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 758 : if ((psOptions->nForcePixels != 0 || psOptions->nForceLines != 0) &&
1507 1675 : (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 892 : 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 891 : 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 891 : if (psOptions->dfErrorThreshold < 0)
1542 : {
1543 : // By default, use approximate transformer unless RPC_DEM is specified
1544 879 : if (psOptions->aosTransformerOptions.FetchNameValue("RPC_DEM") !=
1545 : nullptr)
1546 4 : psOptions->dfErrorThreshold = 0.0;
1547 : else
1548 875 : psOptions->dfErrorThreshold = 0.125;
1549 : }
1550 :
1551 : /* -------------------------------------------------------------------- */
1552 : /* -te_srs option */
1553 : /* -------------------------------------------------------------------- */
1554 891 : 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 891 : return true;
1639 : }
1640 :
1641 : /************************************************************************/
1642 : /* ProcessCutlineOptions() */
1643 : /************************************************************************/
1644 :
1645 891 : static bool ProcessCutlineOptions(int nSrcCount, GDALDatasetH *pahSrcDS,
1646 : GDALWarpAppOptions *psOptions,
1647 : OGRGeometryH &hCutline)
1648 : {
1649 891 : 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 885 : 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 880 : psOptions->aosWarpOptions.FetchNameValue("NUM_THREADS");
1677 880 : 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 880 : return true;
1686 : }
1687 :
1688 : /************************************************************************/
1689 : /* CreateOutput() */
1690 : /************************************************************************/
1691 :
1692 791 : static GDALDatasetH CreateOutput(const char *pszDest, int nSrcCount,
1693 : GDALDatasetH *pahSrcDS,
1694 : GDALWarpAppOptions *psOptions,
1695 : const bool bInitDestSetByUser,
1696 : void *&hUniqueTransformArg)
1697 : {
1698 791 : if (nSrcCount == 1 && !psOptions->bDisableSrcAlpha)
1699 : {
1700 1520 : if (GDALGetRasterCount(pahSrcDS[0]) > 0 &&
1701 760 : 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 791 : auto hDstDS = GDALWarpCreateOutput(
1713 : nSrcCount, pahSrcDS, pszDest, psOptions->osFormat.c_str(),
1714 : psOptions->aosTransformerOptions.List(),
1715 791 : psOptions->aosCreateOptions.List(), psOptions->eOutputType,
1716 791 : &hUniqueTransformArg, psOptions->bSetColorInterpretation, psOptions);
1717 791 : if (hDstDS == nullptr)
1718 : {
1719 10 : return nullptr;
1720 : }
1721 781 : psOptions->bCreateOutput = true;
1722 :
1723 781 : if (!bInitDestSetByUser)
1724 : {
1725 766 : if (psOptions->osDstNodata.empty())
1726 : {
1727 728 : psOptions->aosWarpOptions.SetNameValue("INIT_DEST", "0");
1728 : }
1729 : else
1730 : {
1731 38 : psOptions->aosWarpOptions.SetNameValue("INIT_DEST", "NO_DATA");
1732 : }
1733 : }
1734 :
1735 781 : return hDstDS;
1736 : }
1737 :
1738 : /************************************************************************/
1739 : /* ProcessMetadata() */
1740 : /************************************************************************/
1741 :
1742 889 : static void ProcessMetadata(int iSrc, GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
1743 : GDALWarpAppOptions *psOptions,
1744 : const bool bEnableDstAlpha)
1745 : {
1746 889 : if (psOptions->bCopyMetadata)
1747 : {
1748 889 : const char *pszSrcInfo = nullptr;
1749 889 : GDALRasterBandH hSrcBand = nullptr;
1750 889 : GDALRasterBandH hDstBand = nullptr;
1751 :
1752 : /* copy metadata from first dataset */
1753 889 : if (iSrc == 0)
1754 : {
1755 867 : CPLDebug(
1756 : "WARP",
1757 : "Copying metadata from first source to destination dataset");
1758 : /* copy dataset-level metadata */
1759 867 : char **papszMetadata = GDALGetMetadata(hSrcDS, nullptr);
1760 :
1761 867 : char **papszMetadataNew = nullptr;
1762 1879 : for (int i = 0;
1763 1879 : 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 867 : 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 867 : CSLDestroy(papszMetadataNew);
1791 :
1792 : /* ISIS3 -> ISIS3 special case */
1793 867 : 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 866 : 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 866 : 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 867 : if (GDALGetRasterCount(hSrcDS) == GDALGetRasterCount(hDstDS))
1814 : {
1815 1878 : for (int iBand = 0; iBand < GDALGetRasterCount(hSrcDS); iBand++)
1816 : {
1817 1065 : hSrcBand = GDALGetRasterBand(hSrcDS, iBand + 1);
1818 1065 : hDstBand = GDALGetRasterBand(hDstDS, iBand + 1);
1819 : /* copy metadata, except stats (#5319) */
1820 1065 : papszMetadata = GDALGetMetadata(hSrcBand, nullptr);
1821 1065 : 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 1065 : if (psOptions->bCopyBandInfo)
1838 : {
1839 1065 : pszSrcInfo = GDALGetDescription(hSrcBand);
1840 1065 : if (pszSrcInfo != nullptr && strlen(pszSrcInfo) > 0)
1841 2 : GDALSetDescription(hDstBand, pszSrcInfo);
1842 1065 : pszSrcInfo = GDALGetRasterUnitType(hSrcBand);
1843 1065 : 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 889 : }
1892 :
1893 : /************************************************************************/
1894 : /* SetupNoData() */
1895 : /************************************************************************/
1896 :
1897 886 : 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 913 : if (!psOptions->osSrcNodata.empty() &&
1904 27 : !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 886 : if (psOptions->osSrcNodata.empty())
1968 : {
1969 859 : int bHaveNodata = FALSE;
1970 859 : double dfReal = 0.0;
1971 :
1972 1918 : for (int i = 0; !bHaveNodata && i < psWO->nBandCount; i++)
1973 : {
1974 : GDALRasterBandH hBand =
1975 1059 : GDALGetRasterBand(hWrkSrcDS, psWO->panSrcBands[i]);
1976 1059 : dfReal = GDALGetRasterNoDataValue(hBand, &bHaveNodata);
1977 : }
1978 :
1979 859 : if (bHaveNodata)
1980 : {
1981 53 : 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 53 : psWO->padfSrcNoDataReal = static_cast<double *>(
1993 53 : CPLMalloc(psWO->nBandCount * sizeof(double)));
1994 :
1995 140 : for (int i = 0; i < psWO->nBandCount; i++)
1996 : {
1997 : GDALRasterBandH hBand =
1998 87 : GDALGetRasterBand(hWrkSrcDS, psWO->panSrcBands[i]);
1999 :
2000 87 : dfReal = GDALGetRasterNoDataValue(hBand, &bHaveNodata);
2001 :
2002 87 : if (bHaveNodata)
2003 : {
2004 87 : 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 933 : if (!psOptions->osDstNodata.empty() &&
2019 47 : !EQUAL(psOptions->osDstNodata.c_str(), "none"))
2020 : {
2021 47 : char **papszTokens = CSLTokenizeString(psOptions->osDstNodata.c_str());
2022 47 : const int nTokenCount = CSLCount(papszTokens);
2023 47 : bool bDstNoDataNone = true;
2024 :
2025 47 : psWO->padfDstNoDataReal =
2026 47 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
2027 47 : psWO->padfDstNoDataImag =
2028 47 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
2029 :
2030 94 : for (int i = 0; i < psWO->nBandCount; i++)
2031 : {
2032 47 : psWO->padfDstNoDataReal[i] = -1.1e20;
2033 47 : psWO->padfDstNoDataImag[i] = 0.0;
2034 :
2035 47 : if (i < nTokenCount)
2036 : {
2037 47 : 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 47 : 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 47 : CPLStringToComplex(papszTokens[i], psWO->padfDstNoDataReal + i,
2052 47 : psWO->padfDstNoDataImag + i);
2053 94 : psWO->padfDstNoDataReal[i] =
2054 47 : GDALAdjustNoDataCloseToFloatMax(psWO->padfDstNoDataReal[i]);
2055 94 : psWO->padfDstNoDataImag[i] =
2056 47 : GDALAdjustNoDataCloseToFloatMax(psWO->padfDstNoDataImag[i]);
2057 47 : bDstNoDataNone = false;
2058 47 : CPLDebug("WARP", "dstnodata of band %d set to %f", i,
2059 47 : 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 47 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]);
2079 47 : int bClamped = FALSE;
2080 47 : int bRounded = FALSE;
2081 47 : psWO->padfDstNoDataReal[i] = GDALAdjustValueToDataType(
2082 47 : GDALGetRasterDataType(hBand), psWO->padfDstNoDataReal[i],
2083 : &bClamped, &bRounded);
2084 :
2085 47 : 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 47 : 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 47 : if (psOptions->bCreateOutput && iSrc == 0)
2104 : {
2105 43 : GDALSetRasterNoDataValue(
2106 43 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]),
2107 43 : psWO->padfDstNoDataReal[i]);
2108 : }
2109 : }
2110 :
2111 47 : CSLDestroy(papszTokens);
2112 : }
2113 :
2114 : /* check if the output dataset has already nodata */
2115 886 : 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 1725 : if (psOptions->osDstNodata.empty() && psWO->padfSrcNoDataReal != nullptr &&
2146 64 : psWO->padfDstNoDataReal != nullptr && psOptions->bCreateOutput &&
2147 1725 : 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 884 : else if (psOptions->osDstNodata.empty() &&
2190 837 : psWO->padfSrcNoDataReal != nullptr &&
2191 1721 : psWO->padfDstNoDataReal == nullptr && !bEnableDstAlpha)
2192 : {
2193 48 : psWO->padfDstNoDataReal =
2194 48 : static_cast<double *>(CPLMalloc(psWO->nBandCount * sizeof(double)));
2195 :
2196 48 : if (psWO->padfSrcNoDataImag != nullptr)
2197 : {
2198 1 : psWO->padfDstNoDataImag = static_cast<double *>(
2199 1 : CPLMalloc(psWO->nBandCount * sizeof(double)));
2200 : }
2201 :
2202 48 : if (!psOptions->bQuiet)
2203 15 : printf("Copying nodata values from source %s to destination %s.\n",
2204 : GDALGetDescription(hSrcDS), pszDest);
2205 :
2206 115 : for (int i = 0; i < psWO->nBandCount; i++)
2207 : {
2208 67 : psWO->padfDstNoDataReal[i] = psWO->padfSrcNoDataReal[i];
2209 67 : if (psWO->padfSrcNoDataImag != nullptr)
2210 : {
2211 1 : psWO->padfDstNoDataImag[i] = psWO->padfSrcNoDataImag[i];
2212 : }
2213 67 : CPLDebug("WARP", "srcNoData=%f dstNoData=%f",
2214 67 : psWO->padfSrcNoDataReal[i], psWO->padfDstNoDataReal[i]);
2215 :
2216 67 : if (psOptions->bCreateOutput && iSrc == 0)
2217 : {
2218 67 : CPLDebug("WARP",
2219 : "calling GDALSetRasterNoDataValue() for band#%d", i);
2220 67 : GDALSetRasterNoDataValue(
2221 67 : GDALGetRasterBand(hDstDS, psWO->panDstBands[i]),
2222 67 : psWO->padfDstNoDataReal[i]);
2223 : }
2224 : }
2225 :
2226 48 : 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 45 : psWO->papszWarpOptions =
2232 45 : CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
2233 : }
2234 : }
2235 886 : }
2236 :
2237 : /************************************************************************/
2238 : /* SetupSkipNoSource() */
2239 : /************************************************************************/
2240 :
2241 886 : static void SetupSkipNoSource(int iSrc, GDALDatasetH hDstDS,
2242 : GDALWarpOptions *psWO,
2243 : GDALWarpAppOptions *psOptions)
2244 : {
2245 801 : if (psOptions->bCreateOutput && iSrc == 0 &&
2246 779 : CSLFetchNameValue(psWO->papszWarpOptions, "SKIP_NOSOURCE") == nullptr &&
2247 775 : CSLFetchNameValue(psWO->papszWarpOptions, "STREAMABLE_OUTPUT") ==
2248 1687 : nullptr &&
2249 : // This white list of drivers could potentially be extended.
2250 774 : (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 679 : bool bOKRegardingInitDest = false;
2258 : const char *pszInitDest =
2259 679 : CSLFetchNameValue(psWO->papszWarpOptions, "INIT_DEST");
2260 679 : if (pszInitDest == nullptr || EQUAL(pszInitDest, "NO_DATA"))
2261 : {
2262 71 : 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 71 : if (EQUAL(psOptions->osFormat.c_str(), "MEM"))
2267 : {
2268 83 : for (int i = 0; i < GDALGetRasterCount(hDstDS); i++)
2269 : {
2270 57 : int bHasNoData = false;
2271 57 : double dfDstNoDataVal = GDALGetRasterNoDataValue(
2272 : GDALGetRasterBand(hDstDS, i + 1), &bHasNoData);
2273 57 : if (bHasNoData && dfDstNoDataVal != 0.0)
2274 : {
2275 27 : bOKRegardingInitDest = false;
2276 27 : break;
2277 : }
2278 : }
2279 71 : }
2280 : }
2281 : else
2282 : {
2283 608 : char **papszTokensInitDest = CSLTokenizeString(pszInitDest);
2284 608 : const int nTokenCountInitDest = CSLCount(papszTokensInitDest);
2285 608 : if (nTokenCountInitDest == 1 ||
2286 0 : nTokenCountInitDest == GDALGetRasterCount(hDstDS))
2287 : {
2288 608 : bOKRegardingInitDest = true;
2289 1424 : for (int i = 0; i < GDALGetRasterCount(hDstDS); i++)
2290 : {
2291 822 : double dfInitVal = GDALAdjustNoDataCloseToFloatMax(
2292 822 : CPLAtofM(papszTokensInitDest[std::min(
2293 822 : i, nTokenCountInitDest - 1)]));
2294 822 : int bHasNoData = false;
2295 822 : double dfDstNoDataVal = GDALGetRasterNoDataValue(
2296 : GDALGetRasterBand(hDstDS, i + 1), &bHasNoData);
2297 822 : if (!((bHasNoData && dfInitVal == dfDstNoDataVal) ||
2298 820 : (!bHasNoData && dfInitVal == 0.0)))
2299 : {
2300 5 : bOKRegardingInitDest = false;
2301 6 : break;
2302 : }
2303 817 : if (EQUAL(psOptions->osFormat.c_str(), "MEM") &&
2304 817 : bHasNoData && dfDstNoDataVal != 0.0)
2305 : {
2306 1 : bOKRegardingInitDest = false;
2307 1 : break;
2308 : }
2309 : }
2310 : }
2311 608 : CSLDestroy(papszTokensInitDest);
2312 : }
2313 :
2314 679 : if (bOKRegardingInitDest)
2315 : {
2316 646 : CPLDebug("GDALWARP", "Defining SKIP_NOSOURCE=YES");
2317 646 : psWO->papszWarpOptions =
2318 646 : CSLSetNameValue(psWO->papszWarpOptions, "SKIP_NOSOURCE", "YES");
2319 : }
2320 : }
2321 886 : }
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 886 : 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 886 : if (CPLTestBool(CSLFetchNameValueDef(psWO->papszWarpOptions,
2338 756 : "SKIP_NOSOURCE", "NO")) &&
2339 756 : GDALGetMetadata(hSrcDS, "RPC") != nullptr &&
2340 12 : EQUAL(
2341 : psOptions->aosTransformerOptions.FetchNameValueDef("METHOD", "RPC"),
2342 1642 : "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 885 : return true;
2397 : }
2398 :
2399 : /************************************************************************/
2400 : /* GDALWarpDirect() */
2401 : /************************************************************************/
2402 :
2403 892 : static GDALDatasetH GDALWarpDirect(const char *pszDest, GDALDatasetH hDstDS,
2404 : int nSrcCount, GDALDatasetH *pahSrcDS,
2405 : GDALWarpAppOptions *psOptions,
2406 : int *pbUsageError)
2407 : {
2408 892 : CPLErrorReset();
2409 892 : 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 892 : if (pszDest == nullptr)
2419 86 : pszDest = GDALGetDescription(hDstDS);
2420 :
2421 : #ifdef DEBUG
2422 892 : GDALDataset *poDstDS = GDALDataset::FromHandle(hDstDS);
2423 : const int nExpectedRefCountAtEnd =
2424 892 : (poDstDS != nullptr) ? poDstDS->GetRefCount() : 1;
2425 : (void)nExpectedRefCountAtEnd;
2426 : #endif
2427 892 : const bool bDropDstDSRef = (hDstDS != nullptr);
2428 892 : if (hDstDS != nullptr)
2429 89 : GDALReferenceDataset(hDstDS);
2430 :
2431 : #if defined(USE_PROJ_BASED_VERTICAL_SHIFT_METHOD)
2432 892 : if (psOptions->bNoVShift)
2433 : {
2434 0 : psOptions->aosTransformerOptions.SetNameValue("STRIP_VERT_CS", "YES");
2435 : }
2436 892 : else if (nSrcCount)
2437 : {
2438 886 : bool bSrcHasVertAxis = false;
2439 886 : bool bDstHasVertAxis = false;
2440 1772 : OGRSpatialReference oSRSSrc;
2441 1772 : OGRSpatialReference oSRSDst;
2442 :
2443 886 : 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 892 : bool bVRT = false;
2455 892 : 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 891 : OGRGeometryH hCutline = nullptr;
2466 891 : 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 880 : void *hUniqueTransformArg = nullptr;
2476 : const bool bInitDestSetByUser =
2477 880 : (psOptions->aosWarpOptions.FetchNameValue("INIT_DEST") != nullptr);
2478 :
2479 880 : const bool bFigureoutCorrespondingWindow =
2480 1671 : (hDstDS != nullptr) ||
2481 791 : (((psOptions->nForcePixels != 0 && psOptions->nForceLines != 0) ||
2482 659 : (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 880 : psOptions->aosTransformerOptions.FetchNameValue("METHOD");
2488 17 : if (pszMethod && EQUAL(pszMethod, "GCP_TPS") &&
2489 902 : 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 880 : if (hDstDS == nullptr)
2499 : {
2500 791 : hDstDS = CreateOutput(pszDest, nSrcCount, pahSrcDS, psOptions,
2501 : bInitDestSetByUser, hUniqueTransformArg);
2502 791 : 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 781 : poDstDS = GDALDataset::FromHandle(hDstDS);
2511 781 : 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 870 : bool bEnableDstAlpha = psOptions->bEnableDstAlpha;
2528 812 : if (!bEnableDstAlpha && GDALGetRasterCount(hDstDS) &&
2529 811 : GDALGetRasterColorInterpretation(GDALGetRasterBand(
2530 1682 : 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 19662 : int Do(double dfComplete)
2552 : {
2553 39324 : CPLString osMsg;
2554 : osMsg.Printf("Processing %s [%d/%d]",
2555 19662 : CPLGetFilename(GDALGetDescription(pahSrcDS[iSrc])),
2556 19662 : iSrc + 1, nSrcCount);
2557 19662 : return pfnExternalProgress((iSrc + dfComplete) / nSrcCount,
2558 39324 : 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 870 : oProgress.pfnExternalProgress = psOptions->pfnProgress;
2570 870 : oProgress.pExternalProgressData = psOptions->pProgressData;
2571 870 : oProgress.nSrcCount = nSrcCount;
2572 870 : oProgress.pahSrcDS = pahSrcDS;
2573 :
2574 : /* -------------------------------------------------------------------- */
2575 : /* Loop over all source files, processing each in turn. */
2576 : /* -------------------------------------------------------------------- */
2577 870 : GDALTransformerFunc pfnTransformer = nullptr;
2578 870 : void *hTransformArg = nullptr;
2579 870 : bool bHasGotErr = false;
2580 1667 : for (int iSrc = 0; iSrc < nSrcCount; iSrc++)
2581 : {
2582 : GDALDatasetH hSrcDS;
2583 :
2584 : /* --------------------------------------------------------------------
2585 : */
2586 : /* Open this file. */
2587 : /* --------------------------------------------------------------------
2588 : */
2589 890 : hSrcDS = pahSrcDS[iSrc];
2590 890 : oProgress.iSrc = iSrc;
2591 :
2592 : /* --------------------------------------------------------------------
2593 : */
2594 : /* Check that there's at least one raster band */
2595 : /* --------------------------------------------------------------------
2596 : */
2597 890 : 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 889 : bool bEnableSrcAlpha = psOptions->bEnableSrcAlpha;
2613 889 : if (GDALGetRasterColorInterpretation(GDALGetRasterBand(
2614 59 : hSrcDS, GDALGetRasterCount(hSrcDS))) == GCI_AlphaBand &&
2615 889 : !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 889 : 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 2215 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
2643 1307 : psOptions->eResampleAlg != GRA_Mode &&
2644 418 : 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 892 : 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 889 : if (hUniqueTransformArg)
2693 759 : hTransformArg = hUniqueTransformArg;
2694 : else
2695 130 : hTransformArg = GDALCreateGenImgProjTransformer2(
2696 : hSrcDS, hDstDS, psOptions->aosTransformerOptions.List());
2697 :
2698 889 : if (hTransformArg == nullptr)
2699 : {
2700 0 : OGR_G_DestroyGeometry(hCutline);
2701 0 : GDALReleaseDataset(hDstDS);
2702 0 : return nullptr;
2703 : }
2704 :
2705 889 : pfnTransformer = GDALGenImgProjTransform;
2706 :
2707 : // Check if transformation is inversible
2708 : {
2709 889 : double dfX = GDALGetRasterXSize(hDstDS) / 2.0;
2710 889 : double dfY = GDALGetRasterYSize(hDstDS) / 2.0;
2711 889 : double dfZ = 0;
2712 889 : int bSuccess = false;
2713 889 : const auto nErrorCounterBefore = CPLGetErrorCounter();
2714 889 : pfnTransformer(hTransformArg, TRUE, 1, &dfX, &dfY, &dfZ, &bSuccess);
2715 889 : 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 889 : GDALDataset *poSrcDS = static_cast<GDALDataset *>(hSrcDS);
2732 889 : GDALDataset *poSrcOvrDS = nullptr;
2733 889 : int nOvCount = poSrcDS->GetRasterBand(1)->GetOverviewCount();
2734 889 : 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 4 : std::vector<int> abSuccess(nPoints);
2765 4 : pfnTransformer(hTransformArg, TRUE, nPoints, &adfX[0], &adfY[0],
2766 4 : &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 4 : dfTargetRatioX =
2785 4 : (dfMaxSrcX - dfMinSrcX) / GDALGetRasterXSize(hDstDS);
2786 : }
2787 4 : if (dfMaxSrcY > dfMinSrcY)
2788 : {
2789 4 : dfTargetRatioY =
2790 4 : (dfMaxSrcY - dfMinSrcY) / GDALGetRasterYSize(hDstDS);
2791 : }
2792 : // take the minimum of these ratios #7019
2793 4 : dfTargetRatio = std::min(dfTargetRatioX, dfTargetRatioY);
2794 : }
2795 : else
2796 : {
2797 : /* Compute what the "natural" output resolution (in pixels)
2798 : * would be for this */
2799 : /* input dataset */
2800 : double adfSuggestedGeoTransform[6];
2801 : int nPixels, nLines;
2802 17 : if (GDALSuggestedWarpOutput(
2803 : hSrcDS, pfnTransformer, hTransformArg,
2804 17 : adfSuggestedGeoTransform, &nPixels, &nLines) == CE_None)
2805 : {
2806 17 : dfTargetRatio = 1.0 / adfSuggestedGeoTransform[1];
2807 : }
2808 : }
2809 :
2810 21 : if (dfTargetRatio > 1.0)
2811 : {
2812 : // Note: keep this logic for overview selection in sync between
2813 : // gdalwarp_lib.cpp and rasterio.cpp
2814 15 : const char *pszOversampligThreshold = CPLGetConfigOption(
2815 : "GDALWARP_OVERSAMPLING_THRESHOLD", nullptr);
2816 : const double dfOversamplingThreshold =
2817 15 : pszOversampligThreshold ? CPLAtof(pszOversampligThreshold)
2818 15 : : 1.0;
2819 :
2820 15 : int iBestOvr = -1;
2821 15 : double dfBestRatio = 0;
2822 35 : for (int iOvr = -1; iOvr < nOvCount; iOvr++)
2823 : {
2824 : const double dfOvrRatio =
2825 : iOvr < 0
2826 31 : ? 1.0
2827 16 : : static_cast<double>(poSrcDS->GetRasterXSize()) /
2828 16 : poSrcDS->GetRasterBand(1)
2829 16 : ->GetOverview(iOvr)
2830 16 : ->GetXSize();
2831 :
2832 : // Is it nearly the requested factor and better (lower) than
2833 : // the current best factor?
2834 : // Use an epsilon because of numerical instability.
2835 31 : constexpr double EPSILON = 1e-1;
2836 35 : if (dfOvrRatio >=
2837 31 : dfTargetRatio * dfOversamplingThreshold + EPSILON ||
2838 : dfOvrRatio <= dfBestRatio)
2839 : {
2840 4 : continue;
2841 : }
2842 :
2843 27 : iBestOvr = iOvr;
2844 27 : dfBestRatio = dfOvrRatio;
2845 27 : if (std::abs(dfTargetRatio - dfOvrRatio) < EPSILON)
2846 : {
2847 11 : break;
2848 : }
2849 : }
2850 15 : const int iOvr =
2851 15 : iBestOvr + (psOptions->nOvLevel - OVR_LEVEL_AUTO);
2852 15 : if (iOvr >= 0)
2853 : {
2854 9 : CPLDebug("WARP", "Selecting overview level %d for %s", iOvr,
2855 : GDALGetDescription(hSrcDS));
2856 : poSrcOvrDS =
2857 9 : GDALCreateOverviewDataset(poSrcDS, iOvr,
2858 : /* bThisLevelOnly = */ false);
2859 : }
2860 21 : }
2861 : }
2862 868 : else if (psOptions->nOvLevel >= 0)
2863 : {
2864 6 : poSrcOvrDS = GDALCreateOverviewDataset(poSrcDS, psOptions->nOvLevel,
2865 : /* bThisLevelOnly = */ true);
2866 6 : if (poSrcOvrDS == nullptr)
2867 : {
2868 1 : if (!psOptions->bQuiet)
2869 : {
2870 1 : CPLError(CE_Warning, CPLE_AppDefined,
2871 : "cannot get overview level %d for "
2872 : "dataset %s. Defaulting to level %d",
2873 : psOptions->nOvLevel, GDALGetDescription(hSrcDS),
2874 : nOvCount - 1);
2875 : }
2876 1 : if (nOvCount > 0)
2877 : poSrcOvrDS =
2878 1 : GDALCreateOverviewDataset(poSrcDS, nOvCount - 1,
2879 : /* bThisLevelOnly = */ false);
2880 : }
2881 : else
2882 : {
2883 5 : CPLDebug("WARP", "Selecting overview level %d for %s",
2884 : psOptions->nOvLevel, GDALGetDescription(hSrcDS));
2885 : }
2886 : }
2887 :
2888 889 : if (poSrcOvrDS == nullptr)
2889 874 : GDALReferenceDataset(hSrcDS);
2890 :
2891 889 : GDALDatasetH hWrkSrcDS =
2892 889 : poSrcOvrDS ? static_cast<GDALDatasetH>(poSrcOvrDS) : hSrcDS;
2893 :
2894 : #if !defined(USE_PROJ_BASED_VERTICAL_SHIFT_METHOD)
2895 : if (!psOptions->bNoVShift)
2896 : {
2897 : bool bErrorOccurred = false;
2898 : hWrkSrcDS = ApplyVerticalShiftGrid(
2899 : hWrkSrcDS, psOptions, bVRT ? hDstDS : nullptr, bErrorOccurred);
2900 : if (bErrorOccurred)
2901 : {
2902 : GDALDestroyTransformer(hTransformArg);
2903 : OGR_G_DestroyGeometry(hCutline);
2904 : GDALReleaseDataset(hWrkSrcDS);
2905 : GDALReleaseDataset(hDstDS);
2906 : return nullptr;
2907 : }
2908 : }
2909 : #endif
2910 :
2911 : /* --------------------------------------------------------------------
2912 : */
2913 : /* Clear temporary INIT_DEST settings after the first image. */
2914 : /* --------------------------------------------------------------------
2915 : */
2916 889 : if (psOptions->bCreateOutput && iSrc == 1)
2917 22 : psOptions->aosWarpOptions.SetNameValue("INIT_DEST", nullptr);
2918 :
2919 : /* --------------------------------------------------------------------
2920 : */
2921 : /* Define SKIP_NOSOURCE after the first image (since
2922 : * initialization*/
2923 : /* has already be done). */
2924 : /* --------------------------------------------------------------------
2925 : */
2926 889 : if (iSrc == 1 && psOptions->aosWarpOptions.FetchNameValue(
2927 : "SKIP_NOSOURCE") == nullptr)
2928 : {
2929 22 : CPLDebug("GDALWARP", "Defining SKIP_NOSOURCE=YES");
2930 22 : psOptions->aosWarpOptions.SetNameValue("SKIP_NOSOURCE", "YES");
2931 : }
2932 :
2933 : /* --------------------------------------------------------------------
2934 : */
2935 : /* Setup warp options. */
2936 : /* --------------------------------------------------------------------
2937 : */
2938 889 : GDALWarpOptions *psWO = GDALCreateWarpOptions();
2939 :
2940 889 : psWO->papszWarpOptions = CSLDuplicate(psOptions->aosWarpOptions.List());
2941 889 : psWO->eWorkingDataType = psOptions->eWorkingType;
2942 :
2943 889 : psWO->eResampleAlg = psOptions->eResampleAlg;
2944 :
2945 889 : psWO->hSrcDS = hWrkSrcDS;
2946 889 : psWO->hDstDS = hDstDS;
2947 :
2948 889 : if (!bVRT)
2949 : {
2950 801 : if (psOptions->pfnProgress == GDALDummyProgress)
2951 : {
2952 677 : psWO->pfnProgress = GDALDummyProgress;
2953 677 : psWO->pProgressArg = nullptr;
2954 : }
2955 : else
2956 : {
2957 124 : psWO->pfnProgress = Progress::ProgressFunc;
2958 124 : psWO->pProgressArg = &oProgress;
2959 : }
2960 : }
2961 :
2962 889 : if (psOptions->dfWarpMemoryLimit != 0.0)
2963 31 : psWO->dfWarpMemoryLimit = psOptions->dfWarpMemoryLimit;
2964 :
2965 : /* --------------------------------------------------------------------
2966 : */
2967 : /* Setup band mapping. */
2968 : /* --------------------------------------------------------------------
2969 : */
2970 889 : if (psOptions->anSrcBands.empty())
2971 : {
2972 871 : if (bEnableSrcAlpha)
2973 53 : psWO->nBandCount = GDALGetRasterCount(hWrkSrcDS) - 1;
2974 : else
2975 818 : psWO->nBandCount = GDALGetRasterCount(hWrkSrcDS);
2976 : }
2977 : else
2978 : {
2979 18 : psWO->nBandCount = static_cast<int>(psOptions->anSrcBands.size());
2980 : }
2981 :
2982 889 : const int nNeededDstBands =
2983 889 : psWO->nBandCount + (bEnableDstAlpha ? 1 : 0);
2984 889 : if (nNeededDstBands > GDALGetRasterCount(hDstDS))
2985 : {
2986 1 : CPLError(CE_Failure, CPLE_AppDefined,
2987 : "Destination dataset has %d bands, but at least %d "
2988 : "are needed",
2989 : GDALGetRasterCount(hDstDS), nNeededDstBands);
2990 1 : GDALDestroyTransformer(hTransformArg);
2991 1 : GDALDestroyWarpOptions(psWO);
2992 1 : OGR_G_DestroyGeometry(hCutline);
2993 1 : GDALReleaseDataset(hWrkSrcDS);
2994 1 : GDALReleaseDataset(hDstDS);
2995 1 : return nullptr;
2996 : }
2997 :
2998 888 : psWO->panSrcBands =
2999 888 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
3000 888 : psWO->panDstBands =
3001 888 : static_cast<int *>(CPLMalloc(psWO->nBandCount * sizeof(int)));
3002 888 : if (psOptions->anSrcBands.empty())
3003 : {
3004 1981 : for (int i = 0; i < psWO->nBandCount; i++)
3005 : {
3006 1111 : psWO->panSrcBands[i] = i + 1;
3007 1111 : psWO->panDstBands[i] = i + 1;
3008 : }
3009 : }
3010 : else
3011 : {
3012 45 : for (int i = 0; i < psWO->nBandCount; i++)
3013 : {
3014 58 : if (psOptions->anSrcBands[i] <= 0 ||
3015 29 : psOptions->anSrcBands[i] > GDALGetRasterCount(hSrcDS))
3016 : {
3017 1 : CPLError(CE_Failure, CPLE_AppDefined,
3018 : "-srcband[%d] = %d is invalid", i,
3019 1 : psOptions->anSrcBands[i]);
3020 1 : GDALDestroyTransformer(hTransformArg);
3021 1 : GDALDestroyWarpOptions(psWO);
3022 1 : OGR_G_DestroyGeometry(hCutline);
3023 1 : GDALReleaseDataset(hWrkSrcDS);
3024 1 : GDALReleaseDataset(hDstDS);
3025 1 : return nullptr;
3026 : }
3027 56 : if (psOptions->anDstBands[i] <= 0 ||
3028 28 : psOptions->anDstBands[i] > GDALGetRasterCount(hDstDS))
3029 : {
3030 1 : CPLError(CE_Failure, CPLE_AppDefined,
3031 : "-dstband[%d] = %d is invalid", i,
3032 1 : psOptions->anDstBands[i]);
3033 1 : GDALDestroyTransformer(hTransformArg);
3034 1 : GDALDestroyWarpOptions(psWO);
3035 1 : OGR_G_DestroyGeometry(hCutline);
3036 1 : GDALReleaseDataset(hWrkSrcDS);
3037 1 : GDALReleaseDataset(hDstDS);
3038 1 : return nullptr;
3039 : }
3040 27 : psWO->panSrcBands[i] = psOptions->anSrcBands[i];
3041 27 : psWO->panDstBands[i] = psOptions->anDstBands[i];
3042 : }
3043 : }
3044 :
3045 : /* --------------------------------------------------------------------
3046 : */
3047 : /* Setup alpha bands used if any. */
3048 : /* --------------------------------------------------------------------
3049 : */
3050 886 : if (bEnableSrcAlpha)
3051 56 : psWO->nSrcAlphaBand = GDALGetRasterCount(hWrkSrcDS);
3052 :
3053 886 : if (bEnableDstAlpha)
3054 : {
3055 102 : if (psOptions->anSrcBands.empty())
3056 98 : psWO->nDstAlphaBand = GDALGetRasterCount(hDstDS);
3057 : else
3058 4 : psWO->nDstAlphaBand =
3059 4 : static_cast<int>(psOptions->anDstBands.size()) + 1;
3060 : }
3061 :
3062 : /* --------------------------------------------------------------------
3063 : */
3064 : /* Setup NODATA options. */
3065 : /* --------------------------------------------------------------------
3066 : */
3067 886 : SetupNoData(pszDest, iSrc, hSrcDS, hWrkSrcDS, hDstDS, psWO, psOptions,
3068 : bEnableDstAlpha, bInitDestSetByUser);
3069 :
3070 886 : oProgress.Do(0);
3071 :
3072 : /* --------------------------------------------------------------------
3073 : */
3074 : /* For the first source image of a newly created dataset, decide */
3075 : /* if we can safely enable SKIP_NOSOURCE optimization. */
3076 : /* --------------------------------------------------------------------
3077 : */
3078 886 : SetupSkipNoSource(iSrc, hDstDS, psWO, psOptions);
3079 :
3080 : /* --------------------------------------------------------------------
3081 : */
3082 : /* In some cases, RPC evaluation can find valid input pixel for */
3083 : /* output pixels that are outside the footprint of the source */
3084 : /* dataset, so limit the area we update in the target dataset from
3085 : */
3086 : /* the suggested warp output (only in cases where
3087 : * SKIP_NOSOURCE=YES) */
3088 : /* --------------------------------------------------------------------
3089 : */
3090 886 : int nWarpDstXOff = 0;
3091 886 : int nWarpDstYOff = 0;
3092 886 : int nWarpDstXSize = GDALGetRasterXSize(hDstDS);
3093 886 : int nWarpDstYSize = GDALGetRasterYSize(hDstDS);
3094 :
3095 886 : if (!AdjustOutputExtentForRPC(
3096 : hSrcDS, hDstDS, pfnTransformer, hTransformArg, psWO, psOptions,
3097 : nWarpDstXOff, nWarpDstYOff, nWarpDstXSize, nWarpDstYSize))
3098 : {
3099 1 : GDALDestroyTransformer(hTransformArg);
3100 1 : GDALDestroyWarpOptions(psWO);
3101 1 : GDALReleaseDataset(hWrkSrcDS);
3102 1 : continue;
3103 : }
3104 :
3105 : /* We need to recreate the transform when operating on an overview */
3106 885 : if (poSrcOvrDS != nullptr)
3107 : {
3108 15 : GDALDestroyGenImgProjTransformer(hTransformArg);
3109 15 : hTransformArg = GDALCreateGenImgProjTransformer2(
3110 : hWrkSrcDS, hDstDS, psOptions->aosTransformerOptions.List());
3111 : }
3112 :
3113 885 : bool bUseApproxTransformer = psOptions->dfErrorThreshold != 0.0;
3114 : #ifdef USE_PROJ_BASED_VERTICAL_SHIFT_METHOD
3115 885 : if (!psOptions->bNoVShift)
3116 : {
3117 : // Can modify psWO->papszWarpOptions
3118 885 : if (ApplyVerticalShift(hWrkSrcDS, psOptions, psWO))
3119 : {
3120 19 : bUseApproxTransformer = false;
3121 : }
3122 : }
3123 : #endif
3124 :
3125 : /* --------------------------------------------------------------------
3126 : */
3127 : /* Warp the transformer with a linear approximator unless the */
3128 : /* acceptable error is zero. */
3129 : /* --------------------------------------------------------------------
3130 : */
3131 885 : if (bUseApproxTransformer)
3132 : {
3133 853 : hTransformArg = GDALCreateApproxTransformer(
3134 : GDALGenImgProjTransform, hTransformArg,
3135 : psOptions->dfErrorThreshold);
3136 853 : pfnTransformer = GDALApproxTransform;
3137 853 : GDALApproxTransformerOwnsSubtransformer(hTransformArg, TRUE);
3138 : }
3139 :
3140 885 : psWO->pfnTransformer = pfnTransformer;
3141 885 : psWO->pTransformerArg = hTransformArg;
3142 :
3143 : /* --------------------------------------------------------------------
3144 : */
3145 : /* If we have a cutline, transform it into the source */
3146 : /* pixel/line coordinate system and insert into warp options. */
3147 : /* --------------------------------------------------------------------
3148 : */
3149 885 : if (hCutline != nullptr)
3150 : {
3151 : CPLErr eError;
3152 33 : eError = TransformCutlineToSource(
3153 : GDALDataset::FromHandle(hWrkSrcDS),
3154 : OGRGeometry::FromHandle(hCutline), &(psWO->papszWarpOptions),
3155 33 : psOptions->aosTransformerOptions.List());
3156 33 : if (eError == CE_Failure)
3157 : {
3158 1 : GDALDestroyTransformer(hTransformArg);
3159 1 : GDALDestroyWarpOptions(psWO);
3160 1 : OGR_G_DestroyGeometry(hCutline);
3161 1 : GDALReleaseDataset(hWrkSrcDS);
3162 1 : GDALReleaseDataset(hDstDS);
3163 1 : return nullptr;
3164 : }
3165 : }
3166 :
3167 : /* --------------------------------------------------------------------
3168 : */
3169 : /* If we are producing VRT output, then just initialize it with */
3170 : /* the warp options and write out now rather than proceeding */
3171 : /* with the operations. */
3172 : /* --------------------------------------------------------------------
3173 : */
3174 884 : if (bVRT)
3175 : {
3176 88 : GDALSetMetadataItem(hDstDS, "SrcOvrLevel",
3177 : CPLSPrintf("%d", psOptions->nOvLevel), nullptr);
3178 88 : CPLErr eErr = GDALInitializeWarpedVRT(hDstDS, psWO);
3179 88 : GDALDestroyWarpOptions(psWO);
3180 88 : OGR_G_DestroyGeometry(hCutline);
3181 88 : GDALReleaseDataset(hWrkSrcDS);
3182 88 : if (eErr != CE_None)
3183 : {
3184 1 : GDALDestroyTransformer(hTransformArg);
3185 1 : GDALReleaseDataset(hDstDS);
3186 1 : return nullptr;
3187 : }
3188 : // In case of success, hDstDS has become the owner of hTransformArg
3189 : // so do not free it.
3190 87 : if (!EQUAL(pszDest, ""))
3191 : {
3192 : const bool bWasFailureBefore =
3193 20 : (CPLGetLastErrorType() == CE_Failure);
3194 20 : GDALFlushCache(hDstDS);
3195 20 : if (!bWasFailureBefore && CPLGetLastErrorType() == CE_Failure)
3196 : {
3197 1 : GDALReleaseDataset(hDstDS);
3198 1 : hDstDS = nullptr;
3199 : }
3200 : }
3201 :
3202 87 : if (hDstDS)
3203 86 : oProgress.Do(1);
3204 :
3205 87 : return hDstDS;
3206 : }
3207 :
3208 : /* --------------------------------------------------------------------
3209 : */
3210 : /* Initialize and execute the warp. */
3211 : /* --------------------------------------------------------------------
3212 : */
3213 1592 : GDALWarpOperation oWO;
3214 :
3215 796 : if (oWO.Initialize(psWO) == CE_None)
3216 : {
3217 : CPLErr eErr;
3218 793 : if (psOptions->bMulti)
3219 5 : eErr = oWO.ChunkAndWarpMulti(nWarpDstXOff, nWarpDstYOff,
3220 : nWarpDstXSize, nWarpDstYSize);
3221 : else
3222 788 : eErr = oWO.ChunkAndWarpImage(nWarpDstXOff, nWarpDstYOff,
3223 : nWarpDstXSize, nWarpDstYSize);
3224 793 : if (eErr != CE_None)
3225 2 : bHasGotErr = true;
3226 : }
3227 : else
3228 : {
3229 3 : bHasGotErr = true;
3230 : }
3231 :
3232 : /* --------------------------------------------------------------------
3233 : */
3234 : /* Cleanup */
3235 : /* --------------------------------------------------------------------
3236 : */
3237 796 : GDALDestroyTransformer(hTransformArg);
3238 :
3239 796 : GDALDestroyWarpOptions(psWO);
3240 :
3241 796 : GDALReleaseDataset(hWrkSrcDS);
3242 : }
3243 :
3244 : /* -------------------------------------------------------------------- */
3245 : /* Final Cleanup. */
3246 : /* -------------------------------------------------------------------- */
3247 777 : const bool bWasFailureBefore = (CPLGetLastErrorType() == CE_Failure);
3248 777 : GDALFlushCache(hDstDS);
3249 777 : if (!bWasFailureBefore && CPLGetLastErrorType() == CE_Failure)
3250 : {
3251 1 : bHasGotErr = true;
3252 : }
3253 :
3254 777 : OGR_G_DestroyGeometry(hCutline);
3255 :
3256 777 : if (bHasGotErr || bDropDstDSRef)
3257 93 : GDALReleaseDataset(hDstDS);
3258 :
3259 : #ifdef DEBUG
3260 777 : if (!bHasGotErr || bDropDstDSRef)
3261 : {
3262 771 : CPLAssert(poDstDS->GetRefCount() == nExpectedRefCountAtEnd);
3263 : }
3264 : #endif
3265 :
3266 777 : return bHasGotErr ? nullptr : hDstDS;
3267 : }
3268 :
3269 : /************************************************************************/
3270 : /* ValidateCutline() */
3271 : /* Same as OGR_G_IsValid() except that it processes polygon per polygon*/
3272 : /* without paying attention to MultiPolygon specific validity rules. */
3273 : /************************************************************************/
3274 :
3275 196 : static bool ValidateCutline(const OGRGeometry *poGeom, bool bVerbose)
3276 : {
3277 196 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3278 196 : if (eType == wkbMultiPolygon)
3279 : {
3280 151 : for (const auto *poSubGeom : *(poGeom->toMultiPolygon()))
3281 : {
3282 81 : if (!ValidateCutline(poSubGeom, bVerbose))
3283 5 : return false;
3284 : }
3285 : }
3286 121 : else if (eType == wkbPolygon)
3287 : {
3288 120 : if (OGRGeometryFactory::haveGEOS() && !poGeom->IsValid())
3289 : {
3290 6 : if (!bVerbose)
3291 6 : return false;
3292 :
3293 2 : char *pszWKT = nullptr;
3294 2 : poGeom->exportToWkt(&pszWKT);
3295 2 : CPLDebug("GDALWARP", "WKT = \"%s\"", pszWKT ? pszWKT : "(null)");
3296 : const char *pszFile =
3297 2 : CPLGetConfigOption("GDALWARP_DUMP_WKT_TO_FILE", nullptr);
3298 2 : if (pszFile && pszWKT)
3299 : {
3300 : FILE *f =
3301 0 : EQUAL(pszFile, "stderr") ? stderr : fopen(pszFile, "wb");
3302 0 : if (f)
3303 : {
3304 0 : fprintf(f, "id,WKT\n");
3305 0 : fprintf(f, "1,\"%s\"\n", pszWKT);
3306 0 : if (!EQUAL(pszFile, "stderr"))
3307 0 : fclose(f);
3308 : }
3309 : }
3310 2 : CPLFree(pszWKT);
3311 :
3312 2 : if (CPLTestBool(
3313 : CPLGetConfigOption("GDALWARP_IGNORE_BAD_CUTLINE", "NO")))
3314 0 : CPLError(CE_Warning, CPLE_AppDefined,
3315 : "Cutline polygon is invalid.");
3316 : else
3317 : {
3318 2 : CPLError(CE_Failure, CPLE_AppDefined,
3319 : "Cutline polygon is invalid.");
3320 2 : return false;
3321 : }
3322 : }
3323 : }
3324 : else
3325 : {
3326 1 : if (bVerbose)
3327 : {
3328 1 : CPLError(CE_Failure, CPLE_AppDefined,
3329 : "Cutline not of polygon type.");
3330 : }
3331 1 : return false;
3332 : }
3333 :
3334 184 : return true;
3335 : }
3336 :
3337 : /************************************************************************/
3338 : /* LoadCutline() */
3339 : /* */
3340 : /* Load blend cutline from OGR datasource. */
3341 : /************************************************************************/
3342 :
3343 45 : static CPLErr LoadCutline(const std::string &osCutlineDSNameOrWKT,
3344 : const std::string &osSRS, const std::string &osCLayer,
3345 : const std::string &osCWHERE,
3346 : const std::string &osCSQL, OGRGeometryH *phCutlineRet)
3347 :
3348 : {
3349 45 : if (STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "POLYGON(") ||
3350 45 : STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "POLYGON (") ||
3351 134 : STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "MULTIPOLYGON(") ||
3352 44 : STARTS_WITH_CI(osCutlineDSNameOrWKT.c_str(), "MULTIPOLYGON ("))
3353 : {
3354 0 : std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> poSRS;
3355 1 : if (!osSRS.empty())
3356 : {
3357 1 : poSRS.reset(new OGRSpatialReference());
3358 1 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3359 1 : poSRS->SetFromUserInput(osSRS.c_str());
3360 : }
3361 1 : OGRGeometry *poGeom = nullptr;
3362 1 : OGRGeometryFactory::createFromWkt(osCutlineDSNameOrWKT.c_str(),
3363 1 : poSRS.get(), &poGeom);
3364 1 : *phCutlineRet = OGRGeometry::ToHandle(poGeom);
3365 1 : return *phCutlineRet ? CE_None : CE_Failure;
3366 : }
3367 :
3368 : /* -------------------------------------------------------------------- */
3369 : /* Open source vector dataset. */
3370 : /* -------------------------------------------------------------------- */
3371 : auto poDS = std::unique_ptr<GDALDataset>(
3372 88 : GDALDataset::Open(osCutlineDSNameOrWKT.c_str(), GDAL_OF_VECTOR));
3373 44 : if (poDS == nullptr)
3374 : {
3375 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot open %s.",
3376 : osCutlineDSNameOrWKT.c_str());
3377 1 : return CE_Failure;
3378 : }
3379 :
3380 : /* -------------------------------------------------------------------- */
3381 : /* Get the source layer */
3382 : /* -------------------------------------------------------------------- */
3383 43 : OGRLayer *poLayer = nullptr;
3384 :
3385 43 : if (!osCSQL.empty())
3386 2 : poLayer = poDS->ExecuteSQL(osCSQL.c_str(), nullptr, nullptr);
3387 41 : else if (!osCLayer.empty())
3388 15 : poLayer = poDS->GetLayerByName(osCLayer.c_str());
3389 : else
3390 26 : poLayer = poDS->GetLayer(0);
3391 :
3392 43 : if (poLayer == nullptr)
3393 : {
3394 1 : CPLError(CE_Failure, CPLE_AppDefined,
3395 : "Failed to identify source layer from datasource.");
3396 1 : return CE_Failure;
3397 : }
3398 :
3399 : /* -------------------------------------------------------------------- */
3400 : /* Apply WHERE clause if there is one. */
3401 : /* -------------------------------------------------------------------- */
3402 42 : if (!osCWHERE.empty())
3403 1 : poLayer->SetAttributeFilter(osCWHERE.c_str());
3404 :
3405 : /* -------------------------------------------------------------------- */
3406 : /* Collect the geometries from this layer, and build list of */
3407 : /* burn values. */
3408 : /* -------------------------------------------------------------------- */
3409 84 : auto poMultiPolygon = std::make_unique<OGRMultiPolygon>();
3410 :
3411 81 : for (auto &&poFeature : poLayer)
3412 : {
3413 42 : auto poGeom = std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
3414 42 : if (poGeom == nullptr)
3415 : {
3416 1 : CPLError(CE_Failure, CPLE_AppDefined,
3417 : "Cutline feature without a geometry.");
3418 1 : goto error;
3419 : }
3420 :
3421 41 : if (!ValidateCutline(poGeom.get(), true))
3422 : {
3423 2 : goto error;
3424 : }
3425 :
3426 39 : OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3427 :
3428 39 : if (eType == wkbPolygon)
3429 36 : poMultiPolygon->addGeometry(std::move(poGeom));
3430 3 : else if (eType == wkbMultiPolygon)
3431 : {
3432 7 : for (const auto *poSubGeom : poGeom->toMultiPolygon())
3433 : {
3434 4 : poMultiPolygon->addGeometry(poSubGeom);
3435 : }
3436 : }
3437 : }
3438 :
3439 39 : if (poMultiPolygon->IsEmpty())
3440 : {
3441 1 : CPLError(CE_Failure, CPLE_AppDefined,
3442 : "Did not get any cutline features.");
3443 1 : goto error;
3444 : }
3445 :
3446 : /* -------------------------------------------------------------------- */
3447 : /* Ensure the coordinate system gets set on the geometry. */
3448 : /* -------------------------------------------------------------------- */
3449 38 : if (!osSRS.empty())
3450 : {
3451 : std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> poSRS(
3452 2 : new OGRSpatialReference());
3453 1 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3454 1 : poSRS->SetFromUserInput(osSRS.c_str());
3455 1 : poMultiPolygon->assignSpatialReference(poSRS.get());
3456 : }
3457 : else
3458 : {
3459 37 : poMultiPolygon->assignSpatialReference(poLayer->GetSpatialRef());
3460 : }
3461 :
3462 38 : *phCutlineRet = OGRGeometry::ToHandle(poMultiPolygon.release());
3463 :
3464 : /* -------------------------------------------------------------------- */
3465 : /* Cleanup */
3466 : /* -------------------------------------------------------------------- */
3467 38 : if (!osCSQL.empty())
3468 1 : poDS->ReleaseResultSet(poLayer);
3469 :
3470 38 : return CE_None;
3471 :
3472 4 : error:
3473 4 : if (!osCSQL.empty())
3474 1 : poDS->ReleaseResultSet(poLayer);
3475 :
3476 4 : return CE_Failure;
3477 : }
3478 :
3479 : /************************************************************************/
3480 : /* GDALWarpCreateOutput() */
3481 : /* */
3482 : /* Create the output file based on various command line options, */
3483 : /* and the input file. */
3484 : /* If there's just one source file, then *phTransformArg will be */
3485 : /* set in order them to be reused by main function. This saves */
3486 : /* transform recomputation, which can be expensive in the -tps case*/
3487 : /************************************************************************/
3488 :
3489 794 : static GDALDatasetH GDALWarpCreateOutput(
3490 : int nSrcCount, GDALDatasetH *pahSrcDS, const char *pszFilename,
3491 : const char *pszFormat, char **papszTO, CSLConstList papszCreateOptions,
3492 : GDALDataType eDT, void **phTransformArg, bool bSetColorInterpretation,
3493 : GDALWarpAppOptions *psOptions)
3494 :
3495 : {
3496 : GDALDriverH hDriver;
3497 : GDALDatasetH hDstDS;
3498 : void *hTransformArg;
3499 794 : GDALColorTableH hCT = nullptr;
3500 794 : GDALRasterAttributeTableH hRAT = nullptr;
3501 794 : double dfWrkMinX = 0, dfWrkMaxX = 0, dfWrkMinY = 0, dfWrkMaxY = 0;
3502 794 : double dfWrkResX = 0, dfWrkResY = 0;
3503 794 : int nDstBandCount = 0;
3504 1588 : std::vector<GDALColorInterp> apeColorInterpretations;
3505 794 : bool bVRT = false;
3506 :
3507 794 : if (EQUAL(pszFormat, "VRT"))
3508 92 : bVRT = true;
3509 :
3510 : // Special case for geographic to Mercator (typically EPSG:4326 to EPSG:3857)
3511 : // where latitudes close to 90 go to infinity
3512 : // We clamp latitudes between ~ -85 and ~ 85 degrees.
3513 794 : const char *pszDstSRS = CSLFetchNameValue(papszTO, "DST_SRS");
3514 794 : if (nSrcCount == 1 && pszDstSRS && psOptions->dfMinX == 0.0 &&
3515 95 : psOptions->dfMinY == 0.0 && psOptions->dfMaxX == 0.0 &&
3516 95 : psOptions->dfMaxY == 0.0)
3517 : {
3518 95 : auto hSrcDS = pahSrcDS[0];
3519 190 : const auto osSrcSRS = GetSrcDSProjection(pahSrcDS[0], papszTO);
3520 190 : OGRSpatialReference oSrcSRS;
3521 95 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3522 95 : oSrcSRS.SetFromUserInput(osSrcSRS.c_str());
3523 190 : OGRSpatialReference oDstSRS;
3524 95 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3525 95 : oDstSRS.SetFromUserInput(pszDstSRS);
3526 95 : const char *pszProjection = oDstSRS.GetAttrValue("PROJECTION");
3527 95 : const char *pszMethod = CSLFetchNameValue(papszTO, "METHOD");
3528 95 : double adfSrcGT[6] = {0};
3529 : // This MAX_LAT values is equivalent to the semi_major_axis * PI
3530 : // easting/northing value only for EPSG:3857, but it is also quite
3531 : // reasonable for other Mercator projections
3532 95 : constexpr double MAX_LAT = 85.0511287798066;
3533 95 : constexpr double EPS = 1e-3;
3534 5 : const auto GetMinLon = [&adfSrcGT]() { return adfSrcGT[0]; };
3535 5 : const auto GetMaxLon = [&adfSrcGT, hSrcDS]()
3536 5 : { return adfSrcGT[0] + adfSrcGT[1] * GDALGetRasterXSize(hSrcDS); };
3537 5 : const auto GetMinLat = [&adfSrcGT, hSrcDS]()
3538 5 : { return adfSrcGT[3] + adfSrcGT[5] * GDALGetRasterYSize(hSrcDS); };
3539 6 : const auto GetMaxLat = [&adfSrcGT]() { return adfSrcGT[3]; };
3540 123 : if (oSrcSRS.IsGeographic() && !oSrcSRS.IsDerivedGeographic() &&
3541 25 : pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) &&
3542 3 : oDstSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0) == 0 &&
3543 0 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
3544 3 : CSLFetchNameValue(papszTO, "COORDINATE_OPERATION") == nullptr &&
3545 3 : CSLFetchNameValue(papszTO, "SRC_METHOD") == nullptr &&
3546 3 : CSLFetchNameValue(papszTO, "DST_METHOD") == nullptr &&
3547 3 : GDALGetGeoTransform(hSrcDS, adfSrcGT) == CE_None &&
3548 6 : adfSrcGT[2] == 0 && adfSrcGT[4] == 0 && adfSrcGT[5] < 0 &&
3549 9 : GetMinLon() >= -180 - EPS && GetMaxLon() <= 180 + EPS &&
3550 6 : ((GetMaxLat() > MAX_LAT && GetMinLat() < MAX_LAT) ||
3551 2 : (GetMaxLat() > -MAX_LAT && GetMinLat() < -MAX_LAT)) &&
3552 125 : GDALGetMetadata(hSrcDS, "GEOLOC_ARRAY") == nullptr &&
3553 2 : GDALGetMetadata(hSrcDS, "RPC") == nullptr)
3554 : {
3555 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
3556 4 : OGRCreateCoordinateTransformation(&oSrcSRS, &oDstSRS));
3557 2 : if (poCT)
3558 : {
3559 2 : double xLL = std::max(GetMinLon(), -180.0);
3560 2 : double yLL = std::max(GetMinLat(), -MAX_LAT);
3561 2 : double xUR = std::min(GetMaxLon(), 180.0);
3562 2 : double yUR = std::min(GetMaxLat(), MAX_LAT);
3563 4 : if (poCT->Transform(1, &xLL, &yLL) &&
3564 2 : poCT->Transform(1, &xUR, &yUR))
3565 : {
3566 2 : psOptions->dfMinX = xLL;
3567 2 : psOptions->dfMinY = yLL;
3568 2 : psOptions->dfMaxX = xUR;
3569 2 : psOptions->dfMaxY = yUR;
3570 2 : CPLError(CE_Warning, CPLE_AppDefined,
3571 : "Clamping output bounds to (%f,%f) -> (%f, %f)",
3572 : psOptions->dfMinX, psOptions->dfMinY,
3573 : psOptions->dfMaxX, psOptions->dfMaxY);
3574 : }
3575 : }
3576 : }
3577 : }
3578 :
3579 : /* If (-ts and -te) or (-tr and -te) are specified, we don't need to compute
3580 : * the suggested output extent */
3581 794 : const bool bNeedsSuggestedWarpOutput =
3582 1588 : !(((psOptions->nForcePixels != 0 && psOptions->nForceLines != 0) ||
3583 662 : (psOptions->dfXRes != 0 && psOptions->dfYRes != 0)) &&
3584 170 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
3585 106 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0));
3586 :
3587 : // If -te is specified, not not -tr and -ts
3588 794 : const bool bKnownTargetExtentButNotResolution =
3589 657 : !(psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
3590 655 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0) &&
3591 139 : psOptions->nForcePixels == 0 && psOptions->nForceLines == 0 &&
3592 1588 : psOptions->dfXRes == 0 && psOptions->dfYRes == 0;
3593 :
3594 794 : if (phTransformArg)
3595 791 : *phTransformArg = nullptr;
3596 :
3597 : /* -------------------------------------------------------------------- */
3598 : /* Find the output driver. */
3599 : /* -------------------------------------------------------------------- */
3600 794 : hDriver = GDALGetDriverByName(pszFormat);
3601 1588 : if (hDriver == nullptr ||
3602 794 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) == nullptr &&
3603 0 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) ==
3604 : nullptr))
3605 : {
3606 0 : printf("Output driver `%s' not recognised or does not support\n",
3607 : pszFormat);
3608 0 : printf("direct output file creation or CreateCopy. "
3609 : "The following format drivers are eligible for warp output:\n");
3610 :
3611 0 : for (int iDr = 0; iDr < GDALGetDriverCount(); iDr++)
3612 : {
3613 0 : hDriver = GDALGetDriver(iDr);
3614 :
3615 0 : if (GDALGetMetadataItem(hDriver, GDAL_DCAP_RASTER, nullptr) !=
3616 0 : nullptr &&
3617 0 : (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) !=
3618 0 : nullptr ||
3619 0 : GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, nullptr) !=
3620 : nullptr))
3621 : {
3622 0 : printf(" %s: %s\n", GDALGetDriverShortName(hDriver),
3623 : GDALGetDriverLongName(hDriver));
3624 : }
3625 : }
3626 0 : printf("\n");
3627 0 : return nullptr;
3628 : }
3629 :
3630 : /* -------------------------------------------------------------------- */
3631 : /* For virtual output files, we have to set a special subclass */
3632 : /* of dataset to create. */
3633 : /* -------------------------------------------------------------------- */
3634 1588 : CPLStringList aosCreateOptions(CSLDuplicate(papszCreateOptions));
3635 794 : if (bVRT)
3636 92 : aosCreateOptions.SetNameValue("SUBCLASS", "VRTWarpedDataset");
3637 :
3638 : /* -------------------------------------------------------------------- */
3639 : /* Loop over all input files to collect extents. */
3640 : /* -------------------------------------------------------------------- */
3641 1588 : CPLString osThisTargetSRS;
3642 : {
3643 794 : const char *pszThisTargetSRS = CSLFetchNameValue(papszTO, "DST_SRS");
3644 794 : if (pszThisTargetSRS != nullptr)
3645 178 : osThisTargetSRS = pszThisTargetSRS;
3646 : }
3647 :
3648 1588 : CPLStringList aoTOList(papszTO, FALSE);
3649 :
3650 794 : double dfResFromSourceAndTargetExtent =
3651 : std::numeric_limits<double>::infinity();
3652 :
3653 : /* -------------------------------------------------------------------- */
3654 : /* Establish list of files of output dataset if it already exists. */
3655 : /* -------------------------------------------------------------------- */
3656 1588 : std::set<std::string> oSetExistingDestFiles;
3657 : {
3658 794 : CPLPushErrorHandler(CPLQuietErrorHandler);
3659 794 : const char *const apszAllowedDrivers[] = {pszFormat, nullptr};
3660 : auto poExistingOutputDS = std::unique_ptr<GDALDataset>(
3661 1588 : GDALDataset::Open(pszFilename, GDAL_OF_RASTER, apszAllowedDrivers));
3662 794 : if (poExistingOutputDS)
3663 : {
3664 72 : for (const char *pszFilenameInList :
3665 68 : CPLStringList(poExistingOutputDS->GetFileList()))
3666 : {
3667 : oSetExistingDestFiles.insert(
3668 36 : CPLString(pszFilenameInList).replaceAll('\\', '/'));
3669 : }
3670 : }
3671 794 : CPLPopErrorHandler();
3672 : }
3673 1588 : std::set<std::string> oSetExistingDestFilesFoundInSource;
3674 :
3675 1606 : for (int iSrc = 0; iSrc < nSrcCount; iSrc++)
3676 : {
3677 : /* --------------------------------------------------------------------
3678 : */
3679 : /* Check that there's at least one raster band */
3680 : /* --------------------------------------------------------------------
3681 : */
3682 818 : GDALDatasetH hSrcDS = pahSrcDS[iSrc];
3683 818 : if (GDALGetRasterCount(hSrcDS) == 0)
3684 : {
3685 1 : CPLError(CE_Failure, CPLE_AppDefined,
3686 : "Input file %s has no raster bands.",
3687 : GDALGetDescription(hSrcDS));
3688 1 : if (hCT != nullptr)
3689 1 : GDALDestroyColorTable(hCT);
3690 6 : return nullptr;
3691 : }
3692 :
3693 : // Examine desired overview level and retrieve the corresponding dataset
3694 : // if it exists.
3695 0 : std::unique_ptr<GDALDataset> oDstDSOverview;
3696 817 : if (psOptions->nOvLevel >= 0)
3697 : {
3698 5 : oDstDSOverview.reset(GDALCreateOverviewDataset(
3699 : GDALDataset::FromHandle(hSrcDS), psOptions->nOvLevel,
3700 : /* bThisLevelOnly = */ true));
3701 5 : if (oDstDSOverview)
3702 4 : hSrcDS = oDstDSOverview.get();
3703 : }
3704 :
3705 : /* --------------------------------------------------------------------
3706 : */
3707 : /* Check if the source dataset shares some files with the dest
3708 : * one.*/
3709 : /* --------------------------------------------------------------------
3710 : */
3711 817 : if (!oSetExistingDestFiles.empty())
3712 : {
3713 : // We need to reopen in a temporary dataset for the particular
3714 : // case of overwritten a .tif.ovr file from a .tif
3715 : // If we probe the file list of the .tif, it will then open the
3716 : // .tif.ovr !
3717 35 : auto poSrcDS = GDALDataset::FromHandle(hSrcDS);
3718 35 : const char *const apszAllowedDrivers[] = {
3719 35 : poSrcDS->GetDriver() ? poSrcDS->GetDriver()->GetDescription()
3720 : : nullptr,
3721 35 : nullptr};
3722 : auto poSrcDSTmp = std::unique_ptr<GDALDataset>(GDALDataset::Open(
3723 70 : poSrcDS->GetDescription(), GDAL_OF_RASTER, apszAllowedDrivers));
3724 35 : if (poSrcDSTmp)
3725 : {
3726 37 : for (const char *pszFilenameInList :
3727 72 : CPLStringList(poSrcDSTmp->GetFileList()))
3728 : {
3729 74 : CPLString osFilename(pszFilenameInList);
3730 37 : osFilename.replaceAll('\\', '/');
3731 37 : if (oSetExistingDestFiles.find(osFilename) !=
3732 74 : oSetExistingDestFiles.end())
3733 : {
3734 5 : oSetExistingDestFilesFoundInSource.insert(osFilename);
3735 : }
3736 : }
3737 : }
3738 : }
3739 :
3740 817 : if (eDT == GDT_Unknown)
3741 760 : eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));
3742 :
3743 : /* --------------------------------------------------------------------
3744 : */
3745 : /* If we are processing the first file, and it has a raster */
3746 : /* attribute table, then we will copy it to the destination file.
3747 : */
3748 : /* --------------------------------------------------------------------
3749 : */
3750 817 : if (iSrc == 0)
3751 : {
3752 791 : hRAT = GDALGetDefaultRAT(GDALGetRasterBand(hSrcDS, 1));
3753 791 : if (hRAT != nullptr)
3754 : {
3755 0 : if (psOptions->eResampleAlg != GRA_NearestNeighbour &&
3756 0 : psOptions->eResampleAlg != GRA_Mode &&
3757 0 : GDALRATGetTableType(hRAT) == GRTT_THEMATIC)
3758 : {
3759 0 : if (!psOptions->bQuiet)
3760 : {
3761 0 : CPLError(CE_Warning, CPLE_AppDefined,
3762 : "Warning: Input file %s has a thematic RAT, "
3763 : "which will likely lead "
3764 : "to bad results when using a resampling "
3765 : "method other than nearest neighbour "
3766 : "or mode so we are discarding it.\n",
3767 : GDALGetDescription(hSrcDS));
3768 : }
3769 0 : hRAT = nullptr;
3770 : }
3771 : else
3772 : {
3773 0 : if (!psOptions->bQuiet)
3774 0 : printf("Copying raster attribute table from %s to new "
3775 : "file.\n",
3776 : GDALGetDescription(hSrcDS));
3777 : }
3778 : }
3779 : }
3780 :
3781 : /* --------------------------------------------------------------------
3782 : */
3783 : /* If we are processing the first file, and it has a color */
3784 : /* table, then we will copy it to the destination file. */
3785 : /* --------------------------------------------------------------------
3786 : */
3787 817 : if (iSrc == 0)
3788 : {
3789 791 : hCT = GDALGetRasterColorTable(GDALGetRasterBand(hSrcDS, 1));
3790 791 : if (hCT != nullptr)
3791 : {
3792 5 : hCT = GDALCloneColorTable(hCT);
3793 5 : if (!psOptions->bQuiet)
3794 0 : printf("Copying color table from %s to new file.\n",
3795 : GDALGetDescription(hSrcDS));
3796 : }
3797 :
3798 791 : if (psOptions->anDstBands.empty())
3799 : {
3800 774 : nDstBandCount = GDALGetRasterCount(hSrcDS);
3801 1786 : for (int iBand = 0; iBand < nDstBandCount; iBand++)
3802 : {
3803 1012 : if (psOptions->anDstBands.empty())
3804 : {
3805 : GDALColorInterp eInterp =
3806 1012 : GDALGetRasterColorInterpretation(
3807 1012 : GDALGetRasterBand(hSrcDS, iBand + 1));
3808 1012 : apeColorInterpretations.push_back(eInterp);
3809 : }
3810 : }
3811 :
3812 : // Do we want to generate an alpha band in the output file?
3813 774 : if (psOptions->bEnableSrcAlpha)
3814 16 : nDstBandCount--;
3815 :
3816 774 : if (psOptions->bEnableDstAlpha)
3817 54 : nDstBandCount++;
3818 : }
3819 : else
3820 : {
3821 45 : for (int nSrcBand : psOptions->anSrcBands)
3822 : {
3823 28 : auto hBand = GDALGetRasterBand(hSrcDS, nSrcBand);
3824 : GDALColorInterp eInterp =
3825 28 : hBand ? GDALGetRasterColorInterpretation(hBand)
3826 28 : : GCI_Undefined;
3827 28 : apeColorInterpretations.push_back(eInterp);
3828 : }
3829 17 : nDstBandCount = static_cast<int>(psOptions->anDstBands.size());
3830 17 : if (psOptions->bEnableDstAlpha)
3831 : {
3832 4 : nDstBandCount++;
3833 4 : apeColorInterpretations.push_back(GCI_AlphaBand);
3834 : }
3835 13 : else if (GDALGetRasterCount(hSrcDS) &&
3836 13 : GDALGetRasterColorInterpretation(GDALGetRasterBand(
3837 : hSrcDS, GDALGetRasterCount(hSrcDS))) ==
3838 26 : GCI_AlphaBand &&
3839 1 : !psOptions->bDisableSrcAlpha)
3840 : {
3841 0 : nDstBandCount++;
3842 0 : apeColorInterpretations.push_back(GCI_AlphaBand);
3843 : }
3844 : }
3845 : }
3846 :
3847 : /* --------------------------------------------------------------------
3848 : */
3849 : /* If we are processing the first file, get the source srs from */
3850 : /* dataset, if not set already. */
3851 : /* --------------------------------------------------------------------
3852 : */
3853 817 : const auto osThisSourceSRS = GetSrcDSProjection(hSrcDS, papszTO);
3854 817 : if (iSrc == 0 && osThisTargetSRS.empty())
3855 : {
3856 613 : if (!osThisSourceSRS.empty())
3857 : {
3858 532 : osThisTargetSRS = osThisSourceSRS;
3859 532 : aoTOList.SetNameValue("DST_SRS", osThisSourceSRS);
3860 : }
3861 : }
3862 :
3863 : /* --------------------------------------------------------------------
3864 : */
3865 : /* Create a transformation object from the source to */
3866 : /* destination coordinate system. */
3867 : /* --------------------------------------------------------------------
3868 : */
3869 : hTransformArg =
3870 817 : GDALCreateGenImgProjTransformer2(hSrcDS, nullptr, aoTOList.List());
3871 :
3872 817 : if (hTransformArg == nullptr)
3873 : {
3874 5 : if (hCT != nullptr)
3875 1 : GDALDestroyColorTable(hCT);
3876 5 : return nullptr;
3877 : }
3878 :
3879 812 : GDALTransformerInfo *psInfo =
3880 : static_cast<GDALTransformerInfo *>(hTransformArg);
3881 :
3882 : /* --------------------------------------------------------------------
3883 : */
3884 : /* Get approximate output resolution */
3885 : /* --------------------------------------------------------------------
3886 : */
3887 :
3888 812 : if (bKnownTargetExtentButNotResolution)
3889 : {
3890 : // Sample points along a grid in target CRS
3891 82 : constexpr int nPointsX = 10;
3892 82 : constexpr int nPointsY = 10;
3893 82 : constexpr int nPoints = 3 * nPointsX * nPointsY;
3894 164 : std::vector<double> padfX;
3895 164 : std::vector<double> padfY;
3896 164 : std::vector<double> padfZ(nPoints);
3897 164 : std::vector<int> pabSuccess(nPoints);
3898 : const double dfEps =
3899 164 : std::min(psOptions->dfMaxX - psOptions->dfMinX,
3900 82 : std::abs(psOptions->dfMaxY - psOptions->dfMinY)) /
3901 82 : 1000;
3902 902 : for (int iY = 0; iY < nPointsY; iY++)
3903 : {
3904 9020 : for (int iX = 0; iX < nPointsX; iX++)
3905 : {
3906 8200 : const double dfX =
3907 8200 : psOptions->dfMinX +
3908 8200 : static_cast<double>(iX) *
3909 8200 : (psOptions->dfMaxX - psOptions->dfMinX) /
3910 : (nPointsX - 1);
3911 8200 : const double dfY =
3912 8200 : psOptions->dfMinY +
3913 8200 : static_cast<double>(iY) *
3914 8200 : (psOptions->dfMaxY - psOptions->dfMinY) /
3915 : (nPointsY - 1);
3916 :
3917 : // Reproject each destination sample point and its
3918 : // neighbours at (x+1,y) and (x,y+1), so as to get the local
3919 : // scale.
3920 8200 : padfX.push_back(dfX);
3921 8200 : padfY.push_back(dfY);
3922 :
3923 15580 : padfX.push_back((iX == nPointsX - 1) ? dfX - dfEps
3924 7380 : : dfX + dfEps);
3925 8200 : padfY.push_back(dfY);
3926 :
3927 8200 : padfX.push_back(dfX);
3928 15580 : padfY.push_back((iY == nPointsY - 1) ? dfY - dfEps
3929 7380 : : dfY + dfEps);
3930 : }
3931 : }
3932 :
3933 82 : bool transformedToSrcCRS{false};
3934 :
3935 82 : GDALGenImgProjTransformInfo *psTransformInfo{
3936 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg)};
3937 :
3938 : // If a transformer is available, use an extent that covers the
3939 : // target extent instead of the real source image extent, but also
3940 : // check for target extent compatibility with source CRS extent
3941 82 : if (psTransformInfo && psTransformInfo->pReprojectArg &&
3942 67 : psTransformInfo->pSrcTransformer == nullptr)
3943 : {
3944 67 : const GDALReprojectionTransformInfo *psRTI =
3945 : static_cast<const GDALReprojectionTransformInfo *>(
3946 : psTransformInfo->pReprojectArg);
3947 67 : if (psRTI && psRTI->poReverseTransform)
3948 : {
3949 :
3950 : // Compute new geotransform from transformed target extent
3951 : double adfGeoTransform[6];
3952 67 : if (GDALGetGeoTransform(hSrcDS, adfGeoTransform) ==
3953 67 : CE_None &&
3954 67 : adfGeoTransform[2] == 0 && adfGeoTransform[4] == 0)
3955 : {
3956 :
3957 : // Transform target extent to source CRS
3958 67 : double dfMinX = psOptions->dfMinX;
3959 67 : double dfMinY = psOptions->dfMinY;
3960 :
3961 : // Need this to check if the target extent is compatible with the source extent
3962 67 : double dfMaxX = psOptions->dfMaxX;
3963 67 : double dfMaxY = psOptions->dfMaxY;
3964 :
3965 : // Clone of psRTI->poReverseTransform with CHECK_WITH_INVERT_PROJ set to TRUE
3966 : // to detect out of source CRS bounds destination extent and fall back to original
3967 : // algorithm if needed
3968 : CPLConfigOptionSetter oSetter("CHECK_WITH_INVERT_PROJ",
3969 134 : "TRUE", false);
3970 134 : OGRCoordinateTransformationOptions options;
3971 : auto poReverseTransform =
3972 : std::unique_ptr<OGRCoordinateTransformation>(
3973 : OGRCreateCoordinateTransformation(
3974 67 : psRTI->poReverseTransform->GetSourceCS(),
3975 67 : psRTI->poReverseTransform->GetTargetCS(),
3976 134 : options));
3977 :
3978 67 : if (poReverseTransform)
3979 : {
3980 :
3981 134 : poReverseTransform->Transform(
3982 67 : 1, &dfMinX, &dfMinY, nullptr, &pabSuccess[0]);
3983 :
3984 67 : if (pabSuccess[0])
3985 : {
3986 65 : adfGeoTransform[0] = dfMinX;
3987 65 : adfGeoTransform[3] = dfMinY;
3988 :
3989 130 : poReverseTransform->Transform(1, &dfMaxX,
3990 : &dfMaxY, nullptr,
3991 65 : &pabSuccess[0]);
3992 :
3993 65 : if (pabSuccess[0])
3994 : {
3995 :
3996 : // Reproject to source image CRS
3997 64 : psRTI->poReverseTransform->Transform(
3998 64 : nPoints, &padfX[0], &padfY[0],
3999 64 : &padfZ[0], &pabSuccess[0]);
4000 :
4001 : // Transform back to source image coordinate space using geotransform
4002 19264 : for (size_t i = 0; i < padfX.size(); i++)
4003 : {
4004 38400 : padfX[i] =
4005 19200 : (padfX[i] - adfGeoTransform[0]) /
4006 19200 : adfGeoTransform[1];
4007 19200 : padfY[i] = std::abs(
4008 19200 : (padfY[i] - adfGeoTransform[3]) /
4009 19200 : adfGeoTransform[5]);
4010 : }
4011 :
4012 64 : transformedToSrcCRS = true;
4013 : }
4014 : }
4015 : }
4016 : }
4017 : }
4018 : }
4019 :
4020 82 : if (!transformedToSrcCRS)
4021 : {
4022 : // Transform to source image coordinate space
4023 18 : psInfo->pfnTransform(hTransformArg, TRUE, nPoints, &padfX[0],
4024 18 : &padfY[0], &padfZ[0], &pabSuccess[0]);
4025 : }
4026 :
4027 : // Compute the resolution at sampling points
4028 164 : std::vector<std::pair<double, double>> aoResPairs;
4029 :
4030 15144 : const auto Distance = [](double x, double y)
4031 15144 : { return sqrt(x * x + y * y); };
4032 :
4033 82 : const int nSrcXSize = GDALGetRasterXSize(hSrcDS);
4034 82 : const int nSrcYSize = GDALGetRasterYSize(hSrcDS);
4035 :
4036 8282 : for (int i = 0; i < nPoints; i += 3)
4037 : {
4038 16400 : if (pabSuccess[i] && pabSuccess[i + 1] && pabSuccess[i + 2] &&
4039 17963 : padfX[i] >= 0 && padfY[i] >= 0 &&
4040 1563 : (transformedToSrcCRS ||
4041 1563 : (padfX[i] <= nSrcXSize && padfY[i] <= nSrcYSize)))
4042 : {
4043 : const double dfRes1 =
4044 15144 : std::abs(dfEps) / Distance(padfX[i + 1] - padfX[i],
4045 7572 : padfY[i + 1] - padfY[i]);
4046 : const double dfRes2 =
4047 15144 : std::abs(dfEps) / Distance(padfX[i + 2] - padfX[i],
4048 7572 : padfY[i + 2] - padfY[i]);
4049 7572 : if (std::isfinite(dfRes1) && std::isfinite(dfRes2))
4050 : {
4051 7552 : aoResPairs.push_back(
4052 15104 : std::pair<double, double>(dfRes1, dfRes2));
4053 : }
4054 : }
4055 : }
4056 :
4057 : // Find the minimum resolution that is at least 10 times greater
4058 : // than the median, to remove outliers.
4059 : // Start first by doing that on dfRes1, then dfRes2 and then their
4060 : // average.
4061 82 : std::sort(aoResPairs.begin(), aoResPairs.end(),
4062 40075 : [](const std::pair<double, double> &oPair1,
4063 : const std::pair<double, double> &oPair2)
4064 40075 : { return oPair1.first < oPair2.first; });
4065 :
4066 82 : if (!aoResPairs.empty())
4067 : {
4068 164 : std::vector<std::pair<double, double>> aoResPairsNew;
4069 : const double dfMedian1 =
4070 82 : aoResPairs[aoResPairs.size() / 2].first;
4071 7634 : for (const auto &oPair : aoResPairs)
4072 : {
4073 7552 : if (oPair.first > dfMedian1 / 10)
4074 : {
4075 7532 : aoResPairsNew.push_back(oPair);
4076 : }
4077 : }
4078 :
4079 82 : aoResPairs = std::move(aoResPairsNew);
4080 82 : std::sort(aoResPairs.begin(), aoResPairs.end(),
4081 42866 : [](const std::pair<double, double> &oPair1,
4082 : const std::pair<double, double> &oPair2)
4083 42866 : { return oPair1.second < oPair2.second; });
4084 82 : if (!aoResPairs.empty())
4085 : {
4086 164 : std::vector<double> adfRes;
4087 : const double dfMedian2 =
4088 82 : aoResPairs[aoResPairs.size() / 2].second;
4089 7614 : for (const auto &oPair : aoResPairs)
4090 : {
4091 7532 : if (oPair.second > dfMedian2 / 10)
4092 : {
4093 7532 : adfRes.push_back((oPair.first + oPair.second) / 2);
4094 : }
4095 : }
4096 :
4097 82 : std::sort(adfRes.begin(), adfRes.end());
4098 82 : if (!adfRes.empty())
4099 : {
4100 82 : const double dfMedian = adfRes[adfRes.size() / 2];
4101 82 : for (const double dfRes : adfRes)
4102 : {
4103 82 : if (dfRes > dfMedian / 10)
4104 : {
4105 82 : dfResFromSourceAndTargetExtent = std::min(
4106 82 : dfResFromSourceAndTargetExtent, dfRes);
4107 82 : break;
4108 : }
4109 : }
4110 : }
4111 : }
4112 : }
4113 : }
4114 :
4115 : /* --------------------------------------------------------------------
4116 : */
4117 : /* Get approximate output definition. */
4118 : /* --------------------------------------------------------------------
4119 : */
4120 : double adfThisGeoTransform[6];
4121 : double adfExtent[4];
4122 812 : if (bNeedsSuggestedWarpOutput)
4123 : {
4124 : int nThisPixels, nThisLines;
4125 :
4126 : // For sum, round-up dimension, to be sure that the output extent
4127 : // includes all source pixels, to have the sum preserving property.
4128 1494 : int nOptions = (psOptions->eResampleAlg == GRA_Sum)
4129 747 : ? GDAL_SWO_ROUND_UP_SIZE
4130 : : 0;
4131 747 : if (psOptions->bSquarePixels)
4132 : {
4133 1 : nOptions |= GDAL_SWO_FORCE_SQUARE_PIXEL;
4134 : }
4135 :
4136 747 : if (GDALSuggestedWarpOutput2(hSrcDS, psInfo->pfnTransform,
4137 : hTransformArg, adfThisGeoTransform,
4138 : &nThisPixels, &nThisLines, adfExtent,
4139 747 : nOptions) != CE_None)
4140 : {
4141 0 : if (hCT != nullptr)
4142 0 : GDALDestroyColorTable(hCT);
4143 0 : GDALDestroyGenImgProjTransformer(hTransformArg);
4144 0 : return nullptr;
4145 : }
4146 :
4147 747 : if (CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", nullptr) ==
4148 : nullptr)
4149 : {
4150 747 : double MinX = adfExtent[0];
4151 747 : double MaxX = adfExtent[2];
4152 747 : double MaxY = adfExtent[3];
4153 747 : double MinY = adfExtent[1];
4154 747 : int bSuccess = TRUE;
4155 :
4156 : // +/-180 deg in longitude do not roundtrip sometimes
4157 747 : if (MinX == -180)
4158 35 : MinX += 1e-6;
4159 747 : if (MaxX == 180)
4160 34 : MaxX -= 1e-6;
4161 :
4162 : // +/-90 deg in latitude do not roundtrip sometimes
4163 747 : if (MinY == -90)
4164 24 : MinY += 1e-6;
4165 747 : if (MaxY == 90)
4166 38 : MaxY -= 1e-6;
4167 :
4168 : /* Check that the edges of the target image are in the validity
4169 : * area */
4170 : /* of the target projection */
4171 747 : const int N_STEPS = 20;
4172 15964 : for (int i = 0; i <= N_STEPS && bSuccess; i++)
4173 : {
4174 334378 : for (int j = 0; j <= N_STEPS && bSuccess; j++)
4175 : {
4176 319161 : const double dfRatioI = i * 1.0 / N_STEPS;
4177 319161 : const double dfRatioJ = j * 1.0 / N_STEPS;
4178 319161 : const double expected_x =
4179 319161 : (1 - dfRatioI) * MinX + dfRatioI * MaxX;
4180 319161 : const double expected_y =
4181 319161 : (1 - dfRatioJ) * MinY + dfRatioJ * MaxY;
4182 319161 : double x = expected_x;
4183 319161 : double y = expected_y;
4184 319161 : double z = 0;
4185 : /* Target SRS coordinates to source image pixel
4186 : * coordinates */
4187 319161 : if (!psInfo->pfnTransform(hTransformArg, TRUE, 1, &x,
4188 638308 : &y, &z, &bSuccess) ||
4189 319147 : !bSuccess)
4190 14 : bSuccess = FALSE;
4191 : /* Source image pixel coordinates to target SRS
4192 : * coordinates */
4193 319161 : if (!psInfo->pfnTransform(hTransformArg, FALSE, 1, &x,
4194 638307 : &y, &z, &bSuccess) ||
4195 319146 : !bSuccess)
4196 15 : bSuccess = FALSE;
4197 319161 : if (fabs(x - expected_x) >
4198 319161 : (MaxX - MinX) / nThisPixels ||
4199 319137 : fabs(y - expected_y) > (MaxY - MinY) / nThisLines)
4200 24 : bSuccess = FALSE;
4201 : }
4202 : }
4203 :
4204 : /* If not, retry with CHECK_WITH_INVERT_PROJ=TRUE that forces
4205 : * ogrct.cpp */
4206 : /* to check the consistency of each requested projection result
4207 : * with the */
4208 : /* invert projection */
4209 747 : if (!bSuccess)
4210 : {
4211 24 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ",
4212 : "TRUE");
4213 24 : CPLDebug("WARP", "Recompute out extent with "
4214 : "CHECK_WITH_INVERT_PROJ=TRUE");
4215 :
4216 24 : const CPLErr eErr = GDALSuggestedWarpOutput2(
4217 : hSrcDS, psInfo->pfnTransform, hTransformArg,
4218 : adfThisGeoTransform, &nThisPixels, &nThisLines,
4219 : adfExtent, 0);
4220 24 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ",
4221 : nullptr);
4222 24 : if (eErr != CE_None)
4223 : {
4224 0 : if (hCT != nullptr)
4225 0 : GDALDestroyColorTable(hCT);
4226 0 : GDALDestroyGenImgProjTransformer(hTransformArg);
4227 0 : return nullptr;
4228 : }
4229 : }
4230 : }
4231 : }
4232 :
4233 : // If no reprojection or geometry change is involved, and that the
4234 : // source image is north-up, preserve source resolution instead of
4235 : // forcing square pixels.
4236 812 : const char *pszMethod = CSLFetchNameValue(papszTO, "METHOD");
4237 : double adfThisGeoTransformTmp[6];
4238 811 : if (!psOptions->bSquarePixels && bNeedsSuggestedWarpOutput &&
4239 746 : psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
4240 722 : psOptions->nForcePixels == 0 && psOptions->nForceLines == 0 &&
4241 16 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
4242 619 : CSLFetchNameValue(papszTO, "COORDINATE_OPERATION") == nullptr &&
4243 615 : CSLFetchNameValue(papszTO, "SRC_METHOD") == nullptr &&
4244 605 : CSLFetchNameValue(papszTO, "DST_METHOD") == nullptr &&
4245 603 : GDALGetGeoTransform(hSrcDS, adfThisGeoTransformTmp) == CE_None &&
4246 584 : adfThisGeoTransformTmp[2] == 0 && adfThisGeoTransformTmp[4] == 0 &&
4247 584 : adfThisGeoTransformTmp[5] < 0 &&
4248 2185 : GDALGetMetadata(hSrcDS, "GEOLOC_ARRAY") == nullptr &&
4249 562 : GDALGetMetadata(hSrcDS, "RPC") == nullptr)
4250 : {
4251 562 : bool bIsSameHorizontal = osThisSourceSRS == osThisTargetSRS;
4252 562 : if (!bIsSameHorizontal)
4253 : {
4254 218 : OGRSpatialReference oSrcSRS;
4255 218 : OGRSpatialReference oDstSRS;
4256 218 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
4257 : // DemoteTo2D requires PROJ >= 6.3
4258 109 : if (oSrcSRS.SetFromUserInput(osThisSourceSRS.c_str()) ==
4259 109 : OGRERR_NONE &&
4260 109 : oDstSRS.SetFromUserInput(osThisTargetSRS.c_str()) ==
4261 109 : OGRERR_NONE &&
4262 109 : (oSrcSRS.GetAxesCount() == 3 ||
4263 93 : oDstSRS.GetAxesCount() == 3) &&
4264 238 : oSrcSRS.DemoteTo2D(nullptr) == OGRERR_NONE &&
4265 20 : oDstSRS.DemoteTo2D(nullptr) == OGRERR_NONE)
4266 : {
4267 20 : bIsSameHorizontal = oSrcSRS.IsSame(&oDstSRS);
4268 : }
4269 : }
4270 562 : if (bIsSameHorizontal)
4271 : {
4272 459 : memcpy(adfThisGeoTransform, adfThisGeoTransformTmp,
4273 : 6 * sizeof(double));
4274 459 : adfExtent[0] = adfThisGeoTransform[0];
4275 459 : adfExtent[1] =
4276 918 : adfThisGeoTransform[3] +
4277 459 : GDALGetRasterYSize(hSrcDS) * adfThisGeoTransform[5];
4278 459 : adfExtent[2] =
4279 918 : adfThisGeoTransform[0] +
4280 459 : GDALGetRasterXSize(hSrcDS) * adfThisGeoTransform[1];
4281 459 : adfExtent[3] = adfThisGeoTransform[3];
4282 459 : dfResFromSourceAndTargetExtent =
4283 459 : std::numeric_limits<double>::infinity();
4284 : }
4285 : }
4286 :
4287 812 : if (bNeedsSuggestedWarpOutput)
4288 : {
4289 : /* --------------------------------------------------------------------
4290 : */
4291 : /* Expand the working bounds to include this region, ensure the
4292 : */
4293 : /* working resolution is no more than this resolution. */
4294 : /* --------------------------------------------------------------------
4295 : */
4296 747 : if (dfWrkMaxX == 0.0 && dfWrkMinX == 0.0)
4297 : {
4298 723 : dfWrkMinX = adfExtent[0];
4299 723 : dfWrkMaxX = adfExtent[2];
4300 723 : dfWrkMaxY = adfExtent[3];
4301 723 : dfWrkMinY = adfExtent[1];
4302 723 : dfWrkResX = adfThisGeoTransform[1];
4303 723 : dfWrkResY = std::abs(adfThisGeoTransform[5]);
4304 : }
4305 : else
4306 : {
4307 24 : dfWrkMinX = std::min(dfWrkMinX, adfExtent[0]);
4308 24 : dfWrkMaxX = std::max(dfWrkMaxX, adfExtent[2]);
4309 24 : dfWrkMaxY = std::max(dfWrkMaxY, adfExtent[3]);
4310 24 : dfWrkMinY = std::min(dfWrkMinY, adfExtent[1]);
4311 24 : dfWrkResX = std::min(dfWrkResX, adfThisGeoTransform[1]);
4312 24 : dfWrkResY =
4313 24 : std::min(dfWrkResY, std::abs(adfThisGeoTransform[5]));
4314 : }
4315 : }
4316 :
4317 812 : if (nSrcCount == 1 && phTransformArg)
4318 : {
4319 759 : *phTransformArg = hTransformArg;
4320 : }
4321 : else
4322 : {
4323 53 : GDALDestroyGenImgProjTransformer(hTransformArg);
4324 : }
4325 : }
4326 :
4327 : // If the source file(s) and the dest one share some files in common,
4328 : // only remove the files that are *not* in common
4329 788 : if (!oSetExistingDestFilesFoundInSource.empty())
4330 : {
4331 14 : for (const std::string &osFilename : oSetExistingDestFiles)
4332 : {
4333 9 : if (oSetExistingDestFilesFoundInSource.find(osFilename) ==
4334 18 : oSetExistingDestFilesFoundInSource.end())
4335 : {
4336 4 : VSIUnlink(osFilename.c_str());
4337 : }
4338 : }
4339 : }
4340 :
4341 788 : if (std::isfinite(dfResFromSourceAndTargetExtent))
4342 : {
4343 37 : dfWrkResX = dfResFromSourceAndTargetExtent;
4344 37 : dfWrkResY = dfResFromSourceAndTargetExtent;
4345 : }
4346 :
4347 : /* -------------------------------------------------------------------- */
4348 : /* Did we have any usable sources? */
4349 : /* -------------------------------------------------------------------- */
4350 788 : if (nDstBandCount == 0)
4351 : {
4352 3 : CPLError(CE_Failure, CPLE_AppDefined, "No usable source images.");
4353 3 : if (hCT != nullptr)
4354 0 : GDALDestroyColorTable(hCT);
4355 3 : return nullptr;
4356 : }
4357 :
4358 : /* -------------------------------------------------------------------- */
4359 : /* Turn the suggested region into a geotransform and suggested */
4360 : /* number of pixels and lines. */
4361 : /* -------------------------------------------------------------------- */
4362 785 : double adfDstGeoTransform[6] = {0, 0, 0, 0, 0, 0};
4363 785 : int nPixels = 0;
4364 785 : int nLines = 0;
4365 :
4366 785 : if (bNeedsSuggestedWarpOutput)
4367 : {
4368 721 : adfDstGeoTransform[0] = dfWrkMinX;
4369 721 : adfDstGeoTransform[1] = dfWrkResX;
4370 721 : adfDstGeoTransform[2] = 0.0;
4371 721 : adfDstGeoTransform[3] = dfWrkMaxY;
4372 721 : adfDstGeoTransform[4] = 0.0;
4373 721 : adfDstGeoTransform[5] = -1 * dfWrkResY;
4374 :
4375 721 : nPixels = static_cast<int>((dfWrkMaxX - dfWrkMinX) / dfWrkResX + 0.5);
4376 721 : nLines = static_cast<int>((dfWrkMaxY - dfWrkMinY) / dfWrkResY + 0.5);
4377 : }
4378 :
4379 : /* -------------------------------------------------------------------- */
4380 : /* Did the user override some parameters? */
4381 : /* -------------------------------------------------------------------- */
4382 785 : if (UseTEAndTSAndTRConsistently(psOptions))
4383 : {
4384 25 : adfDstGeoTransform[0] = psOptions->dfMinX;
4385 25 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4386 25 : adfDstGeoTransform[1] = psOptions->dfXRes;
4387 25 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4388 :
4389 25 : nPixels = psOptions->nForcePixels;
4390 25 : nLines = psOptions->nForceLines;
4391 : }
4392 760 : else if (psOptions->dfXRes != 0.0 && psOptions->dfYRes != 0.0)
4393 : {
4394 38 : bool bDetectBlankBorders = false;
4395 :
4396 38 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4397 24 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4398 : {
4399 24 : bDetectBlankBorders = bNeedsSuggestedWarpOutput;
4400 :
4401 24 : psOptions->dfMinX = adfDstGeoTransform[0];
4402 24 : psOptions->dfMaxX =
4403 24 : adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
4404 24 : psOptions->dfMaxY = adfDstGeoTransform[3];
4405 24 : psOptions->dfMinY =
4406 24 : adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;
4407 : }
4408 :
4409 72 : if (psOptions->bTargetAlignedPixels ||
4410 34 : (psOptions->bCropToCutline &&
4411 1 : psOptions->aosWarpOptions.FetchBool("CUTLINE_ALL_TOUCHED", false)))
4412 : {
4413 4 : if ((psOptions->bTargetAlignedPixels &&
4414 10 : bNeedsSuggestedWarpOutput) ||
4415 2 : (psOptions->bCropToCutline &&
4416 1 : psOptions->aosWarpOptions.FetchBool("CUTLINE_ALL_TOUCHED",
4417 : false)))
4418 : {
4419 4 : bDetectBlankBorders = true;
4420 : }
4421 5 : constexpr double EPS = 1e-8;
4422 5 : psOptions->dfMinX =
4423 5 : floor(psOptions->dfMinX / psOptions->dfXRes + EPS) *
4424 5 : psOptions->dfXRes;
4425 5 : psOptions->dfMaxX =
4426 5 : ceil(psOptions->dfMaxX / psOptions->dfXRes - EPS) *
4427 5 : psOptions->dfXRes;
4428 5 : psOptions->dfMinY =
4429 5 : floor(psOptions->dfMinY / psOptions->dfYRes + EPS) *
4430 5 : psOptions->dfYRes;
4431 5 : psOptions->dfMaxY =
4432 5 : ceil(psOptions->dfMaxY / psOptions->dfYRes - EPS) *
4433 5 : psOptions->dfYRes;
4434 : }
4435 :
4436 62 : const auto UpdateGeoTransformandAndPixelLines = [&]()
4437 : {
4438 249 : nPixels = static_cast<int>((psOptions->dfMaxX - psOptions->dfMinX +
4439 62 : (psOptions->dfXRes / 2.0)) /
4440 62 : psOptions->dfXRes);
4441 62 : if (nPixels == 0)
4442 1 : nPixels = 1;
4443 63 : nLines = static_cast<int>(
4444 62 : (std::fabs(psOptions->dfMaxY - psOptions->dfMinY) +
4445 62 : (psOptions->dfYRes / 2.0)) /
4446 62 : psOptions->dfYRes);
4447 62 : if (nLines == 0)
4448 1 : nLines = 1;
4449 124 : adfDstGeoTransform[0] = psOptions->dfMinX;
4450 62 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4451 62 : adfDstGeoTransform[1] = psOptions->dfXRes;
4452 124 : adfDstGeoTransform[5] = (psOptions->dfMaxY > psOptions->dfMinY)
4453 62 : ? -psOptions->dfYRes
4454 0 : : psOptions->dfYRes;
4455 100 : };
4456 :
4457 38 : if (bDetectBlankBorders && nSrcCount == 1 && phTransformArg &&
4458 25 : *phTransformArg != nullptr)
4459 : {
4460 : // Try to detect if the edge of the raster would be blank
4461 : // Cf https://github.com/OSGeo/gdal/issues/7905
4462 27 : while (nPixels > 1 || nLines > 1)
4463 : {
4464 24 : UpdateGeoTransformandAndPixelLines();
4465 :
4466 24 : GDALSetGenImgProjTransformerDstGeoTransform(*phTransformArg,
4467 : adfDstGeoTransform);
4468 :
4469 24 : std::vector<double> adfX(std::max(nPixels, nLines));
4470 24 : std::vector<double> adfY(adfX.size());
4471 24 : std::vector<double> adfZ(adfX.size());
4472 24 : std::vector<int> abSuccess(adfX.size());
4473 :
4474 : const auto DetectBlankBorder =
4475 96 : [&](int nValues,
4476 : std::function<bool(double, double)> funcIsOK)
4477 : {
4478 96 : if (nValues > 3)
4479 : {
4480 : // First try with just a subsample of 3 points
4481 40 : double adf3X[3] = {adfX[0], adfX[nValues / 2],
4482 40 : adfX[nValues - 1]};
4483 40 : double adf3Y[3] = {adfY[0], adfY[nValues / 2],
4484 40 : adfY[nValues - 1]};
4485 40 : double adf3Z[3] = {0};
4486 40 : if (GDALGenImgProjTransform(*phTransformArg, TRUE, 3,
4487 : &adf3X[0], &adf3Y[0],
4488 80 : &adf3Z[0], &abSuccess[0]))
4489 : {
4490 49 : for (int i = 0; i < 3; ++i)
4491 : {
4492 92 : if (abSuccess[i] &&
4493 46 : funcIsOK(adf3X[i], adf3Y[i]))
4494 : {
4495 37 : return false;
4496 : }
4497 : }
4498 : }
4499 : }
4500 :
4501 : // Do on full border to confirm
4502 59 : if (GDALGenImgProjTransform(*phTransformArg, TRUE, nValues,
4503 59 : &adfX[0], &adfY[0], &adfZ[0],
4504 118 : &abSuccess[0]))
4505 : {
4506 125 : for (int i = 0; i < nValues; ++i)
4507 : {
4508 122 : if (abSuccess[i] && funcIsOK(adfX[i], adfY[i]))
4509 : {
4510 56 : return false;
4511 : }
4512 : }
4513 : }
4514 :
4515 3 : return true;
4516 24 : };
4517 :
4518 854 : for (int i = 0; i < nPixels; ++i)
4519 : {
4520 830 : adfX[i] = i + 0.5;
4521 830 : adfY[i] = 0.5;
4522 830 : adfZ[i] = 0;
4523 : }
4524 24 : const bool bTopBlankLine = DetectBlankBorder(
4525 39 : nPixels, [](double, double y) { return y >= 0; });
4526 :
4527 854 : for (int i = 0; i < nPixels; ++i)
4528 : {
4529 830 : adfX[i] = i + 0.5;
4530 830 : adfY[i] = nLines - 0.5;
4531 830 : adfZ[i] = 0;
4532 : }
4533 24 : const int nSrcLines = GDALGetRasterYSize(pahSrcDS[0]);
4534 : const bool bBottomBlankLine =
4535 24 : DetectBlankBorder(nPixels, [nSrcLines](double, double y)
4536 24 : { return y <= nSrcLines; });
4537 :
4538 672 : for (int i = 0; i < nLines; ++i)
4539 : {
4540 648 : adfX[i] = 0.5;
4541 648 : adfY[i] = i + 0.5;
4542 648 : adfZ[i] = 0;
4543 : }
4544 24 : const bool bLeftBlankCol = DetectBlankBorder(
4545 54 : nLines, [](double x, double) { return x >= 0; });
4546 :
4547 672 : for (int i = 0; i < nLines; ++i)
4548 : {
4549 648 : adfX[i] = nPixels - 0.5;
4550 648 : adfY[i] = i + 0.5;
4551 648 : adfZ[i] = 0;
4552 : }
4553 24 : const int nSrcCols = GDALGetRasterXSize(pahSrcDS[0]);
4554 : const bool bRightBlankCol =
4555 24 : DetectBlankBorder(nLines, [nSrcCols](double x, double)
4556 51 : { return x <= nSrcCols; });
4557 :
4558 24 : if (!bTopBlankLine && !bBottomBlankLine && !bLeftBlankCol &&
4559 22 : !bRightBlankCol)
4560 22 : break;
4561 :
4562 2 : if (bTopBlankLine)
4563 : {
4564 1 : if (psOptions->dfMaxY - psOptions->dfMinY <=
4565 1 : 2 * psOptions->dfYRes)
4566 0 : break;
4567 1 : psOptions->dfMaxY -= psOptions->dfYRes;
4568 : }
4569 2 : if (bBottomBlankLine)
4570 : {
4571 0 : if (psOptions->dfMaxY - psOptions->dfMinY <=
4572 0 : 2 * psOptions->dfYRes)
4573 0 : break;
4574 0 : psOptions->dfMinY += psOptions->dfYRes;
4575 : }
4576 2 : if (bLeftBlankCol)
4577 : {
4578 1 : if (psOptions->dfMaxX - psOptions->dfMinX <=
4579 1 : 2 * psOptions->dfXRes)
4580 0 : break;
4581 1 : psOptions->dfMinX += psOptions->dfXRes;
4582 : }
4583 2 : if (bRightBlankCol)
4584 : {
4585 1 : if (psOptions->dfMaxX - psOptions->dfMinX <=
4586 1 : 2 * psOptions->dfXRes)
4587 0 : break;
4588 1 : psOptions->dfMaxX -= psOptions->dfXRes;
4589 : }
4590 : }
4591 : }
4592 :
4593 38 : UpdateGeoTransformandAndPixelLines();
4594 : }
4595 :
4596 722 : else if (psOptions->nForcePixels != 0 && psOptions->nForceLines != 0)
4597 : {
4598 107 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4599 82 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4600 : {
4601 82 : psOptions->dfMinX = dfWrkMinX;
4602 82 : psOptions->dfMaxX = dfWrkMaxX;
4603 82 : psOptions->dfMaxY = dfWrkMaxY;
4604 82 : psOptions->dfMinY = dfWrkMinY;
4605 : }
4606 :
4607 107 : psOptions->dfXRes =
4608 107 : (psOptions->dfMaxX - psOptions->dfMinX) / psOptions->nForcePixels;
4609 107 : psOptions->dfYRes =
4610 107 : (psOptions->dfMaxY - psOptions->dfMinY) / psOptions->nForceLines;
4611 :
4612 107 : adfDstGeoTransform[0] = psOptions->dfMinX;
4613 107 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4614 107 : adfDstGeoTransform[1] = psOptions->dfXRes;
4615 107 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4616 :
4617 107 : nPixels = psOptions->nForcePixels;
4618 107 : nLines = psOptions->nForceLines;
4619 : }
4620 :
4621 615 : else if (psOptions->nForcePixels != 0)
4622 : {
4623 2 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4624 2 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4625 : {
4626 2 : psOptions->dfMinX = dfWrkMinX;
4627 2 : psOptions->dfMaxX = dfWrkMaxX;
4628 2 : psOptions->dfMaxY = dfWrkMaxY;
4629 2 : psOptions->dfMinY = dfWrkMinY;
4630 : }
4631 :
4632 2 : psOptions->dfXRes =
4633 2 : (psOptions->dfMaxX - psOptions->dfMinX) / psOptions->nForcePixels;
4634 2 : psOptions->dfYRes = psOptions->dfXRes;
4635 :
4636 2 : adfDstGeoTransform[0] = psOptions->dfMinX;
4637 2 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4638 2 : adfDstGeoTransform[1] = psOptions->dfXRes;
4639 4 : adfDstGeoTransform[5] = (psOptions->dfMaxY > psOptions->dfMinY)
4640 2 : ? -psOptions->dfYRes
4641 0 : : psOptions->dfYRes;
4642 :
4643 2 : nPixels = psOptions->nForcePixels;
4644 2 : nLines =
4645 2 : static_cast<int>((std::fabs(psOptions->dfMaxY - psOptions->dfMinY) +
4646 2 : (psOptions->dfYRes / 2.0)) /
4647 2 : psOptions->dfYRes);
4648 : }
4649 :
4650 613 : else if (psOptions->nForceLines != 0)
4651 : {
4652 2 : if (psOptions->dfMinX == 0.0 && psOptions->dfMinY == 0.0 &&
4653 2 : psOptions->dfMaxX == 0.0 && psOptions->dfMaxY == 0.0)
4654 : {
4655 2 : psOptions->dfMinX = dfWrkMinX;
4656 2 : psOptions->dfMaxX = dfWrkMaxX;
4657 2 : psOptions->dfMaxY = dfWrkMaxY;
4658 2 : psOptions->dfMinY = dfWrkMinY;
4659 : }
4660 :
4661 2 : psOptions->dfYRes =
4662 2 : (psOptions->dfMaxY - psOptions->dfMinY) / psOptions->nForceLines;
4663 2 : psOptions->dfXRes = std::fabs(psOptions->dfYRes);
4664 :
4665 2 : adfDstGeoTransform[0] = psOptions->dfMinX;
4666 2 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4667 2 : adfDstGeoTransform[1] = psOptions->dfXRes;
4668 2 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4669 :
4670 2 : nPixels = static_cast<int>((psOptions->dfMaxX - psOptions->dfMinX +
4671 2 : (psOptions->dfXRes / 2.0)) /
4672 2 : psOptions->dfXRes);
4673 2 : nLines = psOptions->nForceLines;
4674 : }
4675 :
4676 611 : else if (psOptions->dfMinX != 0.0 || psOptions->dfMinY != 0.0 ||
4677 536 : psOptions->dfMaxX != 0.0 || psOptions->dfMaxY != 0.0)
4678 : {
4679 75 : psOptions->dfXRes = adfDstGeoTransform[1];
4680 75 : psOptions->dfYRes = fabs(adfDstGeoTransform[5]);
4681 :
4682 75 : nPixels = static_cast<int>((psOptions->dfMaxX - psOptions->dfMinX +
4683 75 : (psOptions->dfXRes / 2.0)) /
4684 75 : psOptions->dfXRes);
4685 75 : nLines =
4686 75 : static_cast<int>((std::fabs(psOptions->dfMaxY - psOptions->dfMinY) +
4687 75 : (psOptions->dfYRes / 2.0)) /
4688 75 : psOptions->dfYRes);
4689 :
4690 75 : psOptions->dfXRes = (psOptions->dfMaxX - psOptions->dfMinX) / nPixels;
4691 75 : psOptions->dfYRes = (psOptions->dfMaxY - psOptions->dfMinY) / nLines;
4692 :
4693 75 : adfDstGeoTransform[0] = psOptions->dfMinX;
4694 75 : adfDstGeoTransform[3] = psOptions->dfMaxY;
4695 75 : adfDstGeoTransform[1] = psOptions->dfXRes;
4696 75 : adfDstGeoTransform[5] = -psOptions->dfYRes;
4697 : }
4698 :
4699 785 : if (EQUAL(pszFormat, "GTiff"))
4700 : {
4701 :
4702 : /* --------------------------------------------------------------------
4703 : */
4704 : /* Automatically set PHOTOMETRIC=RGB for GTiff when appropriate */
4705 : /* --------------------------------------------------------------------
4706 : */
4707 420 : if (apeColorInterpretations.size() >= 3 &&
4708 33 : apeColorInterpretations[0] == GCI_RedBand &&
4709 33 : apeColorInterpretations[1] == GCI_GreenBand &&
4710 486 : apeColorInterpretations[2] == GCI_BlueBand &&
4711 33 : aosCreateOptions.FetchNameValue("PHOTOMETRIC") == nullptr)
4712 : {
4713 33 : aosCreateOptions.SetNameValue("PHOTOMETRIC", "RGB");
4714 :
4715 : // Preserve potential ALPHA=PREMULTIPLIED from source alpha band
4716 66 : if (aosCreateOptions.FetchNameValue("ALPHA") == nullptr &&
4717 33 : apeColorInterpretations.size() == 4 &&
4718 71 : apeColorInterpretations[3] == GCI_AlphaBand &&
4719 5 : GDALGetRasterCount(pahSrcDS[0]) == 4)
4720 : {
4721 : const char *pszAlpha =
4722 5 : GDALGetMetadataItem(GDALGetRasterBand(pahSrcDS[0], 4),
4723 : "ALPHA", "IMAGE_STRUCTURE");
4724 5 : if (pszAlpha)
4725 : {
4726 1 : aosCreateOptions.SetNameValue("ALPHA", pszAlpha);
4727 : }
4728 : }
4729 : }
4730 :
4731 : /* The GTiff driver now supports writing band color interpretation */
4732 : /* in the TIFF_GDAL_METADATA tag */
4733 420 : bSetColorInterpretation = true;
4734 : }
4735 :
4736 : /* -------------------------------------------------------------------- */
4737 : /* Create the output file. */
4738 : /* -------------------------------------------------------------------- */
4739 785 : if (!psOptions->bQuiet)
4740 90 : printf("Creating output file that is %dP x %dL.\n", nPixels, nLines);
4741 :
4742 785 : hDstDS = GDALCreate(hDriver, pszFilename, nPixels, nLines, nDstBandCount,
4743 785 : eDT, aosCreateOptions.List());
4744 :
4745 785 : if (hDstDS == nullptr)
4746 : {
4747 1 : if (hCT != nullptr)
4748 1 : GDALDestroyColorTable(hCT);
4749 1 : return nullptr;
4750 : }
4751 :
4752 : /* -------------------------------------------------------------------- */
4753 : /* Write out the projection definition. */
4754 : /* -------------------------------------------------------------------- */
4755 784 : const char *pszDstMethod = CSLFetchNameValue(papszTO, "DST_METHOD");
4756 784 : if (pszDstMethod == nullptr || !EQUAL(pszDstMethod, "NO_GEOTRANSFORM"))
4757 : {
4758 782 : OGRSpatialReference oTargetSRS;
4759 782 : oTargetSRS.SetFromUserInput(osThisTargetSRS);
4760 782 : oTargetSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
4761 :
4762 782 : if (oTargetSRS.IsDynamic())
4763 : {
4764 524 : double dfCoordEpoch = CPLAtof(CSLFetchNameValueDef(
4765 : papszTO, "DST_COORDINATE_EPOCH",
4766 : CSLFetchNameValueDef(papszTO, "COORDINATE_EPOCH", "0")));
4767 524 : if (dfCoordEpoch == 0)
4768 : {
4769 : const OGRSpatialReferenceH hSrcSRS =
4770 524 : GDALGetSpatialRef(pahSrcDS[0]);
4771 524 : const char *pszMethod = CSLFetchNameValue(papszTO, "METHOD");
4772 524 : if (hSrcSRS &&
4773 3 : (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")))
4774 : {
4775 486 : dfCoordEpoch = OSRGetCoordinateEpoch(hSrcSRS);
4776 : }
4777 : }
4778 524 : if (dfCoordEpoch > 0)
4779 1 : oTargetSRS.SetCoordinateEpoch(dfCoordEpoch);
4780 : }
4781 :
4782 782 : if (GDALSetSpatialRef(hDstDS, OGRSpatialReference::ToHandle(
4783 1564 : &oTargetSRS)) == CE_Failure ||
4784 782 : GDALSetGeoTransform(hDstDS, adfDstGeoTransform) == CE_Failure)
4785 : {
4786 0 : if (hCT != nullptr)
4787 0 : GDALDestroyColorTable(hCT);
4788 0 : GDALClose(hDstDS);
4789 0 : return nullptr;
4790 782 : }
4791 : }
4792 : else
4793 : {
4794 2 : adfDstGeoTransform[3] += adfDstGeoTransform[5] * nLines;
4795 2 : adfDstGeoTransform[5] = fabs(adfDstGeoTransform[5]);
4796 : }
4797 :
4798 784 : if (phTransformArg && *phTransformArg != nullptr)
4799 759 : GDALSetGenImgProjTransformerDstGeoTransform(*phTransformArg,
4800 : adfDstGeoTransform);
4801 :
4802 : /* -------------------------------------------------------------------- */
4803 : /* Try to set color interpretation of source bands to target */
4804 : /* dataset. */
4805 : /* FIXME? We should likely do that for other drivers than VRT & */
4806 : /* GTiff but it might create spurious .aux.xml files (at least */
4807 : /* with HFA, and netCDF) */
4808 : /* -------------------------------------------------------------------- */
4809 784 : if (bVRT || bSetColorInterpretation)
4810 : {
4811 508 : int nBandsToCopy = static_cast<int>(apeColorInterpretations.size());
4812 508 : if (psOptions->bEnableSrcAlpha)
4813 10 : nBandsToCopy--;
4814 1140 : for (int iBand = 0; iBand < nBandsToCopy; iBand++)
4815 : {
4816 632 : GDALSetRasterColorInterpretation(
4817 : GDALGetRasterBand(hDstDS, iBand + 1),
4818 632 : apeColorInterpretations[iBand]);
4819 : }
4820 : }
4821 :
4822 : /* -------------------------------------------------------------------- */
4823 : /* Try to set color interpretation of output file alpha band. */
4824 : /* -------------------------------------------------------------------- */
4825 784 : if (psOptions->bEnableDstAlpha)
4826 : {
4827 58 : GDALSetRasterColorInterpretation(
4828 : GDALGetRasterBand(hDstDS, nDstBandCount), GCI_AlphaBand);
4829 : }
4830 :
4831 : /* -------------------------------------------------------------------- */
4832 : /* Copy the raster attribute table, if required. */
4833 : /* -------------------------------------------------------------------- */
4834 784 : if (hRAT != nullptr)
4835 : {
4836 0 : GDALSetDefaultRAT(GDALGetRasterBand(hDstDS, 1), hRAT);
4837 : }
4838 :
4839 : /* -------------------------------------------------------------------- */
4840 : /* Copy the color table, if required. */
4841 : /* -------------------------------------------------------------------- */
4842 784 : if (hCT != nullptr)
4843 : {
4844 2 : GDALSetRasterColorTable(GDALGetRasterBand(hDstDS, 1), hCT);
4845 2 : GDALDestroyColorTable(hCT);
4846 : }
4847 :
4848 : /* -------------------------------------------------------------------- */
4849 : /* Copy scale/offset if found on source */
4850 : /* -------------------------------------------------------------------- */
4851 784 : if (nSrcCount == 1)
4852 : {
4853 760 : GDALDataset *poSrcDS = GDALDataset::FromHandle(pahSrcDS[0]);
4854 760 : GDALDataset *poDstDS = GDALDataset::FromHandle(hDstDS);
4855 :
4856 760 : int nBandsToCopy = nDstBandCount;
4857 760 : if (psOptions->bEnableDstAlpha)
4858 55 : nBandsToCopy--;
4859 760 : nBandsToCopy = std::min(nBandsToCopy, poSrcDS->GetRasterCount());
4860 :
4861 1741 : for (int i = 0; i < nBandsToCopy; i++)
4862 : {
4863 1007 : auto poSrcBand = poSrcDS->GetRasterBand(
4864 981 : psOptions->anSrcBands.empty() ? i + 1
4865 26 : : psOptions->anSrcBands[i]);
4866 1007 : auto poDstBand = poDstDS->GetRasterBand(
4867 981 : psOptions->anDstBands.empty() ? i + 1
4868 26 : : psOptions->anDstBands[i]);
4869 981 : if (poSrcBand && poDstBand)
4870 : {
4871 979 : int bHasScale = FALSE;
4872 979 : const double dfScale = poSrcBand->GetScale(&bHasScale);
4873 979 : if (bHasScale)
4874 20 : poDstBand->SetScale(dfScale);
4875 :
4876 979 : int bHasOffset = FALSE;
4877 979 : const double dfOffset = poSrcBand->GetOffset(&bHasOffset);
4878 979 : if (bHasOffset)
4879 20 : poDstBand->SetOffset(dfOffset);
4880 : }
4881 : }
4882 : }
4883 :
4884 784 : return hDstDS;
4885 : }
4886 :
4887 : /************************************************************************/
4888 : /* GeoTransform_Transformer() */
4889 : /* */
4890 : /* Convert points from georef coordinates to pixel/line based */
4891 : /* on a geotransform. */
4892 : /************************************************************************/
4893 :
4894 : class CutlineTransformer : public OGRCoordinateTransformation
4895 : {
4896 : CPL_DISALLOW_COPY_ASSIGN(CutlineTransformer)
4897 :
4898 : public:
4899 : void *hSrcImageTransformer = nullptr;
4900 :
4901 33 : explicit CutlineTransformer(void *hTransformArg)
4902 33 : : hSrcImageTransformer(hTransformArg)
4903 : {
4904 33 : }
4905 :
4906 0 : virtual const OGRSpatialReference *GetSourceCS() const override
4907 : {
4908 0 : return nullptr;
4909 : }
4910 :
4911 135 : virtual const OGRSpatialReference *GetTargetCS() const override
4912 : {
4913 135 : return nullptr;
4914 : }
4915 :
4916 33 : virtual ~CutlineTransformer()
4917 33 : {
4918 33 : GDALDestroyTransformer(hSrcImageTransformer);
4919 33 : }
4920 :
4921 47 : virtual int Transform(size_t nCount, double *x, double *y, double *z,
4922 : double * /* t */, int *pabSuccess) override
4923 : {
4924 47 : CPLAssert(nCount <=
4925 : static_cast<size_t>(std::numeric_limits<int>::max()));
4926 47 : return GDALGenImgProjTransform(hSrcImageTransformer, TRUE,
4927 : static_cast<int>(nCount), x, y, z,
4928 47 : pabSuccess);
4929 : }
4930 :
4931 0 : virtual OGRCoordinateTransformation *Clone() const override
4932 : {
4933 : return new CutlineTransformer(
4934 0 : GDALCloneTransformer(hSrcImageTransformer));
4935 : }
4936 :
4937 0 : virtual OGRCoordinateTransformation *GetInverse() const override
4938 : {
4939 0 : return nullptr;
4940 : }
4941 : };
4942 :
4943 236 : static double GetMaximumSegmentLength(OGRGeometry *poGeom)
4944 : {
4945 236 : switch (wkbFlatten(poGeom->getGeometryType()))
4946 : {
4947 82 : case wkbLineString:
4948 : {
4949 82 : OGRLineString *poLS = static_cast<OGRLineString *>(poGeom);
4950 82 : double dfMaxSquaredLength = 0.0;
4951 12824 : for (int i = 0; i < poLS->getNumPoints() - 1; i++)
4952 : {
4953 12742 : double dfDeltaX = poLS->getX(i + 1) - poLS->getX(i);
4954 12742 : double dfDeltaY = poLS->getY(i + 1) - poLS->getY(i);
4955 12742 : double dfSquaredLength =
4956 12742 : dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY;
4957 12742 : dfMaxSquaredLength =
4958 12742 : std::max(dfMaxSquaredLength, dfSquaredLength);
4959 : }
4960 82 : return sqrt(dfMaxSquaredLength);
4961 : }
4962 :
4963 80 : case wkbPolygon:
4964 : {
4965 80 : OGRPolygon *poPoly = static_cast<OGRPolygon *>(poGeom);
4966 : double dfMaxLength =
4967 80 : GetMaximumSegmentLength(poPoly->getExteriorRing());
4968 82 : for (int i = 0; i < poPoly->getNumInteriorRings(); i++)
4969 : {
4970 2 : dfMaxLength = std::max(
4971 : dfMaxLength,
4972 2 : GetMaximumSegmentLength(poPoly->getInteriorRing(i)));
4973 : }
4974 80 : return dfMaxLength;
4975 : }
4976 :
4977 74 : case wkbMultiPolygon:
4978 : {
4979 74 : OGRMultiPolygon *poMP = static_cast<OGRMultiPolygon *>(poGeom);
4980 74 : double dfMaxLength = 0.0;
4981 152 : for (int i = 0; i < poMP->getNumGeometries(); i++)
4982 : {
4983 78 : dfMaxLength =
4984 78 : std::max(dfMaxLength,
4985 78 : GetMaximumSegmentLength(poMP->getGeometryRef(i)));
4986 : }
4987 74 : return dfMaxLength;
4988 : }
4989 :
4990 0 : default:
4991 0 : CPLAssert(false);
4992 : return 0.0;
4993 : }
4994 : }
4995 :
4996 : /************************************************************************/
4997 : /* RemoveZeroWidthSlivers() */
4998 : /* */
4999 : /* Such slivers can cause issues after reprojection. */
5000 : /************************************************************************/
5001 :
5002 104 : static void RemoveZeroWidthSlivers(OGRGeometry *poGeom)
5003 : {
5004 104 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
5005 104 : if (eType == wkbMultiPolygon)
5006 : {
5007 32 : auto poMP = poGeom->toMultiPolygon();
5008 32 : int nNumGeometries = poMP->getNumGeometries();
5009 66 : for (int i = 0; i < nNumGeometries; /* incremented in loop */)
5010 : {
5011 34 : auto poPoly = poMP->getGeometryRef(i);
5012 34 : RemoveZeroWidthSlivers(poPoly);
5013 34 : if (poPoly->IsEmpty())
5014 : {
5015 1 : CPLDebug("WARP",
5016 : "RemoveZeroWidthSlivers: removing empty polygon");
5017 1 : poMP->removeGeometry(i, /* bDelete = */ true);
5018 1 : --nNumGeometries;
5019 : }
5020 : else
5021 : {
5022 33 : ++i;
5023 : }
5024 : }
5025 : }
5026 72 : else if (eType == wkbPolygon)
5027 : {
5028 35 : auto poPoly = poGeom->toPolygon();
5029 35 : if (auto poExteriorRing = poPoly->getExteriorRing())
5030 : {
5031 35 : RemoveZeroWidthSlivers(poExteriorRing);
5032 35 : if (poExteriorRing->getNumPoints() < 4)
5033 : {
5034 1 : poPoly->empty();
5035 1 : return;
5036 : }
5037 : }
5038 34 : int nNumInteriorRings = poPoly->getNumInteriorRings();
5039 36 : for (int i = 0; i < nNumInteriorRings; /* incremented in loop */)
5040 : {
5041 2 : auto poRing = poPoly->getInteriorRing(i);
5042 2 : RemoveZeroWidthSlivers(poRing);
5043 2 : if (poRing->getNumPoints() < 4)
5044 : {
5045 1 : CPLDebug(
5046 : "WARP",
5047 : "RemoveZeroWidthSlivers: removing empty interior ring");
5048 1 : constexpr int OFFSET_EXTERIOR_RING = 1;
5049 1 : poPoly->removeRing(i + OFFSET_EXTERIOR_RING,
5050 : /* bDelete = */ true);
5051 1 : --nNumInteriorRings;
5052 : }
5053 : else
5054 : {
5055 1 : ++i;
5056 : }
5057 : }
5058 : }
5059 37 : else if (eType == wkbLineString)
5060 : {
5061 37 : OGRLineString *poLS = poGeom->toLineString();
5062 37 : int numPoints = poLS->getNumPoints();
5063 182 : for (int i = 1; i < numPoints - 1;)
5064 : {
5065 145 : const double x1 = poLS->getX(i - 1);
5066 145 : const double y1 = poLS->getY(i - 1);
5067 145 : const double x2 = poLS->getX(i);
5068 145 : const double y2 = poLS->getY(i);
5069 145 : const double x3 = poLS->getX(i + 1);
5070 145 : const double y3 = poLS->getY(i + 1);
5071 145 : const double dx1 = x2 - x1;
5072 145 : const double dy1 = y2 - y1;
5073 145 : const double dx2 = x3 - x2;
5074 145 : const double dy2 = y3 - y2;
5075 145 : const double scalar_product = dx1 * dx2 + dy1 * dy2;
5076 145 : const double square_scalar_product =
5077 : scalar_product * scalar_product;
5078 145 : const double square_norm1 = dx1 * dx1 + dy1 * dy1;
5079 145 : const double square_norm2 = dx2 * dx2 + dy2 * dy2;
5080 145 : const double square_norm1_mult_norm2 = square_norm1 * square_norm2;
5081 145 : if (scalar_product < 0 &&
5082 20 : fabs(square_scalar_product - square_norm1_mult_norm2) <=
5083 20 : 1e-15 * square_norm1_mult_norm2)
5084 : {
5085 5 : CPLDebug("WARP",
5086 : "RemoveZeroWidthSlivers: removing point %.10g %.10g",
5087 : x2, y2);
5088 5 : poLS->removePoint(i);
5089 5 : numPoints--;
5090 : }
5091 : else
5092 : {
5093 140 : ++i;
5094 : }
5095 : }
5096 : }
5097 : }
5098 :
5099 : /************************************************************************/
5100 : /* TransformCutlineToSource() */
5101 : /* */
5102 : /* Transform cutline from its SRS to source pixel/line coordinates.*/
5103 : /************************************************************************/
5104 33 : static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,
5105 : OGRGeometry *poCutline,
5106 : char ***ppapszWarpOptions,
5107 : CSLConstList papszTO_In)
5108 :
5109 : {
5110 33 : RemoveZeroWidthSlivers(poCutline);
5111 :
5112 66 : auto poMultiPolygon = std::unique_ptr<OGRGeometry>(poCutline->clone());
5113 :
5114 : /* -------------------------------------------------------------------- */
5115 : /* Checkout that if there's a cutline SRS, there's also a raster */
5116 : /* one. */
5117 : /* -------------------------------------------------------------------- */
5118 33 : std::unique_ptr<OGRSpatialReference> poRasterSRS;
5119 : const CPLString osProjection =
5120 66 : GetSrcDSProjection(GDALDataset::ToHandle(poSrcDS), papszTO_In);
5121 33 : if (!osProjection.empty())
5122 : {
5123 31 : poRasterSRS = std::make_unique<OGRSpatialReference>();
5124 31 : poRasterSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
5125 31 : if (poRasterSRS->SetFromUserInput(osProjection) != OGRERR_NONE)
5126 : {
5127 0 : poRasterSRS.reset();
5128 : }
5129 : }
5130 :
5131 33 : std::unique_ptr<OGRSpatialReference> poDstSRS;
5132 33 : const char *pszThisTargetSRS = CSLFetchNameValue(papszTO_In, "DST_SRS");
5133 33 : if (pszThisTargetSRS)
5134 : {
5135 8 : poDstSRS = std::make_unique<OGRSpatialReference>();
5136 8 : poDstSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
5137 8 : if (poDstSRS->SetFromUserInput(pszThisTargetSRS) != OGRERR_NONE)
5138 : {
5139 0 : return CE_Failure;
5140 : }
5141 : }
5142 25 : else if (poRasterSRS)
5143 : {
5144 23 : poDstSRS.reset(poRasterSRS->Clone());
5145 : }
5146 :
5147 : /* -------------------------------------------------------------------- */
5148 : /* Extract the cutline SRS. */
5149 : /* -------------------------------------------------------------------- */
5150 : const OGRSpatialReference *poCutlineSRS =
5151 33 : poMultiPolygon->getSpatialReference();
5152 :
5153 : /* -------------------------------------------------------------------- */
5154 : /* Detect if there's no transform at all involved, in which case */
5155 : /* we can avoid densification. */
5156 : /* -------------------------------------------------------------------- */
5157 33 : bool bMayNeedDensify = true;
5158 64 : if (poRasterSRS && poCutlineSRS && poRasterSRS->IsSame(poCutlineSRS) &&
5159 26 : poSrcDS->GetGCPCount() == 0 && !poSrcDS->GetMetadata("RPC") &&
5160 22 : !poSrcDS->GetMetadata("GEOLOCATION") &&
5161 86 : !CSLFetchNameValue(papszTO_In, "GEOLOC_ARRAY") &&
5162 22 : !CSLFetchNameValue(papszTO_In, "SRC_GEOLOC_ARRAY"))
5163 : {
5164 44 : CPLStringList aosTOTmp(papszTO_In);
5165 22 : aosTOTmp.SetNameValue("SRC_SRS", nullptr);
5166 22 : aosTOTmp.SetNameValue("DST_SRS", nullptr);
5167 22 : if (aosTOTmp.size() == 0)
5168 : {
5169 22 : bMayNeedDensify = false;
5170 : }
5171 : }
5172 :
5173 : /* -------------------------------------------------------------------- */
5174 : /* Compare source raster SRS and cutline SRS */
5175 : /* -------------------------------------------------------------------- */
5176 33 : if (poRasterSRS && poCutlineSRS)
5177 : {
5178 : /* OK, we will reproject */
5179 : }
5180 2 : else if (poRasterSRS && !poCutlineSRS)
5181 : {
5182 0 : CPLError(
5183 : CE_Warning, CPLE_AppDefined,
5184 : "the source raster dataset has a SRS, but the cutline features\n"
5185 : "not. We assume that the cutline coordinates are expressed in the "
5186 : "destination SRS.\n"
5187 : "If not, cutline results may be incorrect.");
5188 : }
5189 2 : else if (!poRasterSRS && poCutlineSRS)
5190 : {
5191 0 : CPLError(CE_Warning, CPLE_AppDefined,
5192 : "the input vector layer has a SRS, but the source raster "
5193 : "dataset does not.\n"
5194 : "Cutline results may be incorrect.");
5195 : }
5196 :
5197 : auto poCTCutlineToSrc = CreateCTCutlineToSrc(
5198 66 : poRasterSRS.get(), poDstSRS.get(), poCutlineSRS, papszTO_In);
5199 :
5200 66 : CPLStringList aosTO(papszTO_In);
5201 :
5202 33 : if (pszThisTargetSRS && !osProjection.empty())
5203 : {
5204 : // Avoid any reprojection when using the GenImgProjTransformer
5205 8 : aosTO.SetNameValue("DST_SRS", osProjection.c_str());
5206 : }
5207 33 : aosTO.SetNameValue("SRC_COORDINATE_EPOCH", nullptr);
5208 33 : aosTO.SetNameValue("DST_COORDINATE_EPOCH", nullptr);
5209 33 : aosTO.SetNameValue("COORDINATE_OPERATION", nullptr);
5210 :
5211 : /* -------------------------------------------------------------------- */
5212 : /* It may be unwise to let the mask geometry be re-wrapped by */
5213 : /* the CENTER_LONG machinery as this can easily screw up world */
5214 : /* spanning masks and invert the mask topology. */
5215 : /* -------------------------------------------------------------------- */
5216 33 : aosTO.SetNameValue("INSERT_CENTER_LONG", "FALSE");
5217 :
5218 : /* -------------------------------------------------------------------- */
5219 : /* Transform the geometry to pixel/line coordinates. */
5220 : /* -------------------------------------------------------------------- */
5221 : /* The cutline transformer will *invert* the hSrcImageTransformer */
5222 : /* so it will convert from the source SRS to the source pixel/line */
5223 : /* coordinates */
5224 : CutlineTransformer oTransformer(GDALCreateGenImgProjTransformer2(
5225 66 : GDALDataset::ToHandle(poSrcDS), nullptr, aosTO.List()));
5226 :
5227 33 : if (oTransformer.hSrcImageTransformer == nullptr)
5228 : {
5229 0 : return CE_Failure;
5230 : }
5231 :
5232 : // Some transforms like RPC can transform a valid geometry into an invalid
5233 : // one if the node density of the input geometry isn't sufficient before
5234 : // reprojection. So after an initial reprojection, we check that the
5235 : // maximum length of a segment is no longer than 1 pixel, and if not,
5236 : // we densify the input geometry before doing a new reprojection
5237 : const double dfMaxLengthInSpatUnits =
5238 33 : GetMaximumSegmentLength(poMultiPolygon.get());
5239 33 : OGRErr eErr = OGRERR_NONE;
5240 33 : if (poCTCutlineToSrc)
5241 : {
5242 10 : poMultiPolygon.reset(OGRGeometryFactory::transformWithOptions(
5243 5 : poMultiPolygon.get(), poCTCutlineToSrc.get(), nullptr));
5244 5 : if (!poMultiPolygon)
5245 : {
5246 0 : eErr = OGRERR_FAILURE;
5247 0 : poMultiPolygon.reset(poCutline->clone());
5248 0 : poMultiPolygon->transform(poCTCutlineToSrc.get());
5249 : }
5250 : }
5251 33 : if (poMultiPolygon->transform(&oTransformer) != OGRERR_NONE)
5252 : {
5253 0 : CPLError(CE_Failure, CPLE_AppDefined,
5254 : "poMultiPolygon->transform(&oTransformer) failed at line %d",
5255 : __LINE__);
5256 0 : eErr = OGRERR_FAILURE;
5257 : }
5258 : const double dfInitialMaxLengthInPixels =
5259 33 : GetMaximumSegmentLength(poMultiPolygon.get());
5260 :
5261 33 : CPLPushErrorHandler(CPLQuietErrorHandler);
5262 : const bool bWasValidInitially =
5263 33 : ValidateCutline(poMultiPolygon.get(), false);
5264 33 : CPLPopErrorHandler();
5265 33 : if (!bWasValidInitially)
5266 : {
5267 3 : CPLDebug("WARP", "Cutline is not valid after initial reprojection");
5268 3 : char *pszWKT = nullptr;
5269 3 : poMultiPolygon->exportToWkt(&pszWKT);
5270 3 : CPLDebug("GDALWARP", "WKT = \"%s\"", pszWKT ? pszWKT : "(null)");
5271 3 : CPLFree(pszWKT);
5272 : }
5273 :
5274 33 : bool bDensify = false;
5275 33 : if (bMayNeedDensify && eErr == OGRERR_NONE &&
5276 : dfInitialMaxLengthInPixels > 1.0)
5277 : {
5278 : const char *pszDensifyCutline =
5279 10 : CPLGetConfigOption("GDALWARP_DENSIFY_CUTLINE", "YES");
5280 10 : if (EQUAL(pszDensifyCutline, "ONLY_IF_INVALID"))
5281 : {
5282 1 : bDensify = (OGRGeometryFactory::haveGEOS() && !bWasValidInitially);
5283 : }
5284 9 : else if (CSLFetchNameValue(*ppapszWarpOptions, "CUTLINE_BLEND_DIST") !=
5285 9 : nullptr &&
5286 0 : CPLGetConfigOption("GDALWARP_DENSIFY_CUTLINE", nullptr) ==
5287 : nullptr)
5288 : {
5289 : // TODO: we should only emit this message if a
5290 : // transform/reprojection will be actually done
5291 0 : CPLDebug("WARP",
5292 : "Densification of cutline could perhaps be useful but as "
5293 : "CUTLINE_BLEND_DIST is used, this could be very slow. So "
5294 : "disabled "
5295 : "unless GDALWARP_DENSIFY_CUTLINE=YES is explicitly "
5296 : "specified as configuration option");
5297 : }
5298 : else
5299 : {
5300 9 : bDensify = CPLTestBool(pszDensifyCutline);
5301 : }
5302 : }
5303 33 : if (bDensify)
5304 : {
5305 9 : CPLDebug("WARP",
5306 : "Cutline maximum segment size was %.0f pixel after "
5307 : "reprojection to source coordinates.",
5308 : dfInitialMaxLengthInPixels);
5309 :
5310 : // Densify and reproject with the aim of having a 1 pixel density
5311 9 : double dfSegmentSize =
5312 : dfMaxLengthInSpatUnits / dfInitialMaxLengthInPixels;
5313 9 : const int MAX_ITERATIONS = 10;
5314 10 : for (int i = 0; i < MAX_ITERATIONS; i++)
5315 : {
5316 10 : poMultiPolygon.reset(poCutline->clone());
5317 10 : poMultiPolygon->segmentize(dfSegmentSize);
5318 10 : if (i == MAX_ITERATIONS - 1)
5319 : {
5320 0 : char *pszWKT = nullptr;
5321 0 : poMultiPolygon->exportToWkt(&pszWKT);
5322 0 : CPLDebug("WARP",
5323 : "WKT of polygon after densification with segment size "
5324 : "= %f: %s",
5325 : dfSegmentSize, pszWKT);
5326 0 : CPLFree(pszWKT);
5327 : }
5328 10 : eErr = OGRERR_NONE;
5329 10 : if (poCTCutlineToSrc)
5330 : {
5331 10 : poMultiPolygon.reset(OGRGeometryFactory::transformWithOptions(
5332 5 : poMultiPolygon.get(), poCTCutlineToSrc.get(), nullptr));
5333 5 : if (!poMultiPolygon)
5334 : {
5335 0 : eErr = OGRERR_FAILURE;
5336 0 : break;
5337 : }
5338 : }
5339 10 : if (poMultiPolygon->transform(&oTransformer) != OGRERR_NONE)
5340 0 : eErr = OGRERR_FAILURE;
5341 10 : if (eErr == OGRERR_NONE)
5342 : {
5343 : const double dfMaxLengthInPixels =
5344 10 : GetMaximumSegmentLength(poMultiPolygon.get());
5345 10 : if (bWasValidInitially)
5346 : {
5347 : // In some cases, the densification itself results in a
5348 : // reprojected invalid polygon due to the non-linearity of
5349 : // RPC DEM transformation, so in those cases, try a less
5350 : // dense cutline
5351 8 : CPLPushErrorHandler(CPLQuietErrorHandler);
5352 : const bool bIsValid =
5353 8 : ValidateCutline(poMultiPolygon.get(), false);
5354 8 : CPLPopErrorHandler();
5355 8 : if (!bIsValid)
5356 : {
5357 1 : if (i == MAX_ITERATIONS - 1)
5358 : {
5359 0 : char *pszWKT = nullptr;
5360 0 : poMultiPolygon->exportToWkt(&pszWKT);
5361 0 : CPLDebug("WARP",
5362 : "After densification, cutline maximum "
5363 : "segment size is now %.0f pixel, "
5364 : "but cutline is invalid. %s",
5365 : dfMaxLengthInPixels, pszWKT);
5366 0 : CPLFree(pszWKT);
5367 0 : break;
5368 : }
5369 1 : CPLDebug("WARP",
5370 : "After densification, cutline maximum segment "
5371 : "size is now %.0f pixel, "
5372 : "but cutline is invalid. So trying a less "
5373 : "dense cutline.",
5374 : dfMaxLengthInPixels);
5375 1 : dfSegmentSize *= 2;
5376 1 : continue;
5377 : }
5378 : }
5379 9 : CPLDebug("WARP",
5380 : "After densification, cutline maximum segment size is "
5381 : "now %.0f pixel.",
5382 : dfMaxLengthInPixels);
5383 : }
5384 9 : break;
5385 : }
5386 : }
5387 :
5388 33 : if (eErr == OGRERR_FAILURE)
5389 : {
5390 0 : if (CPLTestBool(
5391 : CPLGetConfigOption("GDALWARP_IGNORE_BAD_CUTLINE", "NO")))
5392 0 : CPLError(CE_Warning, CPLE_AppDefined,
5393 : "Cutline transformation failed");
5394 : else
5395 : {
5396 0 : CPLError(CE_Failure, CPLE_AppDefined,
5397 : "Cutline transformation failed");
5398 0 : return CE_Failure;
5399 : }
5400 : }
5401 33 : else if (!ValidateCutline(poMultiPolygon.get(), true))
5402 : {
5403 1 : return CE_Failure;
5404 : }
5405 :
5406 : // Optimization: if the cutline contains the footprint of the source
5407 : // dataset, no need to use the cutline.
5408 32 : if (OGRGeometryFactory::haveGEOS()
5409 : #ifdef DEBUG
5410 : // Env var just for debugging purposes
5411 32 : && !CPLTestBool(CPLGetConfigOption(
5412 : "GDALWARP_SKIP_CUTLINE_CONTAINMENT_TEST", "NO"))
5413 : #endif
5414 : )
5415 : {
5416 31 : const double dfCutlineBlendDist = CPLAtof(CSLFetchNameValueDef(
5417 : *ppapszWarpOptions, "CUTLINE_BLEND_DIST", "0"));
5418 31 : auto poRing = std::make_unique<OGRLinearRing>();
5419 31 : poRing->addPoint(-dfCutlineBlendDist, -dfCutlineBlendDist);
5420 31 : poRing->addPoint(-dfCutlineBlendDist,
5421 31 : dfCutlineBlendDist + poSrcDS->GetRasterYSize());
5422 62 : poRing->addPoint(dfCutlineBlendDist + poSrcDS->GetRasterXSize(),
5423 31 : dfCutlineBlendDist + poSrcDS->GetRasterYSize());
5424 31 : poRing->addPoint(dfCutlineBlendDist + poSrcDS->GetRasterXSize(),
5425 : -dfCutlineBlendDist);
5426 31 : poRing->addPoint(-dfCutlineBlendDist, -dfCutlineBlendDist);
5427 31 : OGRPolygon oSrcDSFootprint;
5428 31 : oSrcDSFootprint.addRing(std::move(poRing));
5429 31 : OGREnvelope sSrcDSEnvelope;
5430 31 : oSrcDSFootprint.getEnvelope(&sSrcDSEnvelope);
5431 31 : OGREnvelope sCutlineEnvelope;
5432 31 : poMultiPolygon->getEnvelope(&sCutlineEnvelope);
5433 32 : if (sCutlineEnvelope.Contains(sSrcDSEnvelope) &&
5434 1 : poMultiPolygon->Contains(&oSrcDSFootprint))
5435 : {
5436 1 : CPLDebug("WARP", "Source dataset fully contained within cutline.");
5437 1 : return CE_None;
5438 : }
5439 : }
5440 :
5441 : /* -------------------------------------------------------------------- */
5442 : /* Convert aggregate geometry into WKT. */
5443 : /* -------------------------------------------------------------------- */
5444 31 : char *pszWKT = nullptr;
5445 31 : poMultiPolygon->exportToWkt(&pszWKT);
5446 : // fprintf(stderr, "WKT = \"%s\"\n", pszWKT ? pszWKT : "(null)");
5447 :
5448 31 : *ppapszWarpOptions = CSLSetNameValue(*ppapszWarpOptions, "CUTLINE", pszWKT);
5449 31 : CPLFree(pszWKT);
5450 31 : return CE_None;
5451 : }
5452 :
5453 43 : static void RemoveConflictingMetadata(GDALMajorObjectH hObj,
5454 : CSLConstList papszSrcMetadata,
5455 : const char *pszValueConflict)
5456 : {
5457 43 : if (hObj == nullptr)
5458 0 : return;
5459 :
5460 38 : for (const auto &[pszKey, pszValue] :
5461 81 : cpl::IterateNameValue(papszSrcMetadata))
5462 : {
5463 19 : const char *pszValueComp = GDALGetMetadataItem(hObj, pszKey, nullptr);
5464 19 : if (pszValueComp == nullptr || (!EQUAL(pszValue, pszValueComp) &&
5465 2 : !EQUAL(pszValueComp, pszValueConflict)))
5466 : {
5467 3 : if (STARTS_WITH(pszKey, "STATISTICS_"))
5468 1 : GDALSetMetadataItem(hObj, pszKey, nullptr, nullptr);
5469 : else
5470 2 : GDALSetMetadataItem(hObj, pszKey, pszValueConflict, nullptr);
5471 : }
5472 : }
5473 : }
5474 :
5475 : /************************************************************************/
5476 : /* IsValidSRS */
5477 : /************************************************************************/
5478 :
5479 222 : static bool IsValidSRS(const char *pszUserInput)
5480 :
5481 : {
5482 : OGRSpatialReferenceH hSRS;
5483 222 : bool bRes = true;
5484 :
5485 222 : hSRS = OSRNewSpatialReference(nullptr);
5486 222 : if (OSRSetFromUserInput(hSRS, pszUserInput) != OGRERR_NONE)
5487 : {
5488 0 : bRes = false;
5489 : }
5490 :
5491 222 : OSRDestroySpatialReference(hSRS);
5492 :
5493 222 : return bRes;
5494 : }
5495 :
5496 : /************************************************************************/
5497 : /* GDALWarpAppOptionsGetParser() */
5498 : /************************************************************************/
5499 :
5500 : static std::unique_ptr<GDALArgumentParser>
5501 905 : GDALWarpAppOptionsGetParser(GDALWarpAppOptions *psOptions,
5502 : GDALWarpAppOptionsForBinary *psOptionsForBinary)
5503 : {
5504 : auto argParser = std::make_unique<GDALArgumentParser>(
5505 905 : "gdalwarp", /* bForBinary=*/psOptionsForBinary != nullptr);
5506 :
5507 905 : argParser->add_description(_("Image reprojection and warping utility."));
5508 :
5509 905 : argParser->add_epilog(
5510 905 : _("For more details, consult https://gdal.org/programs/gdalwarp.html"));
5511 :
5512 : argParser->add_quiet_argument(
5513 905 : psOptionsForBinary ? &psOptionsForBinary->bQuiet : nullptr);
5514 :
5515 905 : argParser->add_argument("-overwrite")
5516 905 : .flag()
5517 : .action(
5518 113 : [psOptionsForBinary](const std::string &)
5519 : {
5520 60 : if (psOptionsForBinary)
5521 53 : psOptionsForBinary->bOverwrite = true;
5522 905 : })
5523 905 : .help(_("Overwrite the target dataset if it already exists."));
5524 :
5525 905 : argParser->add_output_format_argument(psOptions->osFormat);
5526 :
5527 905 : argParser->add_argument("-co")
5528 1810 : .metavar("<NAME>=<VALUE>")
5529 905 : .append()
5530 : .action(
5531 235 : [psOptions, psOptionsForBinary](const std::string &s)
5532 : {
5533 114 : psOptions->aosCreateOptions.AddString(s.c_str());
5534 114 : psOptions->bCreateOutput = true;
5535 :
5536 114 : if (psOptionsForBinary)
5537 7 : psOptionsForBinary->aosCreateOptions.AddString(s.c_str());
5538 905 : })
5539 905 : .help(_("Creation option(s)."));
5540 :
5541 905 : argParser->add_argument("-s_srs")
5542 1810 : .metavar("<srs_def>")
5543 : .action(
5544 76 : [psOptions](const std::string &s)
5545 : {
5546 38 : if (!IsValidSRS(s.c_str()))
5547 : {
5548 0 : throw std::invalid_argument("Invalid SRS for -s_srs");
5549 : }
5550 : psOptions->aosTransformerOptions.SetNameValue("SRC_SRS",
5551 38 : s.c_str());
5552 943 : })
5553 905 : .help(_("Set source spatial reference."));
5554 :
5555 905 : argParser->add_argument("-t_srs")
5556 1810 : .metavar("<srs_def>")
5557 : .action(
5558 358 : [psOptions](const std::string &s)
5559 : {
5560 179 : if (!IsValidSRS(s.c_str()))
5561 : {
5562 0 : throw std::invalid_argument("Invalid SRS for -t_srs");
5563 : }
5564 : psOptions->aosTransformerOptions.SetNameValue("DST_SRS",
5565 179 : s.c_str());
5566 1084 : })
5567 905 : .help(_("Set target spatial reference."));
5568 :
5569 : {
5570 905 : auto &group = argParser->add_mutually_exclusive_group();
5571 905 : group.add_argument("-srcalpha")
5572 905 : .flag()
5573 905 : .store_into(psOptions->bEnableSrcAlpha)
5574 : .help(_("Force the last band of a source image to be considered as "
5575 905 : "a source alpha band."));
5576 905 : group.add_argument("-nosrcalpha")
5577 905 : .flag()
5578 905 : .store_into(psOptions->bDisableSrcAlpha)
5579 : .help(_("Prevent the alpha band of a source image to be considered "
5580 905 : "as such."));
5581 : }
5582 :
5583 905 : argParser->add_argument("-dstalpha")
5584 905 : .flag()
5585 905 : .store_into(psOptions->bEnableDstAlpha)
5586 : .help(_("Create an output alpha band to identify nodata "
5587 905 : "(unset/transparent) pixels."));
5588 :
5589 : // Parsing of that option is done in a preprocessing stage
5590 905 : argParser->add_argument("-tr")
5591 1810 : .metavar("<xres> <yres>|square")
5592 905 : .help(_("Target resolution."));
5593 :
5594 905 : argParser->add_argument("-ts")
5595 1810 : .metavar("<width> <height>")
5596 905 : .nargs(2)
5597 905 : .scan<'i', int>()
5598 905 : .help(_("Set output file size in pixels and lines."));
5599 :
5600 905 : argParser->add_argument("-te")
5601 1810 : .metavar("<xmin> <ymin> <xmax> <ymax>")
5602 905 : .nargs(4)
5603 905 : .scan<'g', double>()
5604 905 : .help(_("Set georeferenced extents of output file to be created."));
5605 :
5606 905 : argParser->add_argument("-te_srs")
5607 1810 : .metavar("<srs_def>")
5608 : .action(
5609 9 : [psOptions](const std::string &s)
5610 : {
5611 3 : if (!IsValidSRS(s.c_str()))
5612 : {
5613 0 : throw std::invalid_argument("Invalid SRS for -te_srs");
5614 : }
5615 3 : psOptions->osTE_SRS = s;
5616 3 : psOptions->bCreateOutput = true;
5617 908 : })
5618 905 : .help(_("Set source spatial reference."));
5619 :
5620 905 : argParser->add_argument("-r")
5621 : .metavar("near|bilinear|cubic|cubicspline|lanczos|average|rms|mode|min|"
5622 1810 : "max|med|q1|q3|sum")
5623 : .action(
5624 1019 : [psOptions](const std::string &s)
5625 : {
5626 510 : GetResampleAlg(s.c_str(), psOptions->eResampleAlg,
5627 : /*bThrow=*/true);
5628 509 : psOptions->bResampleAlgSpecifiedByUser = true;
5629 905 : })
5630 905 : .help(_("Resampling method to use."));
5631 :
5632 905 : argParser->add_output_type_argument(psOptions->eOutputType);
5633 :
5634 : ///////////////////////////////////////////////////////////////////////
5635 905 : argParser->add_group("Advanced options");
5636 :
5637 32 : const auto CheckSingleMethod = [psOptions]()
5638 : {
5639 : const char *pszMethod =
5640 16 : psOptions->aosTransformerOptions.FetchNameValue("METHOD");
5641 16 : if (pszMethod)
5642 0 : CPLError(CE_Warning, CPLE_IllegalArg,
5643 : "Warning: only one METHOD can be used. Method %s is "
5644 : "already defined.",
5645 : pszMethod);
5646 : const char *pszMAX_GCP_ORDER =
5647 16 : psOptions->aosTransformerOptions.FetchNameValue("MAX_GCP_ORDER");
5648 16 : if (pszMAX_GCP_ORDER)
5649 0 : CPLError(CE_Warning, CPLE_IllegalArg,
5650 : "Warning: only one METHOD can be used. -order %s "
5651 : "option was specified, so it is likely that "
5652 : "GCP_POLYNOMIAL was implied.",
5653 : pszMAX_GCP_ORDER);
5654 921 : };
5655 :
5656 905 : argParser->add_argument("-wo")
5657 1810 : .metavar("<NAME>=<VALUE>")
5658 905 : .append()
5659 119 : .action([psOptions](const std::string &s)
5660 1024 : { psOptions->aosWarpOptions.AddString(s.c_str()); })
5661 905 : .help(_("Warping option(s)."));
5662 :
5663 905 : argParser->add_argument("-multi")
5664 905 : .flag()
5665 905 : .store_into(psOptions->bMulti)
5666 905 : .help(_("Multithreaded input/output."));
5667 :
5668 905 : argParser->add_argument("-s_coord_epoch")
5669 1810 : .metavar("<epoch>")
5670 : .action(
5671 0 : [psOptions](const std::string &s)
5672 : {
5673 : psOptions->aosTransformerOptions.SetNameValue(
5674 0 : "SRC_COORDINATE_EPOCH", s.c_str());
5675 905 : })
5676 : .help(_("Assign a coordinate epoch, linked with the source SRS when "
5677 905 : "-s_srs is used."));
5678 :
5679 905 : argParser->add_argument("-t_coord_epoch")
5680 1810 : .metavar("<epoch>")
5681 : .action(
5682 0 : [psOptions](const std::string &s)
5683 : {
5684 : psOptions->aosTransformerOptions.SetNameValue(
5685 0 : "DST_COORDINATE_EPOCH", s.c_str());
5686 905 : })
5687 : .help(_("Assign a coordinate epoch, linked with the output SRS when "
5688 905 : "-t_srs is used."));
5689 :
5690 905 : argParser->add_argument("-ct")
5691 1810 : .metavar("<string>")
5692 : .action(
5693 4 : [psOptions](const std::string &s)
5694 : {
5695 : psOptions->aosTransformerOptions.SetNameValue(
5696 4 : "COORDINATE_OPERATION", s.c_str());
5697 905 : })
5698 905 : .help(_("Set a coordinate transformation."));
5699 :
5700 : {
5701 905 : auto &group = argParser->add_mutually_exclusive_group();
5702 905 : group.add_argument("-tps")
5703 905 : .flag()
5704 : .action(
5705 10 : [psOptions, CheckSingleMethod](const std::string &)
5706 : {
5707 5 : CheckSingleMethod();
5708 : psOptions->aosTransformerOptions.SetNameValue("METHOD",
5709 5 : "GCP_TPS");
5710 905 : })
5711 : .help(_("Force use of thin plate spline transformer based on "
5712 905 : "available GCPs."));
5713 :
5714 905 : group.add_argument("-rpc")
5715 905 : .flag()
5716 : .action(
5717 4 : [psOptions, CheckSingleMethod](const std::string &)
5718 : {
5719 2 : CheckSingleMethod();
5720 : psOptions->aosTransformerOptions.SetNameValue("METHOD",
5721 2 : "RPC");
5722 905 : })
5723 905 : .help(_("Force use of RPCs."));
5724 :
5725 905 : group.add_argument("-geoloc")
5726 905 : .flag()
5727 : .action(
5728 18 : [psOptions, CheckSingleMethod](const std::string &)
5729 : {
5730 9 : CheckSingleMethod();
5731 : psOptions->aosTransformerOptions.SetNameValue(
5732 9 : "METHOD", "GEOLOC_ARRAY");
5733 905 : })
5734 905 : .help(_("Force use of Geolocation Arrays."));
5735 : }
5736 :
5737 905 : argParser->add_argument("-order")
5738 1810 : .metavar("<1|2|3>")
5739 905 : .choices("1", "2", "3")
5740 : .action(
5741 0 : [psOptions](const std::string &s)
5742 : {
5743 : const char *pszMethod =
5744 0 : psOptions->aosTransformerOptions.FetchNameValue("METHOD");
5745 0 : if (pszMethod)
5746 0 : CPLError(
5747 : CE_Warning, CPLE_IllegalArg,
5748 : "Warning: only one METHOD can be used. Method %s is "
5749 : "already defined",
5750 : pszMethod);
5751 : psOptions->aosTransformerOptions.SetNameValue("MAX_GCP_ORDER",
5752 0 : s.c_str());
5753 905 : })
5754 905 : .help(_("Order of polynomial used for GCP warping."));
5755 :
5756 : // Parsing of that option is done in a preprocessing stage
5757 905 : argParser->add_argument("-refine_gcps")
5758 1810 : .metavar("<tolerance> [<minimum_gcps>]")
5759 905 : .help(_("Refines the GCPs by automatically eliminating outliers."));
5760 :
5761 905 : argParser->add_argument("-to")
5762 1810 : .metavar("<NAME>=<VALUE>")
5763 905 : .append()
5764 44 : .action([psOptions](const std::string &s)
5765 949 : { psOptions->aosTransformerOptions.AddString(s.c_str()); })
5766 905 : .help(_("Transform option(s)."));
5767 :
5768 905 : argParser->add_argument("-et")
5769 1810 : .metavar("<err_threshold>")
5770 905 : .store_into(psOptions->dfErrorThreshold)
5771 : .action(
5772 12 : [psOptions](const std::string &)
5773 : {
5774 : psOptions->aosWarpOptions.AddString(CPLSPrintf(
5775 12 : "ERROR_THRESHOLD=%.16g", psOptions->dfErrorThreshold));
5776 905 : })
5777 905 : .help(_("Error threshold."));
5778 :
5779 905 : argParser->add_argument("-wm")
5780 1810 : .metavar("<memory_in_mb>")
5781 : .action(
5782 66 : [psOptions](const std::string &s)
5783 : {
5784 34 : bool bUnitSpecified = false;
5785 : GIntBig nBytes;
5786 34 : if (CPLParseMemorySize(s.c_str(), &nBytes, &bUnitSpecified) ==
5787 : CE_None)
5788 : {
5789 32 : if (!bUnitSpecified && nBytes < 10000)
5790 : {
5791 6 : nBytes *= (1024 * 1024);
5792 : }
5793 32 : psOptions->dfWarpMemoryLimit = static_cast<double>(nBytes);
5794 : }
5795 : else
5796 : {
5797 2 : throw std::invalid_argument("Failed to parse value of -wm");
5798 : }
5799 937 : })
5800 905 : .help(_("Set max warp memory."));
5801 :
5802 905 : argParser->add_argument("-srcnodata")
5803 1810 : .metavar("\"<value>[ <value>]...\"")
5804 905 : .store_into(psOptions->osSrcNodata)
5805 905 : .help(_("Nodata masking values for input bands."));
5806 :
5807 905 : argParser->add_argument("-dstnodata")
5808 1810 : .metavar("\"<value>[ <value>]...\"")
5809 905 : .store_into(psOptions->osDstNodata)
5810 905 : .help(_("Nodata masking values for output bands."));
5811 :
5812 905 : argParser->add_argument("-tap")
5813 905 : .flag()
5814 905 : .store_into(psOptions->bTargetAlignedPixels)
5815 905 : .help(_("Force target aligned pixels."));
5816 :
5817 905 : argParser->add_argument("-wt")
5818 1810 : .metavar("Byte|Int8|[U]Int{16|32|64}|CInt{16|32}|[C]Float{32|64}")
5819 : .action(
5820 0 : [psOptions](const std::string &s)
5821 : {
5822 0 : psOptions->eWorkingType = GDALGetDataTypeByName(s.c_str());
5823 0 : if (psOptions->eWorkingType == GDT_Unknown)
5824 : {
5825 : throw std::invalid_argument(
5826 0 : std::string("Unknown output pixel type: ").append(s));
5827 : }
5828 905 : })
5829 905 : .help(_("Working data type."));
5830 :
5831 : // Non-documented alias of -r nearest
5832 905 : argParser->add_argument("-rn")
5833 905 : .flag()
5834 905 : .hidden()
5835 1 : .action([psOptions](const std::string &)
5836 905 : { psOptions->eResampleAlg = GRA_NearestNeighbour; })
5837 905 : .help(_("Nearest neighbour resampling."));
5838 :
5839 : // Non-documented alias of -r bilinear
5840 905 : argParser->add_argument("-rb")
5841 905 : .flag()
5842 905 : .hidden()
5843 2 : .action([psOptions](const std::string &)
5844 905 : { psOptions->eResampleAlg = GRA_Bilinear; })
5845 905 : .help(_("Bilinear resampling."));
5846 :
5847 : // Non-documented alias of -r cubic
5848 905 : argParser->add_argument("-rc")
5849 905 : .flag()
5850 905 : .hidden()
5851 1 : .action([psOptions](const std::string &)
5852 905 : { psOptions->eResampleAlg = GRA_Cubic; })
5853 905 : .help(_("Cubic resampling."));
5854 :
5855 : // Non-documented alias of -r cubicspline
5856 905 : argParser->add_argument("-rcs")
5857 905 : .flag()
5858 905 : .hidden()
5859 1 : .action([psOptions](const std::string &)
5860 905 : { psOptions->eResampleAlg = GRA_CubicSpline; })
5861 905 : .help(_("Cubic spline resampling."));
5862 :
5863 : // Non-documented alias of -r lanczos
5864 905 : argParser->add_argument("-rl")
5865 905 : .flag()
5866 905 : .hidden()
5867 0 : .action([psOptions](const std::string &)
5868 905 : { psOptions->eResampleAlg = GRA_Lanczos; })
5869 905 : .help(_("Lanczos resampling."));
5870 :
5871 : // Non-documented alias of -r average
5872 905 : argParser->add_argument("-ra")
5873 905 : .flag()
5874 905 : .hidden()
5875 0 : .action([psOptions](const std::string &)
5876 905 : { psOptions->eResampleAlg = GRA_Average; })
5877 905 : .help(_("Average resampling."));
5878 :
5879 : // Non-documented alias of -r rms
5880 905 : argParser->add_argument("-rrms")
5881 905 : .flag()
5882 905 : .hidden()
5883 0 : .action([psOptions](const std::string &)
5884 905 : { psOptions->eResampleAlg = GRA_RMS; })
5885 905 : .help(_("RMS resampling."));
5886 :
5887 : // Non-documented alias of -r mode
5888 905 : argParser->add_argument("-rm")
5889 905 : .flag()
5890 905 : .hidden()
5891 0 : .action([psOptions](const std::string &)
5892 905 : { psOptions->eResampleAlg = GRA_Mode; })
5893 905 : .help(_("Mode resampling."));
5894 :
5895 905 : argParser->add_argument("-cutline")
5896 1810 : .metavar("<datasource>|<WKT>")
5897 905 : .store_into(psOptions->osCutlineDSNameOrWKT)
5898 : .help(_("Enable use of a blend cutline from the name of a vector "
5899 905 : "dataset or a WKT geometry."));
5900 :
5901 905 : argParser->add_argument("-cutline_srs")
5902 1810 : .metavar("<srs_def>")
5903 : .action(
5904 4 : [psOptions](const std::string &s)
5905 : {
5906 2 : if (!IsValidSRS(s.c_str()))
5907 : {
5908 0 : throw std::invalid_argument("Invalid SRS for -cutline_srs");
5909 : }
5910 2 : psOptions->osCutlineSRS = s;
5911 907 : })
5912 905 : .help(_("Sets/overrides cutline SRS."));
5913 :
5914 905 : argParser->add_argument("-cwhere")
5915 1810 : .metavar("<expression>")
5916 905 : .store_into(psOptions->osCWHERE)
5917 905 : .help(_("Restrict desired cutline features based on attribute query."));
5918 :
5919 : {
5920 905 : auto &group = argParser->add_mutually_exclusive_group();
5921 905 : group.add_argument("-cl")
5922 1810 : .metavar("<layername>")
5923 905 : .store_into(psOptions->osCLayer)
5924 905 : .help(_("Select the named layer from the cutline datasource."));
5925 :
5926 905 : group.add_argument("-csql")
5927 1810 : .metavar("<query>")
5928 905 : .store_into(psOptions->osCSQL)
5929 905 : .help(_("Select cutline features using an SQL query."));
5930 : }
5931 :
5932 905 : argParser->add_argument("-cblend")
5933 1810 : .metavar("<distance>")
5934 : .action(
5935 0 : [psOptions](const std::string &s) {
5936 : psOptions->aosWarpOptions.SetNameValue("CUTLINE_BLEND_DIST",
5937 0 : s.c_str());
5938 905 : })
5939 : .help(_(
5940 905 : "Set a blend distance to use to blend over cutlines (in pixels)."));
5941 :
5942 905 : argParser->add_argument("-crop_to_cutline")
5943 905 : .flag()
5944 : .action(
5945 18 : [psOptions](const std::string &)
5946 : {
5947 18 : psOptions->bCropToCutline = true;
5948 18 : psOptions->bCreateOutput = true;
5949 905 : })
5950 : .help(_("Crop the extent of the target dataset to the extent of the "
5951 905 : "cutline."));
5952 :
5953 905 : argParser->add_argument("-nomd")
5954 905 : .flag()
5955 : .action(
5956 0 : [psOptions](const std::string &)
5957 : {
5958 0 : psOptions->bCopyMetadata = false;
5959 0 : psOptions->bCopyBandInfo = false;
5960 905 : })
5961 905 : .help(_("Do not copy metadata."));
5962 :
5963 905 : argParser->add_argument("-cvmd")
5964 1810 : .metavar("<meta_conflict_value>")
5965 905 : .store_into(psOptions->osMDConflictValue)
5966 : .help(_("Value to set metadata items that conflict between source "
5967 905 : "datasets."));
5968 :
5969 905 : argParser->add_argument("-setci")
5970 905 : .flag()
5971 905 : .store_into(psOptions->bSetColorInterpretation)
5972 : .help(_("Set the color interpretation of the bands of the target "
5973 905 : "dataset from the source dataset."));
5974 :
5975 : argParser->add_open_options_argument(
5976 905 : psOptionsForBinary ? &(psOptionsForBinary->aosOpenOptions) : nullptr);
5977 :
5978 905 : argParser->add_argument("-doo")
5979 1810 : .metavar("<NAME>=<VALUE>")
5980 905 : .append()
5981 : .action(
5982 0 : [psOptionsForBinary](const std::string &s)
5983 : {
5984 0 : if (psOptionsForBinary)
5985 0 : psOptionsForBinary->aosDestOpenOptions.AddString(s.c_str());
5986 905 : })
5987 905 : .help(_("Open option(s) for output dataset."));
5988 :
5989 905 : argParser->add_argument("-ovr")
5990 1810 : .metavar("<level>|AUTO|AUTO-<n>|NONE")
5991 : .action(
5992 24 : [psOptions](const std::string &s)
5993 : {
5994 12 : const char *pszOvLevel = s.c_str();
5995 12 : if (EQUAL(pszOvLevel, "AUTO"))
5996 1 : psOptions->nOvLevel = OVR_LEVEL_AUTO;
5997 11 : else if (STARTS_WITH_CI(pszOvLevel, "AUTO-"))
5998 1 : psOptions->nOvLevel =
5999 1 : OVR_LEVEL_AUTO - atoi(pszOvLevel + strlen("AUTO-"));
6000 10 : else if (EQUAL(pszOvLevel, "NONE"))
6001 5 : psOptions->nOvLevel = OVR_LEVEL_NONE;
6002 5 : else if (CPLGetValueType(pszOvLevel) == CPL_VALUE_INTEGER)
6003 5 : psOptions->nOvLevel = atoi(pszOvLevel);
6004 : else
6005 : {
6006 : throw std::invalid_argument(CPLSPrintf(
6007 0 : "Invalid value '%s' for -ov option", pszOvLevel));
6008 : }
6009 917 : })
6010 905 : .help(_("Specify which overview level of source files must be used."));
6011 :
6012 : {
6013 905 : auto &group = argParser->add_mutually_exclusive_group();
6014 905 : group.add_argument("-vshift")
6015 905 : .flag()
6016 905 : .store_into(psOptions->bVShift)
6017 905 : .help(_("Force the use of vertical shift."));
6018 905 : group.add_argument("-novshift", "-novshiftgrid")
6019 905 : .flag()
6020 905 : .store_into(psOptions->bNoVShift)
6021 905 : .help(_("Disable the use of vertical shift."));
6022 : }
6023 :
6024 : argParser->add_input_format_argument(
6025 : psOptionsForBinary ? &psOptionsForBinary->aosAllowedInputDrivers
6026 905 : : nullptr);
6027 :
6028 905 : argParser->add_argument("-b", "-srcband")
6029 1810 : .metavar("<band>")
6030 905 : .append()
6031 905 : .store_into(psOptions->anSrcBands)
6032 905 : .help(_("Specify input band(s) number to warp."));
6033 :
6034 905 : argParser->add_argument("-dstband")
6035 1810 : .metavar("<band>")
6036 905 : .append()
6037 905 : .store_into(psOptions->anDstBands)
6038 905 : .help(_("Specify the output band number in which to warp."));
6039 :
6040 905 : if (psOptionsForBinary)
6041 : {
6042 104 : argParser->add_argument("src_dataset_name")
6043 208 : .metavar("<src_dataset_name>")
6044 104 : .nargs(argparse::nargs_pattern::at_least_one)
6045 105 : .action([psOptionsForBinary](const std::string &s)
6046 209 : { psOptionsForBinary->aosSrcFiles.AddString(s.c_str()); })
6047 104 : .help(_("Input dataset(s)."));
6048 :
6049 104 : argParser->add_argument("dst_dataset_name")
6050 208 : .metavar("<dst_dataset_name>")
6051 104 : .store_into(psOptionsForBinary->osDstFilename)
6052 104 : .help(_("Output dataset."));
6053 : }
6054 :
6055 1810 : return argParser;
6056 : }
6057 :
6058 : /************************************************************************/
6059 : /* GDALWarpAppGetParserUsage() */
6060 : /************************************************************************/
6061 :
6062 2 : std::string GDALWarpAppGetParserUsage()
6063 : {
6064 : try
6065 : {
6066 4 : GDALWarpAppOptions sOptions;
6067 4 : GDALWarpAppOptionsForBinary sOptionsForBinary;
6068 : auto argParser =
6069 4 : GDALWarpAppOptionsGetParser(&sOptions, &sOptionsForBinary);
6070 2 : return argParser->usage();
6071 : }
6072 0 : catch (const std::exception &err)
6073 : {
6074 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
6075 0 : err.what());
6076 0 : return std::string();
6077 : }
6078 : }
6079 :
6080 : /************************************************************************/
6081 : /* GDALWarpAppOptionsNew() */
6082 : /************************************************************************/
6083 :
6084 : #ifndef CheckHasEnoughAdditionalArgs_defined
6085 : #define CheckHasEnoughAdditionalArgs_defined
6086 :
6087 64 : static bool CheckHasEnoughAdditionalArgs(CSLConstList papszArgv, int i,
6088 : int nExtraArg, int nArgc)
6089 : {
6090 64 : if (i + nExtraArg >= nArgc)
6091 : {
6092 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6093 0 : "%s option requires %d argument%s", papszArgv[i], nExtraArg,
6094 : nExtraArg == 1 ? "" : "s");
6095 0 : return false;
6096 : }
6097 64 : return true;
6098 : }
6099 : #endif
6100 :
6101 : #define CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(nExtraArg) \
6102 : if (!CheckHasEnoughAdditionalArgs(papszArgv, i, nExtraArg, nArgc)) \
6103 : { \
6104 : return nullptr; \
6105 : }
6106 :
6107 : /**
6108 : * Allocates a GDALWarpAppOptions struct.
6109 : *
6110 : * @param papszArgv NULL terminated list of options (potentially including
6111 : * filename and open options too), or NULL. The accepted options are the ones of
6112 : * the <a href="/programs/gdalwarp.html">gdalwarp</a> utility.
6113 : * @param psOptionsForBinary (output) may be NULL (and should generally be
6114 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
6115 : * GDALWarpAppOptionsForBinaryNew() prior to this
6116 : * function. Will be filled with potentially present filename, open options,...
6117 : * @return pointer to the allocated GDALWarpAppOptions struct. Must be freed
6118 : * with GDALWarpAppOptionsFree().
6119 : *
6120 : * @since GDAL 2.1
6121 : */
6122 :
6123 : GDALWarpAppOptions *
6124 903 : GDALWarpAppOptionsNew(char **papszArgv,
6125 : GDALWarpAppOptionsForBinary *psOptionsForBinary)
6126 : {
6127 1806 : auto psOptions = std::make_unique<GDALWarpAppOptions>();
6128 :
6129 : /* -------------------------------------------------------------------- */
6130 : /* Pre-processing for custom syntax that ArgumentParser does not */
6131 : /* support. */
6132 : /* -------------------------------------------------------------------- */
6133 :
6134 1806 : CPLStringList aosArgv;
6135 903 : const int nArgc = CSLCount(papszArgv);
6136 5726 : for (int i = 0;
6137 5726 : i < nArgc && papszArgv != nullptr && papszArgv[i] != nullptr; i++)
6138 : {
6139 4823 : if (EQUAL(papszArgv[i], "-refine_gcps"))
6140 : {
6141 0 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
6142 0 : psOptions->aosTransformerOptions.SetNameValue("REFINE_TOLERANCE",
6143 0 : papszArgv[++i]);
6144 0 : if (CPLAtof(papszArgv[i]) < 0)
6145 : {
6146 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6147 : "The tolerance for -refine_gcps may not be negative.");
6148 0 : return nullptr;
6149 : }
6150 0 : if (i < nArgc - 1 && atoi(papszArgv[i + 1]) >= 0 &&
6151 0 : isdigit(static_cast<unsigned char>(papszArgv[i + 1][0])))
6152 : {
6153 0 : psOptions->aosTransformerOptions.SetNameValue(
6154 0 : "REFINE_MINIMUM_GCPS", papszArgv[++i]);
6155 : }
6156 : else
6157 : {
6158 0 : psOptions->aosTransformerOptions.SetNameValue(
6159 0 : "REFINE_MINIMUM_GCPS", "-1");
6160 : }
6161 : }
6162 4823 : else if (EQUAL(papszArgv[i], "-tr") && i + 1 < nArgc &&
6163 65 : EQUAL(papszArgv[i + 1], "square"))
6164 : {
6165 1 : ++i;
6166 1 : psOptions->bSquarePixels = true;
6167 1 : psOptions->bCreateOutput = true;
6168 : }
6169 4822 : else if (EQUAL(papszArgv[i], "-tr"))
6170 : {
6171 64 : CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(2);
6172 64 : psOptions->dfXRes = CPLAtofM(papszArgv[++i]);
6173 64 : psOptions->dfYRes = fabs(CPLAtofM(papszArgv[++i]));
6174 64 : if (psOptions->dfXRes == 0 || psOptions->dfYRes == 0)
6175 : {
6176 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6177 : "Wrong value for -tr parameters.");
6178 0 : return nullptr;
6179 : }
6180 64 : psOptions->bCreateOutput = true;
6181 : }
6182 : // argparser will be confused if the value of a string argument
6183 : // starts with a negative sign.
6184 4758 : else if (EQUAL(papszArgv[i], "-srcnodata") && i + 1 < nArgc)
6185 : {
6186 25 : ++i;
6187 25 : psOptions->osSrcNodata = papszArgv[i];
6188 : }
6189 : // argparser will be confused if the value of a string argument
6190 : // starts with a negative sign.
6191 4733 : else if (EQUAL(papszArgv[i], "-dstnodata") && i + 1 < nArgc)
6192 : {
6193 43 : ++i;
6194 43 : psOptions->osDstNodata = papszArgv[i];
6195 : }
6196 : else
6197 : {
6198 4690 : aosArgv.AddString(papszArgv[i]);
6199 : }
6200 : }
6201 :
6202 : try
6203 : {
6204 : auto argParser =
6205 1806 : GDALWarpAppOptionsGetParser(psOptions.get(), psOptionsForBinary);
6206 :
6207 903 : argParser->parse_args_without_binary_name(aosArgv.List());
6208 :
6209 1034 : if (auto oTS = argParser->present<std::vector<int>>("-ts"))
6210 : {
6211 135 : psOptions->nForcePixels = (*oTS)[0];
6212 135 : psOptions->nForceLines = (*oTS)[1];
6213 135 : psOptions->bCreateOutput = true;
6214 : }
6215 :
6216 1023 : if (auto oTE = argParser->present<std::vector<double>>("-te"))
6217 : {
6218 124 : psOptions->dfMinX = (*oTE)[0];
6219 124 : psOptions->dfMinY = (*oTE)[1];
6220 124 : psOptions->dfMaxX = (*oTE)[2];
6221 124 : psOptions->dfMaxY = (*oTE)[3];
6222 124 : psOptions->bCreateOutput = true;
6223 : }
6224 :
6225 904 : if (!psOptions->anDstBands.empty() &&
6226 5 : psOptions->anSrcBands.size() != psOptions->anDstBands.size())
6227 : {
6228 1 : CPLError(
6229 : CE_Failure, CPLE_IllegalArg,
6230 : "-srcband should be specified as many times as -dstband is");
6231 1 : return nullptr;
6232 : }
6233 916 : else if (!psOptions->anSrcBands.empty() &&
6234 18 : psOptions->anDstBands.empty())
6235 : {
6236 37 : for (int i = 0; i < static_cast<int>(psOptions->anSrcBands.size());
6237 : ++i)
6238 : {
6239 23 : psOptions->anDstBands.push_back(i + 1);
6240 : }
6241 : }
6242 :
6243 1370 : if (!psOptions->osFormat.empty() ||
6244 472 : psOptions->eOutputType != GDT_Unknown)
6245 : {
6246 432 : psOptions->bCreateOutput = true;
6247 : }
6248 :
6249 898 : if (psOptionsForBinary)
6250 100 : psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
6251 :
6252 898 : return psOptions.release();
6253 : }
6254 4 : catch (const std::exception &err)
6255 : {
6256 4 : CPLError(CE_Failure, CPLE_AppDefined, "%s", err.what());
6257 4 : return nullptr;
6258 : }
6259 : }
6260 :
6261 : /************************************************************************/
6262 : /* GetResampleAlg() */
6263 : /************************************************************************/
6264 :
6265 513 : static bool GetResampleAlg(const char *pszResampling,
6266 : GDALResampleAlg &eResampleAlg, bool bThrow)
6267 : {
6268 513 : if (STARTS_WITH_CI(pszResampling, "near"))
6269 81 : eResampleAlg = GRA_NearestNeighbour;
6270 432 : else if (EQUAL(pszResampling, "bilinear"))
6271 107 : eResampleAlg = GRA_Bilinear;
6272 325 : else if (EQUAL(pszResampling, "cubic"))
6273 105 : eResampleAlg = GRA_Cubic;
6274 220 : else if (EQUAL(pszResampling, "cubicspline"))
6275 53 : eResampleAlg = GRA_CubicSpline;
6276 167 : else if (EQUAL(pszResampling, "lanczos"))
6277 50 : eResampleAlg = GRA_Lanczos;
6278 117 : else if (EQUAL(pszResampling, "average"))
6279 65 : eResampleAlg = GRA_Average;
6280 52 : else if (EQUAL(pszResampling, "rms"))
6281 3 : eResampleAlg = GRA_RMS;
6282 49 : else if (EQUAL(pszResampling, "mode"))
6283 19 : eResampleAlg = GRA_Mode;
6284 30 : else if (EQUAL(pszResampling, "max"))
6285 2 : eResampleAlg = GRA_Max;
6286 28 : else if (EQUAL(pszResampling, "min"))
6287 2 : eResampleAlg = GRA_Min;
6288 26 : else if (EQUAL(pszResampling, "med"))
6289 2 : eResampleAlg = GRA_Med;
6290 24 : else if (EQUAL(pszResampling, "q1"))
6291 2 : eResampleAlg = GRA_Q1;
6292 22 : else if (EQUAL(pszResampling, "q3"))
6293 2 : eResampleAlg = GRA_Q3;
6294 20 : else if (EQUAL(pszResampling, "sum"))
6295 19 : eResampleAlg = GRA_Sum;
6296 : else
6297 : {
6298 1 : if (bThrow)
6299 : {
6300 1 : throw std::invalid_argument("Unknown resampling method");
6301 : }
6302 : else
6303 : {
6304 0 : CPLError(CE_Failure, CPLE_IllegalArg,
6305 : "Unknown resampling method: %s.", pszResampling);
6306 0 : return false;
6307 : }
6308 : }
6309 512 : return true;
6310 : }
6311 :
6312 : /************************************************************************/
6313 : /* GDALWarpAppOptionsFree() */
6314 : /************************************************************************/
6315 :
6316 : /**
6317 : * Frees the GDALWarpAppOptions struct.
6318 : *
6319 : * @param psOptions the options struct for GDALWarp().
6320 : *
6321 : * @since GDAL 2.1
6322 : */
6323 :
6324 898 : void GDALWarpAppOptionsFree(GDALWarpAppOptions *psOptions)
6325 : {
6326 898 : delete psOptions;
6327 898 : }
6328 :
6329 : /************************************************************************/
6330 : /* GDALWarpAppOptionsSetProgress() */
6331 : /************************************************************************/
6332 :
6333 : /**
6334 : * Set a progress function.
6335 : *
6336 : * @param psOptions the options struct for GDALWarp().
6337 : * @param pfnProgress the progress callback.
6338 : * @param pProgressData the user data for the progress callback.
6339 : *
6340 : * @since GDAL 2.1
6341 : */
6342 :
6343 121 : void GDALWarpAppOptionsSetProgress(GDALWarpAppOptions *psOptions,
6344 : GDALProgressFunc pfnProgress,
6345 : void *pProgressData)
6346 : {
6347 121 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
6348 121 : psOptions->pProgressData = pProgressData;
6349 121 : if (pfnProgress == GDALTermProgress)
6350 0 : psOptions->bQuiet = false;
6351 121 : }
6352 :
6353 : /************************************************************************/
6354 : /* GDALWarpAppOptionsSetQuiet() */
6355 : /************************************************************************/
6356 :
6357 : /**
6358 : * Set a progress function.
6359 : *
6360 : * @param psOptions the options struct for GDALWarp().
6361 : * @param bQuiet whether GDALWarp() should emit messages on stdout.
6362 : *
6363 : * @since GDAL 2.3
6364 : */
6365 :
6366 93 : void GDALWarpAppOptionsSetQuiet(GDALWarpAppOptions *psOptions, int bQuiet)
6367 : {
6368 93 : psOptions->bQuiet = CPL_TO_BOOL(bQuiet);
6369 93 : }
6370 :
6371 : /************************************************************************/
6372 : /* GDALWarpAppOptionsSetWarpOption() */
6373 : /************************************************************************/
6374 :
6375 : /**
6376 : * Set a warp option
6377 : *
6378 : * @param psOptions the options struct for GDALWarp().
6379 : * @param pszKey key.
6380 : * @param pszValue value.
6381 : *
6382 : * @since GDAL 2.1
6383 : */
6384 :
6385 0 : void GDALWarpAppOptionsSetWarpOption(GDALWarpAppOptions *psOptions,
6386 : const char *pszKey, const char *pszValue)
6387 : {
6388 0 : psOptions->aosWarpOptions.SetNameValue(pszKey, pszValue);
6389 0 : }
6390 :
6391 : #undef CHECK_HAS_ENOUGH_ADDITIONAL_ARGS
|
