Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  gdal "raster zonal-stats" subcommand
       5             :  * Author:   Dan Baston
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2025, ISciences LLC
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "gdalalg_raster_zonal_stats.h"
      14             : 
      15             : #include "gdal_alg.h"
      16             : #include "gdal_priv.h"
      17             : #include "gdal_utils.h"
      18             : #include "ogrsf_frmts.h"
      19             : 
      20             : //! @cond Doxygen_Suppress
      21             : 
      22             : #ifndef _
      23             : #define _(x) (x)
      24             : #endif
      25             : 
      26             : /************************************************************************/
      27             : /*    GDALRasterZonalStatsAlgorithm::GDALRasterZonalStatsAlgorithm()    */
      28             : /************************************************************************/
      29             : 
      30         148 : GDALRasterZonalStatsAlgorithm::GDALRasterZonalStatsAlgorithm(bool bStandalone)
      31             :     : GDALPipelineStepAlgorithm(
      32             :           NAME, DESCRIPTION, HELP_URL,
      33           0 :           ConstructorOptions()
      34         148 :               .SetStandaloneStep(bStandalone)
      35         296 :               .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
      36             : {
      37         148 :     if (bStandalone)
      38             :     {
      39         134 :         AddRasterInputArgs(false, false);
      40         134 :         AddVectorOutputArgs(false, false);
      41             :     }
      42             : 
      43         148 :     constexpr const char *ZONES_BAND_OR_LAYER = "BAND_OR_LAYER";
      44             : 
      45         148 :     AddBandArg(&m_bands);
      46         296 :     AddArg("zones", 0, _("Dataset containing zone definitions"), &m_zones)
      47         148 :         .SetRequired();
      48             :     AddArg("zones-band", 0, _("Band from which zones should be read"),
      49         296 :            &m_zonesBand)
      50         148 :         .SetMutualExclusionGroup(ZONES_BAND_OR_LAYER);
      51             :     AddArg("zones-layer", 0, _("Layer from which zones should be read"),
      52         296 :            &m_zonesLayer)
      53         148 :         .SetMutualExclusionGroup(ZONES_BAND_OR_LAYER);
      54         148 :     AddArg("weights", 0, _("Weighting raster dataset"), &m_weights);
      55             :     AddArg("weights-band", 0, _("Band from which weights should be read"),
      56         296 :            &m_weightsBand)
      57         148 :         .SetDefault(1);
      58             :     AddArg(
      59             :         "pixels", 0,
      60             :         _("Method to determine which pixels are included in stat calculation."),
      61         296 :         &m_pixels)
      62         148 :         .SetChoices("default", "fractional", "all-touched");
      63         296 :     AddArg("stat", 0, _("Statistic(s) to compute for each zone"), &m_stats)
      64         148 :         .SetRequired()
      65         148 :         .SetMinValueIncluded(1)
      66             :         .SetChoices("center_x", "center_y", "count", "coverage", "frac", "max",
      67             :                     "max_center_x", "max_center_y", "mean", "median", "min",
      68             :                     "minority", "min_center_x", "min_center_y", "mode", "stdev",
      69             :                     "sum", "unique", "values", "variance", "variety",
      70             :                     "weighted_mean", "weighted_stdev", "weighted_sum",
      71         148 :                     "weighted_variance", "weights");
      72             :     AddArg("include-field", 0,
      73             :            _("Fields from polygon zones to include in output"),
      74         148 :            &m_includeFields);
      75             :     AddArg("strategy", 0,
      76             :            _("For polygon zones, whether to iterate over input features or "
      77             :              "raster chunks"),
      78         296 :            &m_strategy)
      79         148 :         .SetChoices("feature", "raster")
      80         148 :         .SetDefault("feature");
      81             :     AddMemorySizeArg(&m_memoryBytes, &m_memoryStr, "chunk-size",
      82         148 :                      _("Maximum size of raster chunks read into memory"));
      83         148 :     AddProgressArg();
      84         148 : }
      85             : 
      86             : /************************************************************************/
      87             : /*                 GDALRasterZonalStatsAlgorithm::RunStep()             */
      88             : /************************************************************************/
      89             : 
      90         130 : bool GDALRasterZonalStatsAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
      91             :                                             void *pProgressData)
      92             : {
      93         130 :     GDALPipelineStepRunContext stepCtxt;
      94         130 :     stepCtxt.m_pfnProgress = pfnProgress;
      95         130 :     stepCtxt.m_pProgressData = pProgressData;
      96         130 :     return RunPreStepPipelineValidations() && RunStep(stepCtxt);
      97             : }
      98             : 
      99             : template <typename T>
     100         157 : static std::string Join(const T &vec, const std::string &separator)
     101             : {
     102         314 :     std::stringstream ss;
     103         157 :     bool first = true;
     104         406 :     for (const auto &val : vec)
     105             :     {
     106             :         static_assert(!std::is_floating_point_v<decltype(val)>,
     107             :                       "Precision would be lost.");
     108             : 
     109         249 :         if (first)
     110             :         {
     111         157 :             first = false;
     112             :         }
     113             :         else
     114             :         {
     115          92 :             ss << separator;
     116             :         }
     117         249 :         ss << val;
     118             :     }
     119         314 :     return ss.str();
     120             : }
     121             : 
     122         131 : bool GDALRasterZonalStatsAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
     123             : {
     124         131 :     GDALDataset *zones = m_zones.GetDatasetRef();
     125             : 
     126             :     /// Prepare the output dataset.
     127             :     /// This section copy-pasted from gdal raster as-features.
     128             :     /// Avoid duplication.
     129         131 :     GDALDataset *poDstDS = m_outputDataset.GetDatasetRef();
     130         131 :     std::unique_ptr<GDALDataset> poRetDS;
     131         131 :     if (!poDstDS)
     132             :     {
     133         129 :         std::string outputFilename = m_outputDataset.GetName();
     134         129 :         if (m_standaloneStep)
     135             :         {
     136         128 :             if (m_format.empty())
     137             :             {
     138             :                 const auto aosFormats =
     139             :                     CPLStringList(GDALGetOutputDriversForDatasetName(
     140           2 :                         m_outputDataset.GetName().c_str(), GDAL_OF_VECTOR,
     141             :                         /* bSingleMatch = */ true,
     142           2 :                         /* bWarn = */ true));
     143           2 :                 if (aosFormats.size() != 1)
     144             :                 {
     145           0 :                     ReportError(CE_Failure, CPLE_AppDefined,
     146             :                                 "Cannot guess driver for %s",
     147           0 :                                 m_outputDataset.GetName().c_str());
     148           0 :                     return false;
     149             :                 }
     150           2 :                 m_format = aosFormats[0];
     151             :             }
     152             :         }
     153             :         else
     154             :         {
     155           1 :             m_format = "MEM";
     156             :         }
     157             : 
     158             :         auto poDriver =
     159         129 :             GetGDALDriverManager()->GetDriverByName(m_format.c_str());
     160         129 :         if (!poDriver)
     161             :         {
     162             :             // shouldn't happen given checks done in GDALAlgorithm
     163           0 :             ReportError(CE_Failure, CPLE_AppDefined, "Cannot find driver %s",
     164             :                         m_format.c_str());
     165           0 :             return false;
     166             :         }
     167             : 
     168         129 :         poRetDS.reset(
     169             :             poDriver->Create(outputFilename.c_str(), 0, 0, 0, GDT_Unknown,
     170         258 :                              CPLStringList(m_creationOptions).List()));
     171         129 :         if (!poRetDS)
     172           0 :             return false;
     173             : 
     174         129 :         poDstDS = poRetDS.get();
     175             :     }
     176             :     /// End prep of output dataset.
     177             : 
     178         131 :     GDALDataset *src = m_inputDataset.front().GetDatasetRef();
     179             : 
     180         131 :     CPLStringList aosOptions;
     181         131 :     if (!m_bands.empty())
     182             :     {
     183          21 :         aosOptions.AddNameValue("BANDS", Join(m_bands, ",").c_str());
     184             :     }
     185         131 :     if (!m_includeFields.empty())
     186             :     {
     187             :         aosOptions.AddNameValue("INCLUDE_FIELDS",
     188           5 :                                 Join(m_includeFields, ",").c_str());
     189             :     }
     190         131 :     aosOptions.AddNameValue("PIXEL_INTERSECTION", m_pixels.c_str());
     191         131 :     if (m_memoryBytes != 0)
     192             :     {
     193             :         aosOptions.AddNameValue("RASTER_CHUNK_SIZE_BYTES",
     194         131 :                                 std::to_string(m_memoryBytes).c_str());
     195             :     }
     196         131 :     aosOptions.AddNameValue("STATS", Join(m_stats, ",").c_str());
     197         131 :     aosOptions.AddNameValue("STRATEGY", (m_strategy + "_SEQUENTIAL").c_str());
     198         131 :     if (m_weightsBand != 0)
     199             :     {
     200             :         aosOptions.AddNameValue("WEIGHTS_BAND",
     201         131 :                                 std::to_string(m_weightsBand).c_str());
     202             :     }
     203         131 :     if (m_zonesBand != 0)
     204             :     {
     205             :         aosOptions.AddNameValue("ZONES_BAND",
     206           1 :                                 std::to_string(m_zonesBand).c_str());
     207             :     }
     208         131 :     if (!m_zonesLayer.empty())
     209             :     {
     210           1 :         aosOptions.AddNameValue("ZONES_LAYER", m_zonesLayer.c_str());
     211             :     }
     212         149 :     for (const std::string &lco : m_layerCreationOptions)
     213             :     {
     214          18 :         aosOptions.AddString(std::string("LCO_").append(lco));
     215             :     }
     216             : 
     217         131 :     if (poRetDS)
     218             :     {
     219         129 :         m_outputDataset.Set(std::move(poRetDS));
     220             :     }
     221             : 
     222         131 :     return GDALZonalStats(src, m_weights.GetDatasetRef(), zones, poDstDS,
     223         131 :                           aosOptions.List(), ctxt.m_pfnProgress,
     224         131 :                           ctxt.m_pProgressData) == CE_None;
     225             : }
     226             : 
     227             : GDALRasterZonalStatsAlgorithmStandalone::
     228             :     ~GDALRasterZonalStatsAlgorithmStandalone() = default;
     229             : 
     230             : //! @endcond