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