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