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