Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster shift-longitude" subcommand
5 : * Author: Dan Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2026, ISciences LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_shift_longitude.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "gdal_priv.h"
17 : #include "gdal_priv_templates.hpp"
18 : #include "../frmts/vrt/gdal_vrt.h"
19 : #include "../frmts/vrt/vrtdataset.h"
20 :
21 : #include <algorithm>
22 :
23 : //! @cond Doxygen_Suppress
24 :
25 : #ifndef _
26 : #define _(x) (x)
27 : #endif
28 :
29 : /************************************************************************/
30 : /*GDALRasterShiftLongitudeAlgorithm::GDALRasterShiftLongitudeAlgorithm()*/
31 : /************************************************************************/
32 :
33 108 : GDALRasterShiftLongitudeAlgorithm::GDALRasterShiftLongitudeAlgorithm(
34 108 : bool bStandalone)
35 108 : : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL, bStandalone)
36 : {
37 216 : AddArg("min-x", 0, _("Sets the minimum longitude value"), &m_minX)
38 108 : .SetRequired();
39 216 : AddArg("max-x", 0, _("Sets the maximum longitude value"), &m_maxX)
40 108 : .SetRequired();
41 108 : AddNodataArg(&m_nodata, /* noneAllowed = */ true, "output-nodata");
42 108 : }
43 :
44 : /************************************************************************/
45 : /* GDALRasterShiftLongitudeAlgorithm::RunStep() */
46 : /************************************************************************/
47 :
48 15 : static double NormalizeLongitude(double lon)
49 : {
50 15 : while (lon < -180.0)
51 3 : lon += 360.0;
52 14 : while (lon > 180.0)
53 2 : lon -= 360.0;
54 12 : return lon;
55 : }
56 :
57 20 : bool GDALRasterShiftLongitudeAlgorithm::RunStep(GDALPipelineStepRunContext &)
58 : {
59 20 : CPLAssert(!m_outputDataset.GetDatasetRef());
60 :
61 20 : const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
62 :
63 20 : GDALGeoTransform srcGT;
64 20 : if (poSrcDS->GetGeoTransform(srcGT) != CE_None)
65 : {
66 1 : ReportError(CE_Failure, CPLE_AppDefined,
67 : "Input dataset does not have a geotransform");
68 1 : return false;
69 : }
70 :
71 19 : if (!srcGT.IsAxisAligned())
72 : {
73 1 : ReportError(CE_Failure, CPLE_AppDefined,
74 : "Input dataset geotransform cannot have a rotation term");
75 1 : return false;
76 : }
77 :
78 18 : const auto nSrcXSize = poSrcDS->GetRasterXSize();
79 :
80 18 : const double dfSrcMinX = srcGT.xorig;
81 18 : const double dfSrcMaxX = srcGT.xorig + srcGT.xscale * nSrcXSize;
82 :
83 : // Snap m_minX and m_maxX to pixel boundaries of the input raster
84 : {
85 : const double dfXOff =
86 36 : std::ceil(std::abs(dfSrcMinX - m_minX) / srcGT.xscale) *
87 18 : (m_minX < dfSrcMinX ? -1 : 1);
88 18 : if (!GDALIsValueInRange<int>(dfXOff))
89 : {
90 1 : CPLError(CE_Failure, CPLE_AppDefined,
91 : "Failed to compute pixel offset between minimum x "
92 : "coordinate of input and output.");
93 1 : return false;
94 : }
95 17 : const int nXOff = static_cast<int>(dfXOff);
96 17 : m_minX = dfSrcMinX + nXOff * srcGT.xscale;
97 : }
98 17 : const double dfDstXSize = std::ceil((m_maxX - m_minX) / srcGT.xscale);
99 17 : if (!GDALIsValueInRange<int>(dfDstXSize))
100 : {
101 1 : CPLError(CE_Failure, CPLE_AppDefined,
102 : "Size of output raster exceeds maximum allowable size.");
103 1 : return false;
104 : }
105 16 : const int nDstXSize = static_cast<int>(dfDstXSize);
106 :
107 16 : if (nDstXSize <= 0)
108 : {
109 2 : ReportError(CE_Failure, CPLE_AppDefined,
110 : "--max-x must be greater than --min-x");
111 2 : return false;
112 : }
113 :
114 14 : m_maxX = m_minX + nDstXSize * srcGT.xscale;
115 14 : const int nYSize = poSrcDS->GetRasterYSize();
116 :
117 28 : auto poDstDS = std::make_unique<VRTDataset>(nDstXSize, nYSize);
118 : {
119 14 : GDALGeoTransform dstGT = srcGT;
120 14 : dstGT.xorig = m_minX;
121 14 : poDstDS->SetGeoTransform(dstGT);
122 : }
123 14 : poDstDS->SetSpatialRef(poSrcDS->GetSpatialRef());
124 :
125 28 : std::vector<GDALRasterWindow> aosSrcWindows;
126 43 : for (int nDstXOff = 0; nDstXOff < nDstXSize;)
127 : {
128 29 : double dfDstChunkMinX = m_minX + nDstXOff * srcGT.xscale;
129 :
130 38 : while (dfDstChunkMinX < dfSrcMinX)
131 : {
132 9 : dfDstChunkMinX += 360;
133 : }
134 43 : while (dfDstChunkMinX >= dfSrcMaxX)
135 : {
136 14 : dfDstChunkMinX -= 360;
137 : }
138 :
139 : const double dfSrcXOff =
140 29 : std::round((dfDstChunkMinX - dfSrcMinX) / srcGT.xscale);
141 29 : if (!GDALIsValueInRange<int>(dfSrcXOff))
142 : {
143 0 : CPLError(CE_Failure, CPLE_AppDefined,
144 : "Failed to calculate source pixels for longitude range "
145 : "beginning at %g",
146 : dfDstChunkMinX);
147 0 : return false;
148 : }
149 29 : const int nSrcXOff = static_cast<int>(dfSrcXOff);
150 :
151 29 : const int nColumnsWanted = nDstXSize - nDstXOff;
152 : // cppcheck-suppress duplicateExpression
153 53 : const int nColumnsAvailable = nSrcXOff < 0 ? 0
154 : : nSrcXOff > nSrcXSize
155 24 : ? 0
156 : : nSrcXSize - nSrcXOff;
157 :
158 29 : if (nColumnsAvailable == 0)
159 : {
160 : const double dfMissingRangeMinX =
161 6 : NormalizeLongitude(dfDstChunkMinX);
162 : const int nMissingPx =
163 6 : std::min(std::abs(nSrcXOff), nDstXSize - nDstXOff);
164 : double dfMissingRangeMaxX =
165 6 : NormalizeLongitude(dfDstChunkMinX + srcGT.xscale * nMissingPx);
166 7 : while (dfMissingRangeMaxX < dfMissingRangeMinX)
167 : {
168 1 : dfMissingRangeMaxX += 360;
169 : }
170 :
171 6 : CPLError(
172 : CE_Warning, CPLE_AppDefined,
173 : "No source data available for output longitude range %g to %g",
174 : dfMissingRangeMinX, dfMissingRangeMaxX);
175 6 : nDstXOff += nMissingPx;
176 :
177 6 : const GDALRasterWindow chunk{-1, -1, nMissingPx, nYSize};
178 6 : aosSrcWindows.push_back(chunk);
179 6 : continue;
180 : }
181 :
182 23 : const GDALRasterWindow chunk{
183 23 : nSrcXOff, 0, std::min(nColumnsWanted, nColumnsAvailable), nYSize};
184 :
185 23 : CPLAssert(chunk.nXSize > 0);
186 :
187 23 : const double dfDstChunkMaxX =
188 23 : m_minX + (nDstXOff + chunk.nXSize) * srcGT.xscale;
189 23 : const double dfSrcChunkMinX = dfSrcMinX + chunk.nXOff * srcGT.xscale;
190 23 : const double dfSrcChunkMaxX =
191 23 : dfSrcMinX + (nSrcXOff + chunk.nXSize) * srcGT.xscale;
192 :
193 23 : CPLDebug("ShiftLongitude", "Src %g - %g, Dst %g - %g", dfSrcChunkMinX,
194 23 : dfSrcChunkMaxX, m_minX + nDstXOff * srcGT.xscale,
195 : dfDstChunkMaxX);
196 :
197 23 : aosSrcWindows.push_back(chunk);
198 23 : nDstXOff += chunk.nXSize;
199 : }
200 :
201 26 : for (int iBand = 1; iBand <= poSrcDS->GetRasterCount(); ++iBand)
202 : {
203 14 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand);
204 14 : poDstDS->AddBand(poSrcBand->GetRasterDataType());
205 : VRTSourcedRasterBand *poDstBand =
206 14 : cpl::down_cast<VRTSourcedRasterBand *>(
207 14 : poDstDS->GetRasterBand(iBand));
208 14 : poDstBand->CopyCommonInfoFrom(poSrcBand);
209 :
210 : bool bHasNoData;
211 14 : double dfSrcNoData{0};
212 :
213 : {
214 : int bSuccess;
215 14 : dfSrcNoData = poSrcBand->GetNoDataValue(&bSuccess);
216 14 : if (!bSuccess)
217 : {
218 12 : dfSrcNoData = VRT_NODATA_UNSET;
219 : }
220 : }
221 :
222 14 : if (!m_nodata.empty())
223 : {
224 8 : bool inexact = false;
225 8 : if (poDstBand->SetNoDataValueAsString(m_nodata.c_str(), &inexact) !=
226 : CE_None)
227 : {
228 2 : CPLError(CE_Failure, CPLE_AppDefined,
229 : "Invalid NoData value: %s", m_nodata.c_str());
230 2 : return false;
231 : }
232 6 : if (inexact)
233 : {
234 0 : CPLError(CE_Failure, CPLE_AppDefined,
235 : "Specified NoData value cannot be represented in the "
236 : "output data type.");
237 0 : return false;
238 : }
239 6 : bHasNoData = true;
240 : }
241 : else
242 : {
243 6 : bHasNoData = GDALCopyNoDataValue(poDstBand, poSrcBand, nullptr);
244 : }
245 :
246 12 : int nDstXOff = 0;
247 39 : for (const auto &chunk : aosSrcWindows)
248 : {
249 27 : const bool bDataAvailable = chunk.nXOff >= 0;
250 :
251 27 : if (bDataAvailable)
252 : {
253 21 : CPLDebug("ShiftLongitude", "Adding source chunk %d,%d %d %d",
254 21 : chunk.nXOff, chunk.nYOff, chunk.nXSize, chunk.nYSize);
255 :
256 : // Scale and offset are propagated to the output band, not handled by the ComplexSource
257 21 : constexpr double dfSourceScale = 1.0;
258 21 : constexpr double dfSourceOffset = 0.0;
259 21 : constexpr int nSrcYOff = 0;
260 21 : constexpr int nDstYOff = 0;
261 42 : if (poDstBand->AddComplexSource(
262 21 : poSrcBand, chunk.nXOff, nSrcYOff, chunk.nXSize,
263 21 : chunk.nYSize, nDstXOff, nDstYOff, chunk.nXSize,
264 21 : chunk.nYSize, dfSourceOffset, dfSourceScale,
265 21 : dfSrcNoData) != CE_None)
266 : {
267 0 : return false;
268 : }
269 : }
270 6 : else if (!bHasNoData)
271 : {
272 0 : CPLErrorOnce(
273 : CE_Warning, CPLE_AppDefined,
274 : "Some regions of the output dataset have no source data "
275 : "available, but a NoData value has not been defined. You "
276 : "can set a value using --output-nodata");
277 : }
278 27 : nDstXOff += chunk.nXSize;
279 : }
280 : }
281 :
282 12 : m_outputDataset.Set(std::move(poDstDS));
283 :
284 12 : return true;
285 : }
286 :
287 : GDALRasterShiftLongitudeAlgorithmStandalone::
288 : ~GDALRasterShiftLongitudeAlgorithmStandalone() = default;
289 :
290 : //! @endcond
|
