Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  "edit" step of "raster pipeline"
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2024, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_raster_edit.h"
      14             : 
      15             : #include "gdal_priv.h"
      16             : #include "gdal_utils.h"
      17             : #include "ogrsf_frmts.h"
      18             : 
      19             : //! @cond Doxygen_Suppress
      20             : 
      21             : #ifndef _
      22             : #define _(x) (x)
      23             : #endif
      24             : 
      25             : /************************************************************************/
      26             : /*                          GetGCPFilename()                            */
      27             : /************************************************************************/
      28             : 
      29          20 : static std::string GetGCPFilename(const std::vector<std::string> &gcps)
      30             : {
      31          20 :     if (gcps.size() == 1 && !gcps[0].empty() && gcps[0][0] == '@')
      32             :     {
      33           8 :         return gcps[0].substr(1);
      34             :     }
      35          12 :     return std::string();
      36             : }
      37             : 
      38             : /************************************************************************/
      39             : /*          GDALRasterEditAlgorithm::GDALRasterEditAlgorithm()          */
      40             : /************************************************************************/
      41             : 
      42         132 : GDALRasterEditAlgorithm::GDALRasterEditAlgorithm(bool standaloneStep)
      43             :     : GDALRasterPipelineStepAlgorithm(
      44             :           NAME, DESCRIPTION, HELP_URL,
      45         132 :           ConstructorOptions().SetAddDefaultArguments(false))
      46             : {
      47         132 :     if (standaloneStep)
      48             :     {
      49          32 :         AddProgressArg();
      50             : 
      51             :         AddArg("dataset", 0,
      52             :                _("Dataset (to be updated in-place, unless --auxiliary)"),
      53          64 :                &m_dataset, GDAL_OF_RASTER | GDAL_OF_UPDATE)
      54          32 :             .SetPositional()
      55          32 :             .SetRequired();
      56             :         AddArg("auxiliary", 0,
      57             :                _("Ask for an auxiliary .aux.xml file to be edited"),
      58          64 :                &m_readOnly)
      59          64 :             .AddHiddenAlias("ro")
      60          32 :             .AddHiddenAlias(GDAL_ARG_NAME_READ_ONLY);
      61             :     }
      62             : 
      63         264 :     AddArg("crs", 0, _("Override CRS (without reprojection)"), &m_overrideCrs)
      64         264 :         .AddHiddenAlias("a_srs")
      65         264 :         .AddHiddenAlias("srs")
      66         132 :         .SetIsCRSArg(/*noneAllowed=*/true);
      67             : 
      68         132 :     AddBBOXArg(&m_bbox);
      69             : 
      70         132 :     AddNodataArg(&m_nodata, /* noneAllowed = */ true);
      71             : 
      72             :     {
      73             :         auto &arg = AddArg("metadata", 0, _("Add/update dataset metadata item"),
      74         264 :                            &m_metadata)
      75         264 :                         .SetMetaVar("<KEY>=<VALUE>")
      76         132 :                         .SetPackedValuesAllowed(false);
      77          21 :         arg.AddValidationAction([this, &arg]()
      78         153 :                                 { return ParseAndValidateKeyValue(arg); });
      79         132 :         arg.AddHiddenAlias("mo");
      80             :     }
      81             : 
      82             :     AddArg("unset-metadata", 0, _("Remove dataset metadata item(s)"),
      83         264 :            &m_unsetMetadata)
      84         132 :         .SetMetaVar("<KEY>");
      85             : 
      86             :     AddArg("unset-metadata-domain", 0, _("Remove dataset metadata domain(s)"),
      87         264 :            &m_unsetMetadataDomain)
      88         132 :         .SetMetaVar("<DOMAIN>");
      89             : 
      90             :     AddArg("gcp", 0,
      91             :            _("Add ground control point, formatted as "
      92             :              "pixel,line,easting,northing[,elevation], or @filename"),
      93         264 :            &m_gcps)
      94         132 :         .SetPackedValuesAllowed(false)
      95             :         .AddValidationAction(
      96          22 :             [this]()
      97             :             {
      98          12 :                 if (GetGCPFilename(m_gcps).empty())
      99             :                 {
     100          16 :                     for (const std::string &gcp : m_gcps)
     101             :                     {
     102             :                         const CPLStringList aosTokens(
     103          10 :                             CSLTokenizeString2(gcp.c_str(), ",", 0));
     104          10 :                         if (aosTokens.size() != 4 && aosTokens.size() != 5)
     105             :                         {
     106           1 :                             ReportError(CE_Failure, CPLE_IllegalArg,
     107             :                                         "Bad format for %s", gcp.c_str());
     108           1 :                             return false;
     109             :                         }
     110          47 :                         for (int i = 0; i < aosTokens.size(); ++i)
     111             :                         {
     112          39 :                             if (CPLGetValueType(aosTokens[i]) ==
     113             :                                 CPL_VALUE_STRING)
     114             :                             {
     115           1 :                                 ReportError(CE_Failure, CPLE_IllegalArg,
     116             :                                             "Bad format for %s", gcp.c_str());
     117           1 :                                 return false;
     118             :                             }
     119             :                         }
     120             :                     }
     121             :                 }
     122          10 :                 return true;
     123         132 :             });
     124             : 
     125         132 :     if (standaloneStep)
     126             :     {
     127          64 :         AddArg("stats", 0, _("Compute statistics, using all pixels"), &m_stats)
     128          32 :             .SetMutualExclusionGroup("stats");
     129             :         AddArg("approx-stats", 0,
     130             :                _("Compute statistics, using a subset of pixels"),
     131          64 :                &m_approxStats)
     132          32 :             .SetMutualExclusionGroup("stats");
     133          32 :         AddArg("hist", 0, _("Compute histogram"), &m_hist);
     134             :     }
     135         132 : }
     136             : 
     137             : /************************************************************************/
     138             : /*           GDALRasterEditAlgorithm::~GDALRasterEditAlgorithm()        */
     139             : /************************************************************************/
     140             : 
     141             : GDALRasterEditAlgorithm::~GDALRasterEditAlgorithm() = default;
     142             : 
     143             : /************************************************************************/
     144             : /*                              ParseGCPs()                             */
     145             : /************************************************************************/
     146             : 
     147           8 : std::vector<gdal::GCP> GDALRasterEditAlgorithm::ParseGCPs() const
     148             : {
     149           8 :     std::vector<gdal::GCP> ret;
     150          16 :     const std::string osGCPFilename = GetGCPFilename(m_gcps);
     151           8 :     if (!osGCPFilename.empty())
     152             :     {
     153             :         auto poDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
     154           4 :             osGCPFilename.c_str(), GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR));
     155           4 :         if (!poDS)
     156           1 :             return ret;
     157           3 :         if (poDS->GetLayerCount() != 1)
     158             :         {
     159           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     160             :                         "GCPs can only be specified for single-layer datasets");
     161           1 :             return ret;
     162             :         }
     163           2 :         auto poLayer = poDS->GetLayer(0);
     164           2 :         const auto poLayerDefn = poLayer->GetLayerDefn();
     165           2 :         int nIdIdx = -1, nInfoIdx = -1, nColIdx = -1, nLineIdx = -1, nXIdx = -1,
     166           2 :             nYIdx = -1, nZIdx = -1;
     167             : 
     168             :         const struct
     169             :         {
     170             :             int &idx;
     171             :             const char *name;
     172             :             bool required;
     173           2 :         } aFields[] = {
     174             :             {nIdIdx, "id", false},     {nInfoIdx, "info", false},
     175             :             {nColIdx, "column", true}, {nLineIdx, "line", true},
     176             :             {nXIdx, "x", true},        {nYIdx, "y", true},
     177             :             {nZIdx, "z", false},
     178           2 :         };
     179             : 
     180          11 :         for (auto &field : aFields)
     181             :         {
     182          10 :             field.idx = poLayerDefn->GetFieldIndex(field.name);
     183          10 :             if (field.idx < 0 && field.required)
     184             :             {
     185           3 :                 ReportError(CE_Failure, CPLE_AppDefined,
     186           1 :                             "Field '%s' cannot be found in '%s'", field.name,
     187           1 :                             poDS->GetDescription());
     188           1 :                 return ret;
     189             :             }
     190             :         }
     191           2 :         for (auto &&poFeature : poLayer)
     192             :         {
     193           2 :             gdal::GCP gcp;
     194           1 :             if (nIdIdx >= 0)
     195           1 :                 gcp.SetId(poFeature->GetFieldAsString(nIdIdx));
     196           1 :             if (nInfoIdx >= 0)
     197           1 :                 gcp.SetInfo(poFeature->GetFieldAsString(nInfoIdx));
     198           1 :             gcp.Pixel() = poFeature->GetFieldAsDouble(nColIdx);
     199           1 :             gcp.Line() = poFeature->GetFieldAsDouble(nLineIdx);
     200           1 :             gcp.X() = poFeature->GetFieldAsDouble(nXIdx);
     201           1 :             gcp.Y() = poFeature->GetFieldAsDouble(nYIdx);
     202           1 :             if (nZIdx >= 0 && poFeature->IsFieldSetAndNotNull(nZIdx))
     203           1 :                 gcp.Z() = poFeature->GetFieldAsDouble(nZIdx);
     204           1 :             ret.push_back(std::move(gcp));
     205             :         }
     206             :     }
     207             :     else
     208             :     {
     209          10 :         for (const std::string &gcpStr : m_gcps)
     210             :         {
     211             :             const CPLStringList aosTokens(
     212          12 :                 CSLTokenizeString2(gcpStr.c_str(), ",", 0));
     213             :             // Verified by validation action
     214           6 :             CPLAssert(aosTokens.size() == 4 || aosTokens.size() == 5);
     215          12 :             gdal::GCP gcp;
     216           6 :             gcp.Pixel() = CPLAtof(aosTokens[0]);
     217           6 :             gcp.Line() = CPLAtof(aosTokens[1]);
     218           6 :             gcp.X() = CPLAtof(aosTokens[2]);
     219           6 :             gcp.Y() = CPLAtof(aosTokens[3]);
     220           6 :             if (aosTokens.size() == 5)
     221           3 :                 gcp.Z() = CPLAtof(aosTokens[4]);
     222           6 :             ret.push_back(std::move(gcp));
     223             :         }
     224             :     }
     225           5 :     return ret;
     226             : }
     227             : 
     228             : /************************************************************************/
     229             : /*                GDALRasterEditAlgorithm::RunStep()                    */
     230             : /************************************************************************/
     231             : 
     232          41 : bool GDALRasterEditAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     233             : {
     234          41 :     GDALDataset *poDS = m_dataset.GetDatasetRef();
     235          41 :     if (poDS)
     236             :     {
     237          25 :         if (poDS->GetAccess() != GA_Update && !m_readOnly)
     238             :         {
     239           1 :             ReportError(CE_Failure, CPLE_AppDefined,
     240             :                         "Dataset should be opened in update mode unless "
     241             :                         "--auxiliary is set");
     242           1 :             return false;
     243             :         }
     244             :     }
     245             :     else
     246             :     {
     247          16 :         const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
     248          16 :         CPLAssert(poSrcDS);
     249          16 :         CPLAssert(m_outputDataset.GetName().empty());
     250          16 :         CPLAssert(!m_outputDataset.GetDatasetRef());
     251             : 
     252          32 :         CPLStringList aosOptions;
     253          16 :         aosOptions.push_back("-of");
     254          16 :         aosOptions.push_back("VRT");
     255             :         GDALTranslateOptions *psOptions =
     256          16 :             GDALTranslateOptionsNew(aosOptions.List(), nullptr);
     257          16 :         GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
     258          16 :         auto poRetDS = GDALDataset::FromHandle(
     259             :             GDALTranslate("", hSrcDS, psOptions, nullptr));
     260          16 :         GDALTranslateOptionsFree(psOptions);
     261          16 :         m_outputDataset.Set(std::unique_ptr<GDALDataset>(poRetDS));
     262          16 :         poDS = m_outputDataset.GetDatasetRef();
     263             :     }
     264             : 
     265          40 :     bool ret = poDS != nullptr;
     266             : 
     267          40 :     if (poDS)
     268             :     {
     269          40 :         if (m_overrideCrs == "null" || m_overrideCrs == "none")
     270             :         {
     271           3 :             if (poDS->SetSpatialRef(nullptr) != CE_None)
     272             :             {
     273           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     274             :                             "SetSpatialRef(%s) failed", m_overrideCrs.c_str());
     275          10 :                 return false;
     276             :             }
     277             :         }
     278          37 :         else if (!m_overrideCrs.empty() && m_gcps.empty())
     279             :         {
     280           3 :             OGRSpatialReference oSRS;
     281           3 :             oSRS.SetFromUserInput(m_overrideCrs.c_str());
     282           3 :             oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     283           3 :             if (poDS->SetSpatialRef(&oSRS) != CE_None)
     284             :             {
     285           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     286             :                             "SetSpatialRef(%s) failed", m_overrideCrs.c_str());
     287           1 :                 return false;
     288             :             }
     289             :         }
     290             : 
     291          38 :         if (!m_bbox.empty())
     292             :         {
     293           4 :             if (poDS->GetRasterXSize() == 0 || poDS->GetRasterYSize() == 0)
     294             :             {
     295           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     296             :                             "Cannot set extent because one of dataset height "
     297             :                             "or width is null");
     298           2 :                 return false;
     299             :             }
     300           3 :             GDALGeoTransform gt;
     301           3 :             gt[0] = m_bbox[0];
     302           3 :             gt[1] = (m_bbox[2] - m_bbox[0]) / poDS->GetRasterXSize();
     303           3 :             gt[2] = 0;
     304           3 :             gt[3] = m_bbox[3];
     305           3 :             gt[4] = 0;
     306           3 :             gt[5] = -(m_bbox[3] - m_bbox[1]) / poDS->GetRasterYSize();
     307           3 :             if (poDS->SetGeoTransform(gt) != CE_None)
     308             :             {
     309           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     310             :                             "Setting extent failed");
     311           1 :                 return false;
     312             :             }
     313             :         }
     314             : 
     315          36 :         if (!m_nodata.empty())
     316             :         {
     317          18 :             for (int i = 0; i < poDS->GetRasterCount(); ++i)
     318             :             {
     319          10 :                 if (EQUAL(m_nodata.c_str(), "none"))
     320           2 :                     poDS->GetRasterBand(i + 1)->DeleteNoDataValue();
     321             :                 else
     322          16 :                     poDS->GetRasterBand(i + 1)->SetNoDataValue(
     323           8 :                         CPLAtof(m_nodata.c_str()));
     324             :             }
     325             :         }
     326             : 
     327          36 :         const CPLStringList aosMD(m_metadata);
     328          47 :         for (const auto &[key, value] : cpl::IterateNameValue(aosMD))
     329             :         {
     330          12 :             if (poDS->SetMetadataItem(key, value) != CE_None)
     331             :             {
     332           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     333             :                             "SetMetadataItem('%s', '%s') failed", key, value);
     334           1 :                 return false;
     335             :             }
     336             :         }
     337             : 
     338          37 :         for (const std::string &key : m_unsetMetadata)
     339             :         {
     340           3 :             if (poDS->SetMetadataItem(key.c_str(), nullptr) != CE_None)
     341             :             {
     342           1 :                 ReportError(CE_Failure, CPLE_AppDefined,
     343             :                             "SetMetadataItem('%s', NULL) failed", key.c_str());
     344           1 :                 return false;
     345             :             }
     346             :         }
     347             : 
     348          35 :         for (const std::string &domain : m_unsetMetadataDomain)
     349             :         {
     350           1 :             if (poDS->SetMetadata(nullptr, domain.c_str()) != CE_None)
     351             :             {
     352           0 :                 ReportError(CE_Failure, CPLE_AppDefined,
     353             :                             "SetMetadata(NULL, '%s') failed", domain.c_str());
     354           0 :                 return false;
     355             :             }
     356             :         }
     357             : 
     358          34 :         if (!m_gcps.empty())
     359             :         {
     360           8 :             const auto gcps = ParseGCPs();
     361           8 :             if (gcps.empty())
     362           3 :                 return false;  // error already emitted by ParseGCPs()
     363             : 
     364           5 :             OGRSpatialReference oSRS;
     365           5 :             if (!m_overrideCrs.empty())
     366             :             {
     367           1 :                 oSRS.SetFromUserInput(m_overrideCrs.c_str());
     368           1 :                 oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     369             :             }
     370             : 
     371           5 :             if (poDS->SetGCPs(static_cast<int>(gcps.size()), gcps[0].c_ptr(),
     372          10 :                               oSRS.IsEmpty() ? nullptr : &oSRS) != CE_None)
     373             :             {
     374           1 :                 ReportError(CE_Failure, CPLE_AppDefined, "Setting GCPs failed");
     375           1 :                 return false;
     376             :             }
     377             :         }
     378             : 
     379          30 :         const int nBands = poDS->GetRasterCount();
     380          30 :         int nCurProgress = 0;
     381          30 :         const double dfTotalProgress =
     382          30 :             ((m_stats || m_approxStats) ? nBands : 0) + (m_hist ? nBands : 0);
     383          30 :         if (m_stats || m_approxStats)
     384             :         {
     385           4 :             for (int i = 0; (i < nBands) && ret; ++i)
     386             :             {
     387           4 :                 void *pScaledProgress = GDALCreateScaledProgress(
     388             :                     nCurProgress / dfTotalProgress,
     389           2 :                     (nCurProgress + 1) / dfTotalProgress, ctxt.m_pfnProgress,
     390             :                     ctxt.m_pProgressData);
     391           2 :                 ++nCurProgress;
     392           2 :                 double dfMin = 0.0;
     393           2 :                 double dfMax = 0.0;
     394           2 :                 double dfMean = 0.0;
     395           2 :                 double dfStdDev = 0.0;
     396           4 :                 ret = poDS->GetRasterBand(i + 1)->ComputeStatistics(
     397           2 :                           m_approxStats, &dfMin, &dfMax, &dfMean, &dfStdDev,
     398             :                           pScaledProgress ? GDALScaledProgress : nullptr,
     399           2 :                           pScaledProgress) == CE_None;
     400           2 :                 GDALDestroyScaledProgress(pScaledProgress);
     401             :             }
     402             :         }
     403             : 
     404          30 :         if (m_hist)
     405             :         {
     406           2 :             for (int i = 0; (i < nBands) && ret; ++i)
     407             :             {
     408           2 :                 void *pScaledProgress = GDALCreateScaledProgress(
     409             :                     nCurProgress / dfTotalProgress,
     410           1 :                     (nCurProgress + 1) / dfTotalProgress, ctxt.m_pfnProgress,
     411             :                     ctxt.m_pProgressData);
     412           1 :                 ++nCurProgress;
     413           1 :                 double dfMin = 0.0;
     414           1 :                 double dfMax = 0.0;
     415           1 :                 int nBucketCount = 0;
     416           1 :                 GUIntBig *panHistogram = nullptr;
     417           2 :                 ret = poDS->GetRasterBand(i + 1)->GetDefaultHistogram(
     418             :                           &dfMin, &dfMax, &nBucketCount, &panHistogram, TRUE,
     419             :                           pScaledProgress ? GDALScaledProgress : nullptr,
     420           1 :                           pScaledProgress) == CE_None;
     421           1 :                 if (ret)
     422             :                 {
     423           2 :                     ret = poDS->GetRasterBand(i + 1)->SetDefaultHistogram(
     424           1 :                               dfMin, dfMax, nBucketCount, panHistogram) ==
     425             :                           CE_None;
     426             :                 }
     427           1 :                 CPLFree(panHistogram);
     428           1 :                 GDALDestroyScaledProgress(pScaledProgress);
     429             :             }
     430             :         }
     431             :     }
     432             : 
     433          30 :     return poDS != nullptr;
     434             : }
     435             : 
     436             : GDALRasterEditAlgorithmStandalone::~GDALRasterEditAlgorithmStandalone() =
     437             :     default;
     438             : 
     439             : //! @endcond