Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster update" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_update.h"
14 :
15 : #include "cpl_conv.h"
16 :
17 : #include "gdal_priv.h"
18 : #include "gdal_utils.h"
19 : #include "gdalalg_raster_reproject.h" // for GDALRasterReprojectUtils
20 : #include "gdalalg_raster_overview_refresh.h"
21 : #include "ogr_spatialref.h"
22 :
23 : #include <algorithm>
24 : #include <array>
25 : #include <cmath>
26 : #include <tuple>
27 :
28 : //! @cond Doxygen_Suppress
29 :
30 : #ifndef _
31 : #define _(x) (x)
32 : #endif
33 :
34 : /************************************************************************/
35 : /* GDALRasterUpdateAlgorithm::GDALRasterUpdateAlgorithm() */
36 : /************************************************************************/
37 :
38 54 : GDALRasterUpdateAlgorithm::GDALRasterUpdateAlgorithm(bool standaloneStep)
39 : : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
40 0 : ConstructorOptions()
41 54 : .SetStandaloneStep(standaloneStep)
42 54 : .SetInputDatasetMaxCount(1)
43 54 : .SetAddDefaultArguments(false)
44 108 : .SetInputDatasetAlias("dataset"))
45 : {
46 54 : AddProgressArg();
47 :
48 54 : if (standaloneStep)
49 : {
50 18 : AddRasterInputArgs(/* openForMixedRasterVector = */ false,
51 : /* hiddenForCLI = */ false);
52 : }
53 : else
54 : {
55 36 : AddRasterHiddenInputDatasetArg();
56 : }
57 :
58 54 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER)
59 54 : .SetDatasetInputFlags(GADV_NAME | GADV_OBJECT);
60 :
61 54 : m_update = true;
62 54 : AddUpdateArg(&m_update).SetDefault(true).SetHidden();
63 :
64 108 : AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
65 54 : .SetMutualExclusionGroup("bbox-geometry-like");
66 108 : AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
67 108 : .SetIsCRSArg()
68 54 : .AddHiddenAlias("geometry_srs");
69 :
70 54 : GDALRasterReprojectUtils::AddResamplingArg(this, m_resampling);
71 :
72 54 : GDALRasterReprojectUtils::AddWarpOptTransformOptErrorThresholdArg(
73 54 : this, m_warpOptions, m_transformOptions, m_errorThreshold);
74 :
75 : AddArg("no-update-overviews", 0, _("Do not update existing overviews"),
76 54 : &m_noUpdateOverviews);
77 54 : }
78 :
79 : /************************************************************************/
80 : /* GDALRasterUpdateAlgorithm::RunStep() */
81 : /************************************************************************/
82 :
83 11 : bool GDALRasterUpdateAlgorithm::RunStep(GDALPipelineStepRunContext &stepCtxt)
84 : {
85 11 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
86 11 : CPLAssert(poSrcDS);
87 :
88 11 : auto poDstDS = m_outputDataset.GetDatasetRef();
89 11 : CPLAssert(poDstDS);
90 11 : CPLAssert(poDstDS->GetAccess() == GA_Update);
91 :
92 11 : std::unique_ptr<OGRGeometry> poClipGeom;
93 22 : std::string errMsg;
94 11 : if (!m_geometry.empty())
95 : {
96 2 : std::tie(poClipGeom, errMsg) = GetClipGeometry();
97 2 : if (!poClipGeom)
98 : {
99 1 : ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
100 1 : return false;
101 : }
102 : }
103 :
104 10 : auto poSrcDriver = poSrcDS->GetDriver();
105 10 : auto poDstDriver = poDstDS->GetDriver();
106 19 : if (poSrcDS == poDstDS ||
107 9 : (poSrcDriver && poDstDriver &&
108 9 : !EQUAL(poSrcDriver->GetDescription(), "MEM") &&
109 3 : !EQUAL(poDstDriver->GetDescription(), "MEM") &&
110 3 : strcmp(poSrcDS->GetDescription(), poDstDS->GetDescription()) == 0))
111 : {
112 2 : ReportError(CE_Failure, CPLE_NotSupported,
113 : "Source and destination datasets must be different");
114 2 : return false;
115 : }
116 :
117 16 : CPLStringList aosOptions;
118 8 : if (!m_resampling.empty())
119 : {
120 8 : aosOptions.AddString("-r");
121 8 : aosOptions.AddString(m_resampling.c_str());
122 : }
123 9 : for (const std::string &opt : m_warpOptions)
124 : {
125 1 : aosOptions.AddString("-wo");
126 1 : aosOptions.AddString(opt.c_str());
127 : }
128 9 : for (const std::string &opt : m_transformOptions)
129 : {
130 1 : aosOptions.AddString("-to");
131 1 : aosOptions.AddString(opt.c_str());
132 : }
133 8 : if (std::isfinite(m_errorThreshold))
134 : {
135 1 : aosOptions.AddString("-et");
136 1 : aosOptions.AddString(CPLSPrintf("%.17g", m_errorThreshold));
137 : }
138 :
139 8 : if (poClipGeom)
140 : {
141 1 : aosOptions.AddString("-cutline");
142 1 : aosOptions.AddString(poClipGeom->exportToWkt());
143 : }
144 :
145 8 : bool bOvrCanBeUpdated = false;
146 8 : std::vector<double> overviewRefreshBBox;
147 10 : if (poDstDS->GetRasterBand(1)->GetOverviewCount() > 0 &&
148 2 : !m_noUpdateOverviews)
149 : {
150 2 : GDALGeoTransform gt;
151 2 : const auto poSrcCRS = poSrcDS->GetSpatialRef();
152 2 : const auto poDstCRS = poDstDS->GetSpatialRef();
153 2 : const bool bBothCRS = poSrcCRS && poDstCRS;
154 2 : const bool bBothNoCRS = !poSrcCRS && !poDstCRS;
155 2 : if ((bBothCRS || bBothNoCRS) && poSrcDS->GetGeoTransform(gt) == CE_None)
156 : {
157 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
158 1 : bBothCRS ? OGRCreateCoordinateTransformation(poSrcCRS, poDstCRS)
159 5 : : nullptr);
160 2 : if (bBothNoCRS || poCT)
161 : {
162 2 : const double dfTLX = gt[0];
163 2 : const double dfTLY = gt[3];
164 :
165 2 : double dfTRX = 0;
166 2 : double dfTRY = 0;
167 2 : gt.Apply(poSrcDS->GetRasterXSize(), 0, &dfTRX, &dfTRY);
168 :
169 2 : double dfBLX = 0;
170 2 : double dfBLY = 0;
171 2 : gt.Apply(0, poSrcDS->GetRasterYSize(), &dfBLX, &dfBLY);
172 :
173 2 : double dfBRX = 0;
174 2 : double dfBRY = 0;
175 2 : gt.Apply(poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
176 : &dfBRX, &dfBRY);
177 :
178 : const double dfXMin =
179 2 : std::min(std::min(dfTLX, dfTRX), std::min(dfBLX, dfBRX));
180 : const double dfYMin =
181 2 : std::min(std::min(dfTLY, dfTRY), std::min(dfBLY, dfBRY));
182 : const double dfXMax =
183 2 : std::max(std::max(dfTLX, dfTRX), std::max(dfBLX, dfBRX));
184 : const double dfYMax =
185 2 : std::max(std::max(dfTLY, dfTRY), std::max(dfBLY, dfBRY));
186 2 : double dfOutXMin = dfXMin;
187 2 : double dfOutYMin = dfYMin;
188 2 : double dfOutXMax = dfXMax;
189 2 : double dfOutYMax = dfYMax;
190 3 : if (!poCT || poCT->TransformBounds(
191 : dfXMin, dfYMin, dfXMax, dfYMax, &dfOutXMin,
192 1 : &dfOutYMin, &dfOutXMax, &dfOutYMax, 21))
193 : {
194 2 : bOvrCanBeUpdated = true;
195 2 : CPLDebug("update",
196 : "Refresh overviews from (%f,%f) to (%f,%f)",
197 : dfOutXMin, dfOutYMin, dfOutXMax, dfOutYMax);
198 4 : overviewRefreshBBox = std::vector<double>{
199 2 : dfOutXMin, dfOutYMin, dfOutXMax, dfOutYMax};
200 : }
201 : }
202 : }
203 2 : if (!bOvrCanBeUpdated)
204 : {
205 0 : ReportError(CE_Warning, CPLE_AppDefined,
206 : "Overviews can not be updated");
207 : }
208 : }
209 :
210 8 : bool bOK = false;
211 : GDALWarpAppOptions *psOptions =
212 8 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
213 8 : if (psOptions)
214 : {
215 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)> pScaledData(
216 16 : nullptr, GDALDestroyScaledProgress);
217 8 : auto pfnProgress = stepCtxt.m_pfnProgress;
218 8 : void *pProgressData = stepCtxt.m_pProgressData;
219 8 : if (pfnProgress)
220 : {
221 3 : pScaledData.reset(
222 : GDALCreateScaledProgress(0.0, bOvrCanBeUpdated ? 0.75 : 1.0,
223 : pfnProgress, pProgressData));
224 3 : GDALWarpAppOptionsSetProgress(psOptions, GDALScaledProgress,
225 : pScaledData.get());
226 : }
227 :
228 8 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
229 8 : GDALDatasetH hDstDS = GDALDataset::ToHandle(poDstDS);
230 8 : auto poRetDS = GDALDataset::FromHandle(
231 : GDALWarp(nullptr, hDstDS, 1, &hSrcDS, psOptions, nullptr));
232 8 : GDALWarpAppOptionsFree(psOptions);
233 :
234 8 : bOK = poRetDS != nullptr;
235 8 : if (bOK && bOvrCanBeUpdated)
236 : {
237 2 : GDALRasterOverviewAlgorithmRefresh refresh;
238 2 : refresh.GetArg("dataset")->Set(poRetDS);
239 2 : if (!m_resampling.empty())
240 2 : refresh.GetArg("resampling")->Set(m_resampling);
241 2 : refresh.GetArg("bbox")->Set(overviewRefreshBBox);
242 2 : pScaledData.reset(GDALCreateScaledProgress(0.75, 1.0, pfnProgress,
243 : pProgressData));
244 2 : bOK = refresh.Run(pScaledData ? GDALScaledProgress : nullptr,
245 : pScaledData.get());
246 : }
247 8 : if (bOK && pfnProgress)
248 3 : pfnProgress(1.0, "", pProgressData);
249 : }
250 :
251 8 : return bOK;
252 : }
253 :
254 : /************************************************************************/
255 : /* GDALRasterUpdateAlgorithm::RunImpl() */
256 : /************************************************************************/
257 :
258 9 : bool GDALRasterUpdateAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
259 : void *pProgressData)
260 : {
261 9 : GDALPipelineStepRunContext stepCtxt;
262 9 : stepCtxt.m_pfnProgress = pfnProgress;
263 9 : stepCtxt.m_pProgressData = pProgressData;
264 18 : return RunStep(stepCtxt);
265 : }
266 :
267 : GDALRasterUpdateAlgorithmStandalone::~GDALRasterUpdateAlgorithmStandalone() =
268 : default;
269 :
270 : //! @endcond
|
