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 226 : GDALRasterZonalStatsAlgorithm::GDALRasterZonalStatsAlgorithm(bool bStandalone)
31 : : GDALPipelineStepAlgorithm(
32 : NAME, DESCRIPTION, HELP_URL,
33 0 : ConstructorOptions()
34 226 : .SetStandaloneStep(bStandalone)
35 452 : .SetOutputFormatCreateCapability(GDAL_DCAP_CREATE))
36 : {
37 226 : AddRasterInputArgs(false, false);
38 226 : if (bStandalone)
39 : {
40 209 : AddVectorOutputArgs(false, false);
41 209 : AddProgressArg();
42 : }
43 :
44 226 : constexpr const char *ZONES_BAND_OR_LAYER = "BAND_OR_LAYER";
45 :
46 226 : AddBandArg(&m_bands);
47 452 : AddArg("zones", 0, _("Dataset containing zone definitions"), &m_zones)
48 226 : .SetRequired();
49 : AddArg("zones-band", 0, _("Band from which zones should be read"),
50 452 : &m_zonesBand)
51 226 : .SetMutualExclusionGroup(ZONES_BAND_OR_LAYER);
52 : AddArg("zones-layer", 0, _("Layer from which zones should be read"),
53 452 : &m_zonesLayer)
54 226 : .SetMutualExclusionGroup(ZONES_BAND_OR_LAYER);
55 226 : AddArg("weights", 0, _("Weighting raster dataset"), &m_weights);
56 : AddArg("weights-band", 0, _("Band from which weights should be read"),
57 452 : &m_weightsBand)
58 226 : .SetDefault(1);
59 : AddArg(
60 : "pixels", 0,
61 : _("Method to determine which pixels are included in stat calculation."),
62 452 : &m_pixels)
63 226 : .SetChoices("default", "fractional", "all-touched");
64 452 : AddArg("stat", 0, _("Statistic(s) to compute for each zone"), &m_stats)
65 226 : .SetRequired()
66 226 : .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 226 : "weighted_variance", "weights");
73 : AddArg("include-field", 0,
74 : _("Fields from polygon zones to include in output"),
75 226 : &m_includeFields);
76 : AddArg("include-geom", 0, _("Include polygon zone geometry in the output"),
77 226 : &m_includeZoneGeom);
78 : AddArg("strategy", 0,
79 : _("For polygon zones, whether to iterate over input features or "
80 : "raster chunks"),
81 452 : &m_strategy)
82 226 : .SetChoices("feature", "raster")
83 226 : .SetDefault("feature");
84 : AddMemorySizeArg(&m_memoryBytes, &m_memoryStr, "chunk-size",
85 226 : _("Maximum size of raster chunks read into memory"));
86 226 : }
87 :
88 : /************************************************************************/
89 : /* GDALRasterZonalStatsAlgorithm::RunStep() */
90 : /************************************************************************/
91 :
92 159 : bool GDALRasterZonalStatsAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
93 : void *pProgressData)
94 : {
95 159 : GDALPipelineStepRunContext stepCtxt;
96 159 : stepCtxt.m_pfnProgress = pfnProgress;
97 159 : stepCtxt.m_pProgressData = pProgressData;
98 159 : return RunPreStepPipelineValidations() && RunStep(stepCtxt);
99 : }
100 :
101 : template <typename T>
102 191 : static std::string Join(const T &vec, const std::string &separator)
103 : {
104 382 : std::stringstream ss;
105 191 : bool first = true;
106 512 : for (const auto &val : vec)
107 : {
108 : static_assert(!std::is_floating_point_v<decltype(val)>,
109 : "Precision would be lost.");
110 :
111 321 : if (first)
112 : {
113 191 : first = false;
114 : }
115 : else
116 : {
117 130 : ss << separator;
118 : }
119 321 : ss << val;
120 : }
121 382 : return ss.str();
122 : }
123 :
124 161 : bool GDALRasterZonalStatsAlgorithm::RunStep(GDALPipelineStepRunContext &ctxt)
125 : {
126 161 : 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 161 : GDALDataset *poDstDS = m_outputDataset.GetDatasetRef();
132 161 : std::unique_ptr<GDALDataset> poRetDS;
133 161 : if (!poDstDS)
134 : {
135 161 : std::string outputFilename = m_outputDataset.GetName();
136 161 : if (m_standaloneStep)
137 : {
138 159 : 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 161 : GetGDALDriverManager()->GetDriverByName(m_format.c_str());
162 161 : 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 161 : poRetDS.reset(
171 : poDriver->Create(outputFilename.c_str(), 0, 0, 0, GDT_Unknown,
172 322 : CPLStringList(m_creationOptions).List()));
173 161 : if (!poRetDS)
174 0 : return false;
175 :
176 161 : poDstDS = poRetDS.get();
177 : }
178 : /// End prep of output dataset.
179 :
180 161 : GDALDataset *src = m_inputDataset.front().GetDatasetRef();
181 :
182 161 : CPLStringList aosOptions;
183 161 : if (!m_bands.empty())
184 : {
185 21 : aosOptions.AddNameValue("BANDS", Join(m_bands, ",").c_str());
186 : }
187 161 : if (!m_includeFields.empty())
188 : {
189 : aosOptions.AddNameValue("INCLUDE_FIELDS",
190 9 : Join(m_includeFields, ",").c_str());
191 : }
192 161 : if (m_includeZoneGeom)
193 : {
194 3 : aosOptions.AddNameValue("INCLUDE_GEOM", "YES");
195 : }
196 161 : if (!m_outputLayerName.empty())
197 : {
198 2 : aosOptions.AddNameValue("OUTPUT_LAYER", m_outputLayerName.c_str());
199 : }
200 161 : aosOptions.AddNameValue("PIXEL_INTERSECTION", m_pixels.c_str());
201 161 : if (m_memoryBytes != 0)
202 : {
203 : aosOptions.AddNameValue("RASTER_CHUNK_SIZE_BYTES",
204 161 : std::to_string(m_memoryBytes).c_str());
205 : }
206 161 : aosOptions.AddNameValue("STATS", Join(m_stats, ",").c_str());
207 161 : aosOptions.AddNameValue("STRATEGY", (m_strategy + "_SEQUENTIAL").c_str());
208 161 : if (m_weightsBand != 0)
209 : {
210 : aosOptions.AddNameValue("WEIGHTS_BAND",
211 161 : std::to_string(m_weightsBand).c_str());
212 : }
213 161 : if (m_zonesBand != 0)
214 : {
215 : aosOptions.AddNameValue("ZONES_BAND",
216 1 : std::to_string(m_zonesBand).c_str());
217 : }
218 161 : if (!m_zonesLayer.empty())
219 : {
220 1 : aosOptions.AddNameValue("ZONES_LAYER", m_zonesLayer.c_str());
221 : }
222 179 : for (const std::string &lco : m_layerCreationOptions)
223 : {
224 18 : aosOptions.AddString(std::string("LCO_").append(lco));
225 : }
226 :
227 161 : if (poRetDS)
228 : {
229 161 : m_outputDataset.Set(std::move(poRetDS));
230 : }
231 :
232 161 : return GDALZonalStats(src, m_weights.GetDatasetRef(), zones, poDstDS,
233 161 : aosOptions.List(), ctxt.m_pfnProgress,
234 161 : ctxt.m_pProgressData) == CE_None;
235 : }
236 :
237 : GDALRasterZonalStatsAlgorithmStandalone::
238 : ~GDALRasterZonalStatsAlgorithmStandalone() = default;
239 :
240 : //! @endcond
|
