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          99 : GDALRasterShiftLongitudeAlgorithm::GDALRasterShiftLongitudeAlgorithm(
      34          99 :     bool bStandalone)
      35          99 :     : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL, bStandalone)
      36             : {
      37         198 :     AddArg("min-x", 0, _("Sets the minimum longitude value"), &m_minX)
      38          99 :         .SetRequired();
      39         198 :     AddArg("max-x", 0, _("Sets the maximum longitude value"), &m_maxX)
      40          99 :         .SetRequired();
      41          99 :     AddNodataArg(&m_nodata, /* noneAllowed = */ true, "output-nodata");
      42          99 : }
      43             : 
      44             : /************************************************************************/
      45             : /*             GDALRasterShiftLongitudeAlgorithm::RunStep()             */
      46             : /************************************************************************/
      47             : 
      48          11 : static double NormalizeLongitude(double lon)
      49             : {
      50          11 :     while (lon < -180.0)
      51           3 :         lon += 360.0;
      52           8 :     while (lon > 180.0)
      53           0 :         lon -= 360.0;
      54           8 :     return lon;
      55             : }
      56             : 
      57          19 : bool GDALRasterShiftLongitudeAlgorithm::RunStep(GDALPipelineStepRunContext &)
      58             : {
      59          19 :     CPLAssert(!m_outputDataset.GetDatasetRef());
      60             : 
      61          19 :     const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
      62             : 
      63          19 :     GDALGeoTransform srcGT;
      64          19 :     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          18 :     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          17 :     const auto nSrcXSize = poSrcDS->GetRasterXSize();
      79             : 
      80          17 :     const double dfSrcMinX = srcGT.xorig;
      81          17 :     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          34 :             std::ceil(std::abs(dfSrcMinX - m_minX) / srcGT.xscale) *
      87          17 :             (m_minX < dfSrcMinX ? -1 : 1);
      88          17 :         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          16 :         const int nXOff = static_cast<int>(dfXOff);
      96          16 :         m_minX = dfSrcMinX + nXOff * srcGT.xscale;
      97             :     }
      98          16 :     const double dfDstXSize = std::ceil((m_maxX - m_minX) / srcGT.xscale);
      99          16 :     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          15 :     const int nDstXSize = static_cast<int>(dfDstXSize);
     106             : 
     107          15 :     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          13 :     m_maxX = m_minX + nDstXSize * srcGT.xscale;
     115          13 :     const int nYSize = poSrcDS->GetRasterYSize();
     116             : 
     117          26 :     auto poDstDS = std::make_unique<VRTDataset>(nDstXSize, nYSize);
     118             :     {
     119          13 :         GDALGeoTransform dstGT = srcGT;
     120          13 :         dstGT.xorig = m_minX;
     121          13 :         poDstDS->SetGeoTransform(dstGT);
     122             :     }
     123             : 
     124          26 :     std::vector<GDALRasterWindow> aosSrcWindows;
     125          38 :     for (int nDstXOff = 0; nDstXOff < nDstXSize;)
     126             :     {
     127          25 :         double dfDstChunkMinX = m_minX + nDstXOff * srcGT.xscale;
     128             : 
     129          31 :         while (dfDstChunkMinX < dfSrcMinX)
     130             :         {
     131           6 :             dfDstChunkMinX += 360;
     132             :         }
     133          38 :         while (dfDstChunkMinX >= dfSrcMaxX)
     134             :         {
     135          13 :             dfDstChunkMinX -= 360;
     136             :         }
     137             : 
     138             :         const double dfSrcXOff =
     139          25 :             std::round((dfDstChunkMinX - dfSrcMinX) / srcGT.xscale);
     140          25 :         if (!GDALIsValueInRange<int>(dfSrcXOff))
     141             :         {
     142           0 :             CPLError(CE_Failure, CPLE_AppDefined,
     143             :                      "Failed to calculate source pixels for longitude range "
     144             :                      "beginning at %g",
     145             :                      dfDstChunkMinX);
     146           0 :             return false;
     147             :         }
     148          25 :         const int nSrcXOff = static_cast<int>(dfSrcXOff);
     149             : 
     150          25 :         if (nSrcXOff < 0)
     151             :         {
     152             :             const double dfMissingRangeMinX =
     153           4 :                 NormalizeLongitude(dfDstChunkMinX);
     154             :             const int nMissingPx =
     155           4 :                 std::min(std::abs(nSrcXOff), nDstXSize - nDstXOff);
     156             :             double dfMissingRangeMaxX =
     157           4 :                 NormalizeLongitude(dfDstChunkMinX + srcGT.xscale * nMissingPx);
     158           5 :             while (dfMissingRangeMaxX < dfMissingRangeMinX)
     159             :             {
     160           1 :                 dfMissingRangeMaxX += 360;
     161             :             }
     162             : 
     163           4 :             CPLError(
     164             :                 CE_Warning, CPLE_AppDefined,
     165             :                 "No source data available for output longitude range %g to %g",
     166             :                 dfMissingRangeMinX, dfMissingRangeMaxX);
     167           4 :             nDstXOff += nMissingPx;
     168             : 
     169           4 :             const GDALRasterWindow chunk{-1, -1, nMissingPx, nYSize};
     170           4 :             aosSrcWindows.push_back(chunk);
     171           4 :             continue;
     172             :         }
     173             : 
     174          21 :         const int nColumnsWanted = nSrcXSize - nSrcXOff;
     175          21 :         const int nColumnsAvailable = nDstXSize - nDstXOff;
     176             : 
     177          21 :         const GDALRasterWindow chunk{
     178          21 :             nSrcXOff, 0, std::min(nColumnsWanted, nColumnsAvailable), nYSize};
     179             : 
     180          21 :         CPLAssert(chunk.nXSize > 0);
     181             : 
     182          21 :         const double dstChunkMaxX = dfDstChunkMinX + chunk.nXSize;
     183          21 :         CPLDebug("ShiftLongitude", "Src %g - %g, Dst %g - %g", dfDstChunkMinX,
     184          21 :                  dstChunkMaxX, m_minX + chunk.nXOff * srcGT.xscale,
     185          21 :                  m_minX + (chunk.nXOff + chunk.nXSize) * srcGT.xscale);
     186             : 
     187          21 :         aosSrcWindows.push_back(chunk);
     188          21 :         nDstXOff += chunk.nXSize;
     189             :     }
     190             : 
     191          24 :     for (int iBand = 1; iBand <= poSrcDS->GetRasterCount(); ++iBand)
     192             :     {
     193          13 :         GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand);
     194          13 :         poDstDS->AddBand(poSrcBand->GetRasterDataType());
     195             :         VRTSourcedRasterBand *poDstBand =
     196          13 :             cpl::down_cast<VRTSourcedRasterBand *>(
     197          13 :                 poDstDS->GetRasterBand(iBand));
     198          13 :         poDstBand->CopyCommonInfoFrom(poSrcBand);
     199             : 
     200             :         bool bHasNoData;
     201          13 :         double dfSrcNoData{0};
     202             : 
     203             :         {
     204             :             int bSuccess;
     205          13 :             dfSrcNoData = poSrcBand->GetNoDataValue(&bSuccess);
     206          13 :             if (!bSuccess)
     207             :             {
     208          11 :                 dfSrcNoData = VRT_NODATA_UNSET;
     209             :             }
     210             :         }
     211             : 
     212          13 :         if (!m_nodata.empty())
     213             :         {
     214           7 :             bool inexact = false;
     215           7 :             if (poDstBand->SetNoDataValueAsString(m_nodata.c_str(), &inexact) !=
     216             :                 CE_None)
     217             :             {
     218           2 :                 CPLError(CE_Failure, CPLE_AppDefined,
     219             :                          "Invalid NoData value: %s", m_nodata.c_str());
     220           2 :                 return false;
     221             :             }
     222           5 :             if (inexact)
     223             :             {
     224           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     225             :                          "Specified NoData value cannot be represented in the "
     226             :                          "output data type.");
     227           0 :                 return false;
     228             :             }
     229           5 :             bHasNoData = true;
     230             :         }
     231             :         else
     232             :         {
     233           6 :             bHasNoData = GDALCopyNoDataValue(poDstBand, poSrcBand, nullptr);
     234             :         }
     235             : 
     236          11 :         int nDstXOff = 0;
     237          34 :         for (const auto &chunk : aosSrcWindows)
     238             :         {
     239          23 :             const bool bDataAvailable = chunk.nXOff >= 0;
     240             : 
     241          23 :             if (bDataAvailable)
     242             :             {
     243          19 :                 CPLDebug("ShiftLongitude", "Adding source chunk %d,%d %d %d",
     244          19 :                          chunk.nXOff, chunk.nYOff, chunk.nXSize, chunk.nYSize);
     245             : 
     246             :                 // Scale and offset are propagated to the output band, not handled by the ComplexSource
     247          19 :                 constexpr double dfSourceScale = 1.0;
     248          19 :                 constexpr double dfSourceOffset = 0.0;
     249          19 :                 constexpr int nSrcYOff = 0;
     250          19 :                 constexpr int nDstYOff = 0;
     251          38 :                 if (poDstBand->AddComplexSource(
     252          19 :                         poSrcBand, chunk.nXOff, nSrcYOff, chunk.nXSize,
     253          19 :                         chunk.nYSize, nDstXOff, nDstYOff, chunk.nXSize,
     254          19 :                         chunk.nYSize, dfSourceOffset, dfSourceScale,
     255          19 :                         dfSrcNoData) != CE_None)
     256             :                 {
     257           0 :                     return false;
     258             :                 }
     259             :             }
     260           4 :             else if (!bHasNoData)
     261             :             {
     262           0 :                 CPLErrorOnce(
     263             :                     CE_Warning, CPLE_AppDefined,
     264             :                     "Some regions of the output dataset have no source data "
     265             :                     "available, but a NoData value has not been defined. You "
     266             :                     "can set a value using --output-nodata");
     267             :             }
     268          23 :             nDstXOff += chunk.nXSize;
     269             :         }
     270             :     }
     271             : 
     272          11 :     m_outputDataset.Set(std::move(poDstDS));
     273             : 
     274          11 :     return true;
     275             : }
     276             : 
     277             : GDALRasterShiftLongitudeAlgorithmStandalone::
     278             :     ~GDALRasterShiftLongitudeAlgorithmStandalone() = default;
     279             : 
     280             : //! @endcond