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