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