Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL Utilities
4 : * Purpose: Rasterize OGR shapes into a GDAL raster.
5 : * Author: Frank Warmerdam <warmerdam@pobox.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2008-2015, Even Rouault <even dot rouault at spatialys dot com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_utils.h"
16 : #include "gdal_utils_priv.h"
17 :
18 : #include <cmath>
19 : #include <cstdio>
20 : #include <cstdlib>
21 : #include <cstring>
22 : #include <algorithm>
23 : #include <limits>
24 : #include <vector>
25 :
26 : #include "commonutils.h"
27 : #include "cpl_conv.h"
28 : #include "cpl_error.h"
29 : #include "cpl_progress.h"
30 : #include "cpl_string.h"
31 : #include "gdal.h"
32 : #include "gdal_alg.h"
33 : #include "gdal_priv.h"
34 : #include "ogr_api.h"
35 : #include "ogr_core.h"
36 : #include "ogr_srs_api.h"
37 : #include "gdalargumentparser.h"
38 :
39 : /************************************************************************/
40 : /* GDALRasterizeOptions() */
41 : /************************************************************************/
42 :
43 : struct GDALRasterizeOptions
44 : {
45 : std::vector<int> anBandList{};
46 : std::vector<double> adfBurnValues{};
47 : bool bInverse = false;
48 : std::string osFormat{};
49 : bool b3D = false;
50 : GDALProgressFunc pfnProgress = GDALDummyProgress;
51 : void *pProgressData = nullptr;
52 : std::vector<std::string> aosLayers{};
53 : std::string osSQL{};
54 : std::string osDialect{};
55 : std::string osBurnAttribute{};
56 : std::string osWHERE{};
57 : CPLStringList aosRasterizeOptions{};
58 : CPLStringList aosTO{};
59 : double dfXRes = 0;
60 : double dfYRes = 0;
61 : CPLStringList aosCreationOptions{};
62 : GDALDataType eOutputType = GDT_Unknown;
63 : std::vector<double> adfInitVals{};
64 : std::string osNoData{};
65 : OGREnvelope sEnvelop{};
66 : int nXSize = 0;
67 : int nYSize = 0;
68 : OGRSpatialReference oOutputSRS{};
69 :
70 : bool bTargetAlignedPixels = false;
71 : bool bCreateOutput = false;
72 : };
73 :
74 : /************************************************************************/
75 : /* GDALRasterizeOptionsGetParser() */
76 : /************************************************************************/
77 :
78 : static std::unique_ptr<GDALArgumentParser>
79 29 : GDALRasterizeOptionsGetParser(GDALRasterizeOptions *psOptions,
80 : GDALRasterizeOptionsForBinary *psOptionsForBinary)
81 : {
82 : auto argParser = std::make_unique<GDALArgumentParser>(
83 29 : "gdal_rasterize", /* bForBinary=*/psOptionsForBinary != nullptr);
84 :
85 29 : argParser->add_description(_("Burns vector geometries into a raster."));
86 :
87 29 : argParser->add_epilog(
88 : _("This program burns vector geometries (points, lines, and polygons) "
89 29 : "into the raster band(s) of a raster image."));
90 :
91 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
92 29 : argParser->add_argument("-b")
93 58 : .metavar("<band>")
94 29 : .append()
95 29 : .scan<'i', int>()
96 : //.nargs(argparse::nargs_pattern::at_least_one)
97 29 : .help(_("The band(s) to burn values into."));
98 :
99 29 : argParser->add_argument("-i")
100 29 : .flag()
101 29 : .store_into(psOptions->bInverse)
102 29 : .help(_("Invert rasterization."));
103 :
104 29 : argParser->add_argument("-at")
105 29 : .flag()
106 : .action(
107 1 : [psOptions](const std::string &) {
108 : psOptions->aosRasterizeOptions.SetNameValue("ALL_TOUCHED",
109 1 : "TRUE");
110 29 : })
111 29 : .help(_("Enables the ALL_TOUCHED rasterization option."));
112 :
113 : // Mutually exclusive options: -burn, -3d, -a
114 : {
115 : // Required if options for binary
116 29 : auto &group = argParser->add_mutually_exclusive_group(
117 29 : psOptionsForBinary != nullptr);
118 :
119 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
120 29 : group.add_argument("-burn")
121 58 : .metavar("<value>")
122 29 : .scan<'g', double>()
123 29 : .append()
124 : //.nargs(argparse::nargs_pattern::at_least_one)
125 29 : .help(_("A fixed value to burn into the raster band(s)."));
126 :
127 29 : group.add_argument("-a")
128 58 : .metavar("<attribute_name>")
129 29 : .store_into(psOptions->osBurnAttribute)
130 : .help(_("Name of the field in the input layer to get the burn "
131 29 : "values from."));
132 :
133 29 : group.add_argument("-3d")
134 29 : .flag()
135 29 : .store_into(psOptions->b3D)
136 : .action(
137 3 : [psOptions](const std::string &) {
138 : psOptions->aosRasterizeOptions.SetNameValue(
139 3 : "BURN_VALUE_FROM", "Z");
140 29 : })
141 : .help(_("Indicates that a burn value should be extracted from the "
142 29 : "\"Z\" values of the feature."));
143 : }
144 :
145 29 : argParser->add_argument("-add")
146 29 : .flag()
147 : .action(
148 0 : [psOptions](const std::string &) {
149 0 : psOptions->aosRasterizeOptions.SetNameValue("MERGE_ALG", "ADD");
150 29 : })
151 : .help(_("Instead of burning a new value, this adds the new value to "
152 29 : "the existing raster."));
153 :
154 : // Undocumented
155 29 : argParser->add_argument("-chunkysize")
156 29 : .flag()
157 29 : .hidden()
158 : .action(
159 0 : [psOptions](const std::string &s) {
160 : psOptions->aosRasterizeOptions.SetNameValue("CHUNKYSIZE",
161 0 : s.c_str());
162 29 : });
163 :
164 : // Mutually exclusive -l, -sql
165 : {
166 29 : auto &group = argParser->add_mutually_exclusive_group(false);
167 :
168 29 : group.add_argument("-l")
169 58 : .metavar("<layer_name>")
170 29 : .append()
171 29 : .store_into(psOptions->aosLayers)
172 29 : .help(_("Name of the layer(s) to process."));
173 :
174 29 : group.add_argument("-sql")
175 58 : .metavar("<sql_statement>")
176 29 : .store_into(psOptions->osSQL)
177 : .action(
178 3 : [psOptions](const std::string &sql)
179 : {
180 2 : GByte *pabyRet = nullptr;
181 3 : if (!sql.empty() && sql.at(0) == '@' &&
182 3 : VSIIngestFile(nullptr, sql.substr(1).c_str(), &pabyRet,
183 : nullptr, 1024 * 1024))
184 : {
185 1 : GDALRemoveBOM(pabyRet);
186 1 : char *pszSQLStatement =
187 : reinterpret_cast<char *>(pabyRet);
188 1 : psOptions->osSQL = CPLStrdup(
189 2 : CPLRemoveSQLComments(pszSQLStatement).c_str());
190 1 : VSIFree(pszSQLStatement);
191 : }
192 31 : })
193 : .help(
194 : _("An SQL statement to be evaluated against the datasource to "
195 29 : "produce a virtual layer of features to be burned in."));
196 : }
197 :
198 29 : argParser->add_argument("-where")
199 58 : .metavar("<expression>")
200 29 : .store_into(psOptions->osWHERE)
201 : .help(_("An optional SQL WHERE style query expression to be applied to "
202 : "select features "
203 29 : "to burn in from the input layer(s)."));
204 :
205 29 : argParser->add_argument("-dialect")
206 58 : .metavar("<sql_dialect>")
207 29 : .store_into(psOptions->osDialect)
208 29 : .help(_("The SQL dialect to use for the SQL expression."));
209 :
210 : // Store later
211 29 : argParser->add_argument("-a_nodata")
212 58 : .metavar("<value>")
213 29 : .help(_("Assign a specified nodata value to output bands."));
214 :
215 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
216 29 : argParser->add_argument("-init")
217 58 : .metavar("<value>")
218 29 : .append()
219 : //.nargs(argparse::nargs_pattern::at_least_one)
220 29 : .scan<'g', double>()
221 29 : .help(_("Initialize the output bands to the specified value."));
222 :
223 29 : argParser->add_argument("-a_srs")
224 58 : .metavar("<srs_def>")
225 : .action(
226 2 : [psOptions](const std::string &osOutputSRSDef)
227 : {
228 1 : if (psOptions->oOutputSRS.SetFromUserInput(
229 1 : osOutputSRSDef.c_str()) != OGRERR_NONE)
230 : {
231 : throw std::invalid_argument(
232 0 : std::string("Failed to process SRS definition: ")
233 0 : .append(osOutputSRSDef));
234 : }
235 1 : psOptions->bCreateOutput = true;
236 30 : })
237 29 : .help(_("The spatial reference system to use for the output raster."));
238 :
239 29 : argParser->add_argument("-to")
240 58 : .metavar("<NAME>=<VALUE>")
241 29 : .append()
242 1 : .action([psOptions](const std::string &s)
243 30 : { psOptions->aosTO.AddString(s.c_str()); })
244 29 : .help(_("Set a transformer option."));
245 :
246 : // Store later
247 29 : argParser->add_argument("-te")
248 58 : .metavar("<xmin> <ymin> <xmax> <ymax>")
249 29 : .nargs(4)
250 29 : .scan<'g', double>()
251 29 : .help(_("Set georeferenced extents of output file to be created."));
252 :
253 : // Mutex with tr
254 : {
255 29 : auto &group = argParser->add_mutually_exclusive_group(false);
256 :
257 : // Store later
258 29 : group.add_argument("-tr")
259 58 : .metavar("<xres> <yres>")
260 29 : .nargs(2)
261 29 : .scan<'g', double>()
262 : .help(
263 29 : _("Set output file resolution in target georeferenced units."));
264 :
265 : // Store later
266 29 : auto &arg = group.add_argument("-ts")
267 58 : .metavar("<width> <height>")
268 29 : .nargs(2)
269 29 : .scan<'i', int>()
270 29 : .help(_("Set output file size in pixels and lines."));
271 :
272 29 : argParser->add_hidden_alias_for(arg, "-outsize");
273 : }
274 :
275 29 : argParser->add_argument("-tap")
276 29 : .flag()
277 29 : .store_into(psOptions->bTargetAlignedPixels)
278 0 : .action([psOptions](const std::string &)
279 29 : { psOptions->bCreateOutput = true; })
280 : .help(_("Align the coordinates of the extent to the values of the "
281 29 : "output raster."));
282 :
283 29 : argParser->add_argument("-optim")
284 58 : .metavar("AUTO|VECTOR|RASTER")
285 : .action(
286 3 : [psOptions](const std::string &s) {
287 3 : psOptions->aosRasterizeOptions.SetNameValue("OPTIM", s.c_str());
288 29 : })
289 29 : .help(_("Force the algorithm used."));
290 :
291 29 : argParser->add_creation_options_argument(psOptions->aosCreationOptions)
292 0 : .action([psOptions](const std::string &)
293 29 : { psOptions->bCreateOutput = true; });
294 :
295 29 : argParser->add_output_type_argument(psOptions->eOutputType)
296 2 : .action([psOptions](const std::string &)
297 29 : { psOptions->bCreateOutput = true; });
298 :
299 29 : argParser->add_output_format_argument(psOptions->osFormat)
300 6 : .action([psOptions](const std::string &)
301 29 : { psOptions->bCreateOutput = true; });
302 :
303 : // Written that way so that in library mode, users can still use the -q
304 : // switch, even if it has no effect
305 : argParser->add_quiet_argument(
306 29 : psOptionsForBinary ? &(psOptionsForBinary->bQuiet) : nullptr);
307 :
308 29 : if (psOptionsForBinary)
309 : {
310 :
311 : argParser->add_open_options_argument(
312 9 : psOptionsForBinary->aosOpenOptions);
313 :
314 9 : argParser->add_argument("src_datasource")
315 18 : .metavar("<src_datasource>")
316 9 : .store_into(psOptionsForBinary->osSource)
317 9 : .help(_("Any vector supported readable datasource."));
318 :
319 9 : argParser->add_argument("dst_filename")
320 18 : .metavar("<dst_filename>")
321 9 : .store_into(psOptionsForBinary->osDest)
322 9 : .help(_("The GDAL raster supported output file."));
323 : }
324 :
325 29 : return argParser;
326 : }
327 :
328 : /************************************************************************/
329 : /* GDALRasterizeAppGetParserUsage() */
330 : /************************************************************************/
331 :
332 0 : std::string GDALRasterizeAppGetParserUsage()
333 : {
334 : try
335 : {
336 0 : GDALRasterizeOptions sOptions;
337 0 : GDALRasterizeOptionsForBinary sOptionsForBinary;
338 : auto argParser =
339 0 : GDALRasterizeOptionsGetParser(&sOptions, &sOptionsForBinary);
340 0 : return argParser->usage();
341 : }
342 0 : catch (const std::exception &err)
343 : {
344 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
345 0 : err.what());
346 0 : return std::string();
347 : }
348 : }
349 :
350 : /************************************************************************/
351 : /* InvertGeometries() */
352 : /************************************************************************/
353 :
354 2 : static void InvertGeometries(GDALDatasetH hDstDS,
355 : std::vector<OGRGeometryH> &ahGeometries)
356 :
357 : {
358 2 : OGRMultiPolygon *poInvertMP = new OGRMultiPolygon();
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* Create a ring that is a bit outside the raster dataset. */
362 : /* -------------------------------------------------------------------- */
363 2 : const int brx = GDALGetRasterXSize(hDstDS) + 2;
364 2 : const int bry = GDALGetRasterYSize(hDstDS) + 2;
365 :
366 2 : double adfGeoTransform[6] = {};
367 2 : GDALGetGeoTransform(hDstDS, adfGeoTransform);
368 :
369 2 : auto poUniverseRing = std::make_unique<OGRLinearRing>();
370 :
371 2 : poUniverseRing->addPoint(
372 2 : adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
373 2 : adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);
374 :
375 2 : poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
376 2 : -2 * adfGeoTransform[2],
377 2 : adfGeoTransform[3] + brx * adfGeoTransform[4] +
378 2 : -2 * adfGeoTransform[5]);
379 :
380 2 : poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
381 2 : bry * adfGeoTransform[2],
382 2 : adfGeoTransform[3] + brx * adfGeoTransform[4] +
383 2 : bry * adfGeoTransform[5]);
384 :
385 2 : poUniverseRing->addPoint(adfGeoTransform[0] + -2 * adfGeoTransform[1] +
386 2 : bry * adfGeoTransform[2],
387 2 : adfGeoTransform[3] + -2 * adfGeoTransform[4] +
388 2 : bry * adfGeoTransform[5]);
389 :
390 2 : poUniverseRing->addPoint(
391 2 : adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
392 2 : adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);
393 :
394 2 : auto poUniversePoly = std::make_unique<OGRPolygon>();
395 2 : poUniversePoly->addRing(std::move(poUniverseRing));
396 2 : poInvertMP->addGeometry(std::move(poUniversePoly));
397 :
398 2 : bool bFoundNonPoly = false;
399 : // If we have GEOS, use it to "subtract" each polygon from the universe
400 : // multipolygon
401 2 : if (OGRGeometryFactory::haveGEOS())
402 : {
403 2 : OGRGeometry *poInvertMPAsGeom = poInvertMP;
404 2 : poInvertMP = nullptr;
405 2 : CPL_IGNORE_RET_VAL(poInvertMP);
406 6 : for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
407 : {
408 4 : auto poGeom = OGRGeometry::FromHandle(ahGeometries[iGeom]);
409 4 : const auto eGType = OGR_GT_Flatten(poGeom->getGeometryType());
410 4 : if (eGType != wkbPolygon && eGType != wkbMultiPolygon)
411 : {
412 1 : if (!bFoundNonPoly)
413 : {
414 1 : bFoundNonPoly = true;
415 1 : CPLError(CE_Warning, CPLE_AppDefined,
416 : "Ignoring non-polygon geometries in -i mode");
417 : }
418 : }
419 : else
420 : {
421 3 : auto poNewGeom = poInvertMPAsGeom->Difference(poGeom);
422 3 : if (poNewGeom)
423 : {
424 3 : delete poInvertMPAsGeom;
425 3 : poInvertMPAsGeom = poNewGeom;
426 : }
427 : }
428 :
429 4 : delete poGeom;
430 : }
431 :
432 2 : ahGeometries.resize(1);
433 2 : ahGeometries[0] = OGRGeometry::ToHandle(poInvertMPAsGeom);
434 2 : return;
435 : }
436 :
437 : OGRPolygon &hUniversePoly =
438 0 : *poInvertMP->getGeometryRef(poInvertMP->getNumGeometries() - 1);
439 :
440 : /* -------------------------------------------------------------------- */
441 : /* If we don't have GEOS, add outer rings of polygons as inner */
442 : /* rings of poUniversePoly and inner rings as sub-polygons. Note */
443 : /* that this only works properly if the polygons are disjoint, in */
444 : /* the sense that the outer ring of any polygon is not inside the */
445 : /* outer ring of another one. So the scenario of */
446 : /* https://github.com/OSGeo/gdal/issues/8689 with an "island" in */
447 : /* the middle of a hole will not work properly. */
448 : /* -------------------------------------------------------------------- */
449 0 : for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
450 : {
451 : const auto eGType =
452 0 : OGR_GT_Flatten(OGR_G_GetGeometryType(ahGeometries[iGeom]));
453 0 : if (eGType != wkbPolygon && eGType != wkbMultiPolygon)
454 : {
455 0 : if (!bFoundNonPoly)
456 : {
457 0 : bFoundNonPoly = true;
458 0 : CPLError(CE_Warning, CPLE_AppDefined,
459 : "Ignoring non-polygon geometries in -i mode");
460 : }
461 0 : OGR_G_DestroyGeometry(ahGeometries[iGeom]);
462 0 : continue;
463 : }
464 :
465 : const auto ProcessPoly =
466 0 : [&hUniversePoly, poInvertMP](OGRPolygon *poPoly)
467 : {
468 0 : for (int i = poPoly->getNumInteriorRings() - 1; i >= 0; --i)
469 : {
470 0 : auto poNewPoly = std::make_unique<OGRPolygon>();
471 : std::unique_ptr<OGRLinearRing> poRing(
472 0 : poPoly->stealInteriorRing(i));
473 0 : poNewPoly->addRing(std::move(poRing));
474 0 : poInvertMP->addGeometry(std::move(poNewPoly));
475 : }
476 0 : std::unique_ptr<OGRLinearRing> poShell(poPoly->stealExteriorRing());
477 0 : hUniversePoly.addRing(std::move(poShell));
478 0 : };
479 :
480 0 : if (eGType == wkbPolygon)
481 : {
482 : auto poPoly =
483 0 : OGRGeometry::FromHandle(ahGeometries[iGeom])->toPolygon();
484 0 : ProcessPoly(poPoly);
485 0 : delete poPoly;
486 : }
487 : else
488 : {
489 : auto poMulti =
490 0 : OGRGeometry::FromHandle(ahGeometries[iGeom])->toMultiPolygon();
491 0 : for (auto *poPoly : *poMulti)
492 : {
493 0 : ProcessPoly(poPoly);
494 : }
495 0 : delete poMulti;
496 : }
497 : }
498 :
499 0 : ahGeometries.resize(1);
500 0 : ahGeometries[0] = OGRGeometry::ToHandle(poInvertMP);
501 : }
502 :
503 : /************************************************************************/
504 : /* ProcessLayer() */
505 : /* */
506 : /* Process all the features in a layer selection, collecting */
507 : /* geometries and burn values. */
508 : /************************************************************************/
509 :
510 24 : static CPLErr ProcessLayer(OGRLayerH hSrcLayer, bool bSRSIsSet,
511 : GDALDatasetH hDstDS,
512 : const std::vector<int> &anBandList,
513 : const std::vector<double> &adfBurnValues, bool b3D,
514 : bool bInverse, const std::string &osBurnAttribute,
515 : CSLConstList papszRasterizeOptions,
516 : CSLConstList papszTO, GDALProgressFunc pfnProgress,
517 : void *pProgressData)
518 :
519 : {
520 : /* -------------------------------------------------------------------- */
521 : /* Checkout that SRS are the same. */
522 : /* If -a_srs is specified, skip the test */
523 : /* -------------------------------------------------------------------- */
524 24 : OGRCoordinateTransformationH hCT = nullptr;
525 24 : if (!bSRSIsSet)
526 : {
527 23 : OGRSpatialReferenceH hDstSRS = GDALGetSpatialRef(hDstDS);
528 :
529 23 : if (hDstSRS)
530 11 : hDstSRS = OSRClone(hDstSRS);
531 12 : else if (GDALGetMetadata(hDstDS, "RPC") != nullptr)
532 : {
533 2 : hDstSRS = OSRNewSpatialReference(nullptr);
534 2 : CPL_IGNORE_RET_VAL(
535 2 : OSRSetFromUserInput(hDstSRS, SRS_WKT_WGS84_LAT_LONG));
536 2 : OSRSetAxisMappingStrategy(hDstSRS, OAMS_TRADITIONAL_GIS_ORDER);
537 : }
538 :
539 23 : OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
540 23 : if (hDstSRS != nullptr && hSrcSRS != nullptr)
541 : {
542 13 : if (OSRIsSame(hSrcSRS, hDstSRS) == FALSE)
543 : {
544 1 : hCT = OCTNewCoordinateTransformation(hSrcSRS, hDstSRS);
545 1 : if (hCT == nullptr)
546 : {
547 0 : CPLError(CE_Warning, CPLE_AppDefined,
548 : "The output raster dataset and the input vector "
549 : "layer do not have the same SRS.\n"
550 : "And reprojection of input data did not work. "
551 : "Results might be incorrect.");
552 : }
553 : }
554 : }
555 10 : else if (hDstSRS != nullptr && hSrcSRS == nullptr)
556 : {
557 0 : CPLError(CE_Warning, CPLE_AppDefined,
558 : "The output raster dataset has a SRS, but the input "
559 : "vector layer SRS is unknown.\n"
560 : "Ensure input vector has the same SRS, otherwise results "
561 : "might be incorrect.");
562 : }
563 10 : else if (hDstSRS == nullptr && hSrcSRS != nullptr)
564 : {
565 1 : CPLError(CE_Warning, CPLE_AppDefined,
566 : "The input vector layer has a SRS, but the output raster "
567 : "dataset SRS is unknown.\n"
568 : "Ensure output raster dataset has the same SRS, otherwise "
569 : "results might be incorrect.");
570 : }
571 :
572 23 : if (hDstSRS != nullptr)
573 : {
574 13 : OSRDestroySpatialReference(hDstSRS);
575 : }
576 : }
577 :
578 : /* -------------------------------------------------------------------- */
579 : /* Get field index, and check. */
580 : /* -------------------------------------------------------------------- */
581 24 : int iBurnField = -1;
582 24 : bool bUseInt64 = false;
583 24 : if (!osBurnAttribute.empty())
584 : {
585 3 : OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hSrcLayer);
586 3 : iBurnField = OGR_FD_GetFieldIndex(hLayerDefn, osBurnAttribute.c_str());
587 3 : if (iBurnField == -1)
588 : {
589 0 : CPLError(CE_Failure, CPLE_AppDefined,
590 : "Failed to find field %s on layer %s.",
591 : osBurnAttribute.c_str(),
592 : OGR_FD_GetName(OGR_L_GetLayerDefn(hSrcLayer)));
593 0 : if (hCT != nullptr)
594 0 : OCTDestroyCoordinateTransformation(hCT);
595 0 : return CE_Failure;
596 : }
597 3 : if (OGR_Fld_GetType(OGR_FD_GetFieldDefn(hLayerDefn, iBurnField)) ==
598 : OFTInteger64)
599 : {
600 1 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, anBandList[0]);
601 1 : if (hBand && GDALGetRasterDataType(hBand) == GDT_Int64)
602 : {
603 1 : bUseInt64 = true;
604 : }
605 : }
606 : }
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* Collect the geometries from this layer, and build list of */
610 : /* burn values. */
611 : /* -------------------------------------------------------------------- */
612 24 : OGRFeatureH hFeat = nullptr;
613 48 : std::vector<OGRGeometryH> ahGeometries;
614 48 : std::vector<double> adfFullBurnValues;
615 24 : std::vector<int64_t> anFullBurnValues;
616 :
617 24 : OGR_L_ResetReading(hSrcLayer);
618 :
619 1926 : while ((hFeat = OGR_L_GetNextFeature(hSrcLayer)) != nullptr)
620 : {
621 1902 : OGRGeometryH hGeom = OGR_F_StealGeometry(hFeat);
622 1902 : if (hGeom == nullptr)
623 : {
624 5 : OGR_F_Destroy(hFeat);
625 5 : continue;
626 : }
627 :
628 1897 : if (hCT != nullptr)
629 : {
630 1 : if (OGR_G_Transform(hGeom, hCT) != OGRERR_NONE)
631 : {
632 0 : OGR_F_Destroy(hFeat);
633 0 : OGR_G_DestroyGeometry(hGeom);
634 0 : continue;
635 : }
636 : }
637 1897 : ahGeometries.push_back(hGeom);
638 :
639 3820 : for (unsigned int iBand = 0; iBand < anBandList.size(); iBand++)
640 : {
641 1923 : if (!adfBurnValues.empty())
642 62 : adfFullBurnValues.push_back(adfBurnValues[std::min(
643 : iBand,
644 124 : static_cast<unsigned int>(adfBurnValues.size()) - 1)]);
645 1861 : else if (!osBurnAttribute.empty())
646 : {
647 7 : if (bUseInt64)
648 1 : anFullBurnValues.push_back(
649 1 : OGR_F_GetFieldAsInteger64(hFeat, iBurnField));
650 : else
651 6 : adfFullBurnValues.push_back(
652 6 : OGR_F_GetFieldAsDouble(hFeat, iBurnField));
653 : }
654 1854 : else if (b3D)
655 : {
656 : /* Points and Lines will have their "z" values collected at the
657 : point and line levels respectively. Not implemented for
658 : polygons */
659 1854 : adfFullBurnValues.push_back(0.0);
660 : }
661 : }
662 :
663 1897 : OGR_F_Destroy(hFeat);
664 : }
665 :
666 24 : if (hCT != nullptr)
667 1 : OCTDestroyCoordinateTransformation(hCT);
668 :
669 : /* -------------------------------------------------------------------- */
670 : /* If we are in inverse mode, we add one extra ring around the */
671 : /* whole dataset to invert the concept of insideness and then */
672 : /* merge everything into one geometry collection. */
673 : /* -------------------------------------------------------------------- */
674 24 : if (bInverse)
675 : {
676 2 : if (ahGeometries.empty())
677 : {
678 0 : for (unsigned int iBand = 0; iBand < anBandList.size(); iBand++)
679 : {
680 0 : if (!adfBurnValues.empty())
681 0 : adfFullBurnValues.push_back(adfBurnValues[std::min(
682 : iBand,
683 0 : static_cast<unsigned int>(adfBurnValues.size()) - 1)]);
684 : else /* FIXME? Not sure what to do exactly in the else case, but
685 : we must insert a value */
686 : {
687 0 : adfFullBurnValues.push_back(0.0);
688 0 : anFullBurnValues.push_back(0);
689 : }
690 : }
691 : }
692 :
693 2 : InvertGeometries(hDstDS, ahGeometries);
694 : }
695 :
696 : /* -------------------------------------------------------------------- */
697 : /* If we have transformer options, create the transformer here */
698 : /* Coordinate transformation to the target SRS has already been */
699 : /* done, so we just need to convert to target raster space. */
700 : /* Note: this is somewhat identical to what is done in */
701 : /* GDALRasterizeGeometries() itself, except we can pass transformer*/
702 : /* options. */
703 : /* -------------------------------------------------------------------- */
704 :
705 24 : void *pTransformArg = nullptr;
706 24 : GDALTransformerFunc pfnTransformer = nullptr;
707 24 : CPLErr eErr = CE_None;
708 24 : if (papszTO != nullptr)
709 : {
710 1 : GDALDataset *poDS = GDALDataset::FromHandle(hDstDS);
711 1 : char **papszTransformerOptions = CSLDuplicate(papszTO);
712 1 : double adfGeoTransform[6] = {0.0};
713 1 : if (poDS->GetGeoTransform(adfGeoTransform) != CE_None &&
714 1 : poDS->GetGCPCount() == 0 && poDS->GetMetadata("RPC") == nullptr)
715 : {
716 0 : papszTransformerOptions = CSLSetNameValue(
717 : papszTransformerOptions, "DST_METHOD", "NO_GEOTRANSFORM");
718 : }
719 :
720 1 : pTransformArg = GDALCreateGenImgProjTransformer2(
721 : nullptr, hDstDS, papszTransformerOptions);
722 1 : CSLDestroy(papszTransformerOptions);
723 :
724 1 : pfnTransformer = GDALGenImgProjTransform;
725 1 : if (pTransformArg == nullptr)
726 : {
727 0 : eErr = CE_Failure;
728 : }
729 : }
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Perform the burn. */
733 : /* -------------------------------------------------------------------- */
734 24 : if (eErr == CE_None)
735 : {
736 24 : if (bUseInt64)
737 : {
738 2 : eErr = GDALRasterizeGeometriesInt64(
739 1 : hDstDS, static_cast<int>(anBandList.size()), anBandList.data(),
740 1 : static_cast<int>(ahGeometries.size()), ahGeometries.data(),
741 1 : pfnTransformer, pTransformArg, anFullBurnValues.data(),
742 : papszRasterizeOptions, pfnProgress, pProgressData);
743 : }
744 : else
745 : {
746 46 : eErr = GDALRasterizeGeometries(
747 23 : hDstDS, static_cast<int>(anBandList.size()), anBandList.data(),
748 23 : static_cast<int>(ahGeometries.size()), ahGeometries.data(),
749 23 : pfnTransformer, pTransformArg, adfFullBurnValues.data(),
750 : papszRasterizeOptions, pfnProgress, pProgressData);
751 : }
752 : }
753 :
754 : /* -------------------------------------------------------------------- */
755 : /* Cleanup */
756 : /* -------------------------------------------------------------------- */
757 :
758 24 : if (pTransformArg)
759 1 : GDALDestroyTransformer(pTransformArg);
760 :
761 1919 : for (int iGeom = static_cast<int>(ahGeometries.size()) - 1; iGeom >= 0;
762 : iGeom--)
763 1895 : OGR_G_DestroyGeometry(ahGeometries[iGeom]);
764 :
765 24 : return eErr;
766 : }
767 :
768 : /************************************************************************/
769 : /* CreateOutputDataset() */
770 : /************************************************************************/
771 :
772 12 : static GDALDatasetH CreateOutputDataset(
773 : const std::vector<OGRLayerH> &ahLayers, OGRSpatialReferenceH hSRS,
774 : OGREnvelope sEnvelop, GDALDriverH hDriver, const char *pszDest, int nXSize,
775 : int nYSize, double dfXRes, double dfYRes, bool bTargetAlignedPixels,
776 : int nBandCount, GDALDataType eOutputType, CSLConstList papszCreationOptions,
777 : const std::vector<double> &adfInitVals, const char *pszNoData)
778 : {
779 12 : bool bFirstLayer = true;
780 12 : char *pszWKT = nullptr;
781 12 : const bool bBoundsSpecifiedByUser = sEnvelop.IsInit();
782 :
783 23 : for (unsigned int i = 0; i < ahLayers.size(); i++)
784 : {
785 12 : OGRLayerH hLayer = ahLayers[i];
786 :
787 12 : if (!bBoundsSpecifiedByUser)
788 : {
789 10 : OGREnvelope sLayerEnvelop;
790 :
791 10 : if (OGR_L_GetExtent(hLayer, &sLayerEnvelop, TRUE) != OGRERR_NONE)
792 : {
793 1 : CPLError(CE_Failure, CPLE_AppDefined,
794 : "Cannot get layer extent");
795 1 : return nullptr;
796 : }
797 :
798 : /* Voluntarily increase the extent by a half-pixel size to avoid */
799 : /* missing points on the border */
800 9 : if (!bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)
801 : {
802 5 : sLayerEnvelop.MinX -= dfXRes / 2;
803 5 : sLayerEnvelop.MaxX += dfXRes / 2;
804 5 : sLayerEnvelop.MinY -= dfYRes / 2;
805 5 : sLayerEnvelop.MaxY += dfYRes / 2;
806 : }
807 :
808 9 : sEnvelop.Merge(sLayerEnvelop);
809 : }
810 :
811 11 : if (bFirstLayer)
812 : {
813 11 : if (hSRS == nullptr)
814 10 : hSRS = OGR_L_GetSpatialRef(hLayer);
815 :
816 11 : bFirstLayer = false;
817 : }
818 : }
819 :
820 11 : if (!sEnvelop.IsInit())
821 : {
822 0 : CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
823 0 : return nullptr;
824 : }
825 :
826 11 : if (dfXRes == 0 && dfYRes == 0)
827 : {
828 4 : if (nXSize == 0 || nYSize == 0)
829 : {
830 1 : CPLError(CE_Failure, CPLE_AppDefined,
831 : "Size and resolutions are missing");
832 1 : return nullptr;
833 : }
834 3 : dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;
835 3 : dfYRes = (sEnvelop.MaxY - sEnvelop.MinY) / nYSize;
836 : }
837 7 : else if (bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)
838 : {
839 0 : sEnvelop.MinX = floor(sEnvelop.MinX / dfXRes) * dfXRes;
840 0 : sEnvelop.MaxX = ceil(sEnvelop.MaxX / dfXRes) * dfXRes;
841 0 : sEnvelop.MinY = floor(sEnvelop.MinY / dfYRes) * dfYRes;
842 0 : sEnvelop.MaxY = ceil(sEnvelop.MaxY / dfYRes) * dfYRes;
843 : }
844 :
845 10 : if (dfXRes == 0 || dfYRes == 0)
846 : {
847 0 : CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
848 0 : return nullptr;
849 : }
850 :
851 10 : double adfProjection[6] = {sEnvelop.MinX, dfXRes, 0.0,
852 10 : sEnvelop.MaxY, 0.0, -dfYRes};
853 :
854 10 : if (nXSize == 0 && nYSize == 0)
855 : {
856 : // coverity[divide_by_zero]
857 7 : const double dfXSize = 0.5 + (sEnvelop.MaxX - sEnvelop.MinX) / dfXRes;
858 : // coverity[divide_by_zero]
859 7 : const double dfYSize = 0.5 + (sEnvelop.MaxY - sEnvelop.MinY) / dfYRes;
860 7 : if (dfXSize > std::numeric_limits<int>::max() ||
861 6 : dfXSize < std::numeric_limits<int>::min() ||
862 18 : dfYSize > std::numeric_limits<int>::max() ||
863 5 : dfYSize < std::numeric_limits<int>::min())
864 : {
865 2 : CPLError(CE_Failure, CPLE_AppDefined,
866 : "Invalid computed output raster size: %f x %f", dfXSize,
867 : dfYSize);
868 2 : return nullptr;
869 : }
870 5 : nXSize = static_cast<int>(dfXSize);
871 5 : nYSize = static_cast<int>(dfYSize);
872 : }
873 :
874 : GDALDatasetH hDstDS =
875 8 : GDALCreate(hDriver, pszDest, nXSize, nYSize, nBandCount, eOutputType,
876 : papszCreationOptions);
877 8 : if (hDstDS == nullptr)
878 : {
879 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszDest);
880 0 : return nullptr;
881 : }
882 :
883 8 : GDALSetGeoTransform(hDstDS, adfProjection);
884 :
885 8 : if (hSRS)
886 3 : OSRExportToWkt(hSRS, &pszWKT);
887 8 : if (pszWKT)
888 3 : GDALSetProjection(hDstDS, pszWKT);
889 8 : CPLFree(pszWKT);
890 :
891 : /*if( nBandCount == 3 || nBandCount == 4 )
892 : {
893 : for( int iBand = 0; iBand < nBandCount; iBand++ )
894 : {
895 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
896 : GDALSetRasterColorInterpretation(hBand,
897 : (GDALColorInterp)(GCI_RedBand + iBand));
898 : }
899 : }*/
900 :
901 8 : if (pszNoData)
902 : {
903 16 : for (int iBand = 0; iBand < nBandCount; iBand++)
904 : {
905 8 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
906 8 : if (GDALGetRasterDataType(hBand) == GDT_Int64)
907 1 : GDALSetRasterNoDataValueAsInt64(hBand,
908 1 : CPLAtoGIntBig(pszNoData));
909 : else
910 7 : GDALSetRasterNoDataValue(hBand, CPLAtof(pszNoData));
911 : }
912 : }
913 :
914 8 : if (!adfInitVals.empty())
915 : {
916 6 : for (int iBand = 0;
917 6 : iBand < std::min(nBandCount, static_cast<int>(adfInitVals.size()));
918 : iBand++)
919 : {
920 3 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
921 3 : GDALFillRaster(hBand, adfInitVals[iBand], 0);
922 : }
923 : }
924 :
925 8 : return hDstDS;
926 : }
927 :
928 : /************************************************************************/
929 : /* GDALRasterize() */
930 : /************************************************************************/
931 :
932 : /* clang-format off */
933 : /**
934 : * Burns vector geometries into a raster
935 : *
936 : * This is the equivalent of the
937 : * <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
938 : *
939 : * GDALRasterizeOptions* must be allocated and freed with
940 : * GDALRasterizeOptionsNew() and GDALRasterizeOptionsFree() respectively.
941 : * pszDest and hDstDS cannot be used at the same time.
942 : *
943 : * @param pszDest the destination dataset path or NULL.
944 : * @param hDstDS the destination dataset or NULL.
945 : * @param hSrcDataset the source dataset handle.
946 : * @param psOptionsIn the options struct returned by GDALRasterizeOptionsNew()
947 : * or NULL.
948 : * @param pbUsageError pointer to a integer output variable to store if any
949 : * usage error has occurred or NULL.
950 : * @return the output dataset (new dataset that must be closed using
951 : * GDALClose(), or hDstDS is not NULL) or NULL in case of error.
952 : *
953 : * @since GDAL 2.1
954 : */
955 : /* clang-format on */
956 :
957 29 : GDALDatasetH GDALRasterize(const char *pszDest, GDALDatasetH hDstDS,
958 : GDALDatasetH hSrcDataset,
959 : const GDALRasterizeOptions *psOptionsIn,
960 : int *pbUsageError)
961 : {
962 29 : if (pszDest == nullptr && hDstDS == nullptr)
963 : {
964 0 : CPLError(CE_Failure, CPLE_AppDefined,
965 : "pszDest == NULL && hDstDS == NULL");
966 :
967 0 : if (pbUsageError)
968 0 : *pbUsageError = TRUE;
969 0 : return nullptr;
970 : }
971 29 : if (hSrcDataset == nullptr)
972 : {
973 0 : CPLError(CE_Failure, CPLE_AppDefined, "hSrcDataset== NULL");
974 :
975 0 : if (pbUsageError)
976 0 : *pbUsageError = TRUE;
977 0 : return nullptr;
978 : }
979 29 : if (hDstDS != nullptr && psOptionsIn && psOptionsIn->bCreateOutput)
980 : {
981 0 : CPLError(CE_Failure, CPLE_AppDefined,
982 : "hDstDS != NULL but options that imply creating a new dataset "
983 : "have been set.");
984 :
985 0 : if (pbUsageError)
986 0 : *pbUsageError = TRUE;
987 0 : return nullptr;
988 : }
989 :
990 29 : GDALRasterizeOptions *psOptionsToFree = nullptr;
991 29 : const GDALRasterizeOptions *psOptions = psOptionsIn;
992 29 : if (psOptions == nullptr)
993 : {
994 0 : psOptionsToFree = GDALRasterizeOptionsNew(nullptr, nullptr);
995 0 : psOptions = psOptionsToFree;
996 : }
997 :
998 29 : const bool bCloseOutDSOnError = hDstDS == nullptr;
999 29 : if (pszDest == nullptr)
1000 13 : pszDest = GDALGetDescription(hDstDS);
1001 :
1002 44 : if (psOptions->osSQL.empty() && psOptions->aosLayers.empty() &&
1003 15 : GDALDatasetGetLayerCount(hSrcDataset) != 1)
1004 : {
1005 0 : CPLError(CE_Failure, CPLE_NotSupported,
1006 : "Neither -sql nor -l are specified, but the source dataset "
1007 : "has not one single layer.");
1008 0 : if (pbUsageError)
1009 0 : *pbUsageError = TRUE;
1010 0 : GDALRasterizeOptionsFree(psOptionsToFree);
1011 0 : return nullptr;
1012 : }
1013 :
1014 : /* -------------------------------------------------------------------- */
1015 : /* Open target raster file. Eventually we will add optional */
1016 : /* creation. */
1017 : /* -------------------------------------------------------------------- */
1018 29 : const bool bCreateOutput = psOptions->bCreateOutput || hDstDS == nullptr;
1019 :
1020 29 : GDALDriverH hDriver = nullptr;
1021 29 : if (bCreateOutput)
1022 : {
1023 13 : CPLString osFormat;
1024 13 : if (psOptions->osFormat.empty())
1025 : {
1026 7 : osFormat = GetOutputDriverForRaster(pszDest);
1027 7 : if (osFormat.empty())
1028 : {
1029 0 : GDALRasterizeOptionsFree(psOptionsToFree);
1030 0 : return nullptr;
1031 : }
1032 : }
1033 : else
1034 : {
1035 6 : osFormat = psOptions->osFormat;
1036 : }
1037 :
1038 : /* --------------------------------------------------------------------
1039 : */
1040 : /* Find the output driver. */
1041 : /* --------------------------------------------------------------------
1042 : */
1043 13 : hDriver = GDALGetDriverByName(osFormat);
1044 : char **papszDriverMD =
1045 13 : hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;
1046 39 : if (hDriver == nullptr ||
1047 13 : !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_RASTER,
1048 26 : "FALSE")) ||
1049 13 : !CPLTestBool(
1050 : CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))
1051 : {
1052 0 : CPLError(CE_Failure, CPLE_NotSupported,
1053 : "Output driver `%s' not recognised or does not support "
1054 : "direct output file creation.",
1055 : osFormat.c_str());
1056 0 : GDALRasterizeOptionsFree(psOptionsToFree);
1057 0 : return nullptr;
1058 : }
1059 : }
1060 :
1061 10 : const auto GetOutputDataType = [&](OGRLayerH hLayer)
1062 : {
1063 10 : CPLAssert(bCreateOutput);
1064 10 : CPLAssert(hDriver);
1065 21 : GDALDataType eOutputType = psOptions->eOutputType;
1066 10 : if (eOutputType == GDT_Unknown && !psOptions->osBurnAttribute.empty())
1067 : {
1068 1 : OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hLayer);
1069 1 : const int iBurnField = OGR_FD_GetFieldIndex(
1070 1 : hLayerDefn, psOptions->osBurnAttribute.c_str());
1071 1 : if (iBurnField >= 0 && OGR_Fld_GetType(OGR_FD_GetFieldDefn(
1072 : hLayerDefn, iBurnField)) == OFTInteger64)
1073 : {
1074 1 : const char *pszMD = GDALGetMetadataItem(
1075 : hDriver, GDAL_DMD_CREATIONDATATYPES, nullptr);
1076 2 : if (pszMD && CPLStringList(CSLTokenizeString2(pszMD, " ", 0))
1077 1 : .FindString("Int64") >= 0)
1078 : {
1079 1 : eOutputType = GDT_Int64;
1080 : }
1081 : }
1082 : }
1083 10 : if (eOutputType == GDT_Unknown)
1084 : {
1085 9 : eOutputType = GDT_Float64;
1086 : }
1087 10 : return eOutputType;
1088 29 : };
1089 :
1090 : // Store SRS handle
1091 : OGRSpatialReferenceH hSRS =
1092 29 : psOptions->oOutputSRS.IsEmpty()
1093 29 : ? nullptr
1094 1 : : OGRSpatialReference::ToHandle(
1095 1 : const_cast<OGRSpatialReference *>(&psOptions->oOutputSRS));
1096 :
1097 : /* -------------------------------------------------------------------- */
1098 : /* Process SQL request. */
1099 : /* -------------------------------------------------------------------- */
1100 29 : CPLErr eErr = CE_Failure;
1101 :
1102 29 : if (!psOptions->osSQL.empty())
1103 : {
1104 : OGRLayerH hLayer =
1105 2 : GDALDatasetExecuteSQL(hSrcDataset, psOptions->osSQL.c_str(),
1106 2 : nullptr, psOptions->osDialect.c_str());
1107 2 : if (hLayer != nullptr)
1108 : {
1109 2 : if (bCreateOutput)
1110 : {
1111 2 : std::vector<OGRLayerH> ahLayers;
1112 2 : ahLayers.push_back(hLayer);
1113 :
1114 2 : const GDALDataType eOutputType = GetOutputDataType(hLayer);
1115 2 : hDstDS = CreateOutputDataset(
1116 2 : ahLayers, hSRS, psOptions->sEnvelop, hDriver, pszDest,
1117 2 : psOptions->nXSize, psOptions->nYSize, psOptions->dfXRes,
1118 2 : psOptions->dfYRes, psOptions->bTargetAlignedPixels,
1119 2 : static_cast<int>(psOptions->anBandList.size()), eOutputType,
1120 2 : psOptions->aosCreationOptions, psOptions->adfInitVals,
1121 2 : psOptions->osNoData.c_str());
1122 2 : if (hDstDS == nullptr)
1123 : {
1124 0 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1125 0 : GDALRasterizeOptionsFree(psOptionsToFree);
1126 0 : return nullptr;
1127 : }
1128 : }
1129 :
1130 6 : eErr = ProcessLayer(
1131 2 : hLayer, hSRS != nullptr, hDstDS, psOptions->anBandList,
1132 2 : psOptions->adfBurnValues, psOptions->b3D, psOptions->bInverse,
1133 2 : psOptions->osBurnAttribute.c_str(),
1134 2 : psOptions->aosRasterizeOptions, psOptions->aosTO,
1135 2 : psOptions->pfnProgress, psOptions->pProgressData);
1136 :
1137 2 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1138 : }
1139 : }
1140 :
1141 : /* -------------------------------------------------------------------- */
1142 : /* Create output file if necessary. */
1143 : /* -------------------------------------------------------------------- */
1144 : const int nLayerCount =
1145 56 : (psOptions->osSQL.empty() && psOptions->aosLayers.empty())
1146 58 : ? 1
1147 14 : : static_cast<int>(psOptions->aosLayers.size());
1148 :
1149 29 : if (bCreateOutput && hDstDS == nullptr)
1150 : {
1151 11 : std::vector<OGRLayerH> ahLayers;
1152 :
1153 11 : GDALDataType eOutputType = psOptions->eOutputType;
1154 :
1155 21 : for (int i = 0; i < nLayerCount; i++)
1156 : {
1157 : OGRLayerH hLayer;
1158 11 : if (psOptions->aosLayers.size() > static_cast<size_t>(i))
1159 5 : hLayer = GDALDatasetGetLayerByName(
1160 5 : hSrcDataset, psOptions->aosLayers[i].c_str());
1161 : else
1162 6 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1163 11 : if (hLayer == nullptr)
1164 : {
1165 1 : CPLError(CE_Failure, CPLE_AppDefined,
1166 : "Unable to find layer \"%s\".",
1167 1 : psOptions->aosLayers.size() > static_cast<size_t>(i)
1168 1 : ? psOptions->aosLayers[i].c_str()
1169 : : "0");
1170 1 : GDALRasterizeOptionsFree(psOptionsToFree);
1171 1 : return nullptr;
1172 : }
1173 10 : if (eOutputType == GDT_Unknown)
1174 : {
1175 8 : if (GetOutputDataType(hLayer) == GDT_Int64)
1176 1 : eOutputType = GDT_Int64;
1177 : }
1178 :
1179 10 : ahLayers.push_back(hLayer);
1180 : }
1181 :
1182 10 : if (eOutputType == GDT_Unknown)
1183 : {
1184 7 : eOutputType = GDT_Float64;
1185 : }
1186 :
1187 10 : hDstDS = CreateOutputDataset(
1188 10 : ahLayers, hSRS, psOptions->sEnvelop, hDriver, pszDest,
1189 10 : psOptions->nXSize, psOptions->nYSize, psOptions->dfXRes,
1190 10 : psOptions->dfYRes, psOptions->bTargetAlignedPixels,
1191 10 : static_cast<int>(psOptions->anBandList.size()), eOutputType,
1192 10 : psOptions->aosCreationOptions, psOptions->adfInitVals,
1193 10 : psOptions->osNoData.c_str());
1194 10 : if (hDstDS == nullptr)
1195 : {
1196 4 : GDALRasterizeOptionsFree(psOptionsToFree);
1197 4 : return nullptr;
1198 : }
1199 : }
1200 :
1201 : /* -------------------------------------------------------------------- */
1202 : /* Process each layer. */
1203 : /* -------------------------------------------------------------------- */
1204 :
1205 46 : for (int i = 0; i < nLayerCount; i++)
1206 : {
1207 : OGRLayerH hLayer;
1208 22 : if (psOptions->aosLayers.size() > static_cast<size_t>(i))
1209 11 : hLayer = GDALDatasetGetLayerByName(hSrcDataset,
1210 11 : psOptions->aosLayers[i].c_str());
1211 : else
1212 11 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1213 22 : if (hLayer == nullptr)
1214 : {
1215 0 : CPLError(CE_Failure, CPLE_AppDefined,
1216 : "Unable to find layer \"%s\".",
1217 0 : psOptions->aosLayers.size() > static_cast<size_t>(i)
1218 0 : ? psOptions->aosLayers[i].c_str()
1219 : : "0");
1220 0 : eErr = CE_Failure;
1221 0 : break;
1222 : }
1223 :
1224 22 : if (!psOptions->osWHERE.empty())
1225 : {
1226 0 : if (OGR_L_SetAttributeFilter(hLayer, psOptions->osWHERE.c_str()) !=
1227 : OGRERR_NONE)
1228 : {
1229 0 : eErr = CE_Failure;
1230 0 : break;
1231 : }
1232 : }
1233 :
1234 44 : void *pScaledProgress = GDALCreateScaledProgress(
1235 22 : 0.0, 1.0 * (i + 1) / nLayerCount, psOptions->pfnProgress,
1236 22 : psOptions->pProgressData);
1237 :
1238 66 : eErr = ProcessLayer(hLayer, !psOptions->oOutputSRS.IsEmpty(), hDstDS,
1239 22 : psOptions->anBandList, psOptions->adfBurnValues,
1240 22 : psOptions->b3D, psOptions->bInverse,
1241 22 : psOptions->osBurnAttribute.c_str(),
1242 22 : psOptions->aosRasterizeOptions, psOptions->aosTO,
1243 : GDALScaledProgress, pScaledProgress);
1244 :
1245 22 : GDALDestroyScaledProgress(pScaledProgress);
1246 22 : if (eErr != CE_None)
1247 0 : break;
1248 : }
1249 :
1250 24 : GDALRasterizeOptionsFree(psOptionsToFree);
1251 :
1252 24 : if (eErr != CE_None)
1253 : {
1254 0 : if (bCloseOutDSOnError)
1255 0 : GDALClose(hDstDS);
1256 0 : return nullptr;
1257 : }
1258 :
1259 24 : return hDstDS;
1260 : }
1261 :
1262 : /************************************************************************/
1263 : /* ArgIsNumericRasterize() */
1264 : /************************************************************************/
1265 :
1266 91 : static bool ArgIsNumericRasterize(const char *pszArg)
1267 :
1268 : {
1269 91 : char *pszEnd = nullptr;
1270 91 : CPLStrtod(pszArg, &pszEnd);
1271 91 : return pszEnd != nullptr && pszEnd[0] == '\0';
1272 : }
1273 :
1274 : /************************************************************************/
1275 : /* GDALRasterizeOptionsNew() */
1276 : /************************************************************************/
1277 :
1278 : /**
1279 : * Allocates a GDALRasterizeOptions struct.
1280 : *
1281 : * @param papszArgv NULL terminated list of options (potentially including
1282 : * filename and open options too), or NULL. The accepted options are the ones of
1283 : * the <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
1284 : * @param psOptionsForBinary (output) may be NULL (and should generally be
1285 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
1286 : * GDALRasterizeOptionsForBinaryNew() prior to this
1287 : * function. Will be filled with potentially present filename, open options,...
1288 : * @return pointer to the allocated GDALRasterizeOptions struct. Must be freed
1289 : * with GDALRasterizeOptionsFree().
1290 : *
1291 : * @since GDAL 2.1
1292 : */
1293 :
1294 : GDALRasterizeOptions *
1295 29 : GDALRasterizeOptionsNew(char **papszArgv,
1296 : GDALRasterizeOptionsForBinary *psOptionsForBinary)
1297 : {
1298 :
1299 58 : auto psOptions = std::make_unique<GDALRasterizeOptions>();
1300 :
1301 : /*-------------------------------------------------------------------- */
1302 : /* Parse arguments. */
1303 : /*-------------------------------------------------------------------- */
1304 :
1305 58 : CPLStringList aosArgv;
1306 :
1307 : /* -------------------------------------------------------------------- */
1308 : /* Pre-processing for custom syntax that ArgumentParser does not */
1309 : /* support. */
1310 : /* -------------------------------------------------------------------- */
1311 29 : const int argc = CSLCount(papszArgv);
1312 222 : for (int i = 0; i < argc && papszArgv != nullptr && papszArgv[i] != nullptr;
1313 : i++)
1314 : {
1315 : // argparser will be confused if the value of a string argument
1316 : // starts with a negative sign.
1317 193 : if (EQUAL(papszArgv[i], "-a_nodata") && papszArgv[i + 1])
1318 : {
1319 4 : ++i;
1320 4 : const std::string s = papszArgv[i];
1321 4 : psOptions->osNoData = s;
1322 4 : psOptions->bCreateOutput = true;
1323 : }
1324 :
1325 : // argparser is confused by arguments that have at_least_one
1326 : // cardinality, if they immediately precede positional arguments.
1327 189 : else if (EQUAL(papszArgv[i], "-burn") && papszArgv[i + 1])
1328 : {
1329 29 : if (strchr(papszArgv[i + 1], ' '))
1330 : {
1331 : const CPLStringList aosTokens(
1332 0 : CSLTokenizeString(papszArgv[i + 1]));
1333 0 : for (const char *pszToken : aosTokens)
1334 : {
1335 0 : psOptions->adfBurnValues.push_back(CPLAtof(pszToken));
1336 : }
1337 0 : i += 1;
1338 : }
1339 : else
1340 : {
1341 58 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1342 : {
1343 29 : psOptions->adfBurnValues.push_back(
1344 29 : CPLAtof(papszArgv[i + 1]));
1345 29 : i += 1;
1346 : }
1347 : }
1348 :
1349 : // Dummy value to make argparse happy, as at least one of
1350 : // -burn, -a or -3d is required
1351 29 : aosArgv.AddString("-burn");
1352 29 : aosArgv.AddString("0");
1353 : }
1354 160 : else if (EQUAL(papszArgv[i], "-init") && papszArgv[i + 1])
1355 : {
1356 3 : if (strchr(papszArgv[i + 1], ' '))
1357 : {
1358 : const CPLStringList aosTokens(
1359 0 : CSLTokenizeString(papszArgv[i + 1]));
1360 0 : for (const char *pszToken : aosTokens)
1361 : {
1362 0 : psOptions->adfInitVals.push_back(CPLAtof(pszToken));
1363 : }
1364 0 : i += 1;
1365 : }
1366 : else
1367 : {
1368 6 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1369 : {
1370 3 : psOptions->adfInitVals.push_back(CPLAtof(papszArgv[i + 1]));
1371 3 : i += 1;
1372 : }
1373 : }
1374 3 : psOptions->bCreateOutput = true;
1375 : }
1376 157 : else if (EQUAL(papszArgv[i], "-b") && papszArgv[i + 1])
1377 : {
1378 18 : if (strchr(papszArgv[i + 1], ' '))
1379 : {
1380 : const CPLStringList aosTokens(
1381 0 : CSLTokenizeString(papszArgv[i + 1]));
1382 0 : for (const char *pszToken : aosTokens)
1383 : {
1384 0 : psOptions->anBandList.push_back(atoi(pszToken));
1385 : }
1386 0 : i += 1;
1387 : }
1388 : else
1389 : {
1390 36 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1391 : {
1392 18 : psOptions->anBandList.push_back(atoi(papszArgv[i + 1]));
1393 18 : i += 1;
1394 : }
1395 18 : }
1396 : }
1397 : else
1398 : {
1399 139 : aosArgv.AddString(papszArgv[i]);
1400 : }
1401 : }
1402 :
1403 : try
1404 : {
1405 : auto argParser =
1406 29 : GDALRasterizeOptionsGetParser(psOptions.get(), psOptionsForBinary);
1407 29 : argParser->parse_args_without_binary_name(aosArgv.List());
1408 :
1409 : // Check all no store_into args
1410 31 : if (auto oTe = argParser->present<std::vector<double>>("-te"))
1411 : {
1412 2 : psOptions->sEnvelop.MinX = oTe.value()[0];
1413 2 : psOptions->sEnvelop.MinY = oTe.value()[1];
1414 2 : psOptions->sEnvelop.MaxX = oTe.value()[2];
1415 2 : psOptions->sEnvelop.MaxY = oTe.value()[3];
1416 2 : psOptions->bCreateOutput = true;
1417 : }
1418 :
1419 29 : if (auto oTr = argParser->present<std::vector<double>>("-tr"))
1420 : {
1421 7 : psOptions->dfXRes = oTr.value()[0];
1422 7 : psOptions->dfYRes = oTr.value()[1];
1423 :
1424 7 : if (psOptions->dfXRes <= 0 || psOptions->dfYRes <= 0)
1425 : {
1426 0 : CPLError(CE_Failure, CPLE_AppDefined,
1427 : "Wrong value for -tr parameter.");
1428 0 : return nullptr;
1429 : }
1430 :
1431 7 : psOptions->bCreateOutput = true;
1432 : }
1433 :
1434 29 : if (auto oTs = argParser->present<std::vector<int>>("-ts"))
1435 : {
1436 5 : psOptions->nXSize = oTs.value()[0];
1437 5 : psOptions->nYSize = oTs.value()[1];
1438 :
1439 5 : if (psOptions->nXSize <= 0 || psOptions->nYSize <= 0)
1440 : {
1441 0 : CPLError(CE_Failure, CPLE_AppDefined,
1442 : "Wrong value for -ts parameter.");
1443 0 : return nullptr;
1444 : }
1445 :
1446 5 : psOptions->bCreateOutput = true;
1447 : }
1448 :
1449 29 : if (psOptions->bCreateOutput)
1450 : {
1451 17 : if (psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
1452 17 : psOptions->nXSize == 0 && psOptions->nYSize == 0)
1453 : {
1454 0 : CPLError(CE_Failure, CPLE_NotSupported,
1455 : "'-tr xres yres' or '-ts xsize ysize' is required.");
1456 0 : return nullptr;
1457 : }
1458 :
1459 12 : if (psOptions->bTargetAlignedPixels && psOptions->dfXRes == 0 &&
1460 0 : psOptions->dfYRes == 0)
1461 : {
1462 0 : CPLError(CE_Failure, CPLE_NotSupported,
1463 : "-tap option cannot be used without using -tr.");
1464 0 : return nullptr;
1465 : }
1466 :
1467 12 : if (!psOptions->anBandList.empty())
1468 : {
1469 0 : CPLError(
1470 : CE_Failure, CPLE_NotSupported,
1471 : "-b option cannot be used when creating a GDAL dataset.");
1472 0 : return nullptr;
1473 : }
1474 :
1475 12 : int nBandCount = 1;
1476 :
1477 12 : if (!psOptions->adfBurnValues.empty())
1478 3 : nBandCount = static_cast<int>(psOptions->adfBurnValues.size());
1479 :
1480 12 : if (static_cast<int>(psOptions->adfInitVals.size()) > nBandCount)
1481 0 : nBandCount = static_cast<int>(psOptions->adfInitVals.size());
1482 :
1483 12 : if (psOptions->adfInitVals.size() == 1)
1484 : {
1485 3 : for (int i = 1; i <= nBandCount - 1; i++)
1486 0 : psOptions->adfInitVals.push_back(psOptions->adfInitVals[0]);
1487 : }
1488 :
1489 24 : for (int i = 1; i <= nBandCount; i++)
1490 12 : psOptions->anBandList.push_back(i);
1491 : }
1492 : else
1493 : {
1494 17 : if (psOptions->anBandList.empty())
1495 11 : psOptions->anBandList.push_back(1);
1496 : }
1497 :
1498 29 : if (!psOptions->osDialect.empty() && !psOptions->osWHERE.empty() &&
1499 0 : !psOptions->osSQL.empty())
1500 : {
1501 0 : CPLError(CE_Warning, CPLE_AppDefined,
1502 : "-dialect is ignored with -where. Use -sql instead");
1503 : }
1504 :
1505 29 : if (psOptionsForBinary)
1506 : {
1507 9 : psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
1508 9 : if (!psOptions->osFormat.empty())
1509 0 : psOptionsForBinary->osFormat = psOptions->osFormat;
1510 : }
1511 28 : else if (psOptions->adfBurnValues.empty() &&
1512 28 : psOptions->osBurnAttribute.empty() && !psOptions->b3D)
1513 : {
1514 6 : psOptions->adfBurnValues.push_back(255);
1515 : }
1516 : }
1517 0 : catch (const std::exception &e)
1518 : {
1519 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
1520 0 : return nullptr;
1521 : }
1522 :
1523 29 : return psOptions.release();
1524 : }
1525 :
1526 : /************************************************************************/
1527 : /* GDALRasterizeOptionsFree() */
1528 : /************************************************************************/
1529 :
1530 : /**
1531 : * Frees the GDALRasterizeOptions struct.
1532 : *
1533 : * @param psOptions the options struct for GDALRasterize().
1534 : *
1535 : * @since GDAL 2.1
1536 : */
1537 :
1538 58 : void GDALRasterizeOptionsFree(GDALRasterizeOptions *psOptions)
1539 : {
1540 58 : delete psOptions;
1541 58 : }
1542 :
1543 : /************************************************************************/
1544 : /* GDALRasterizeOptionsSetProgress() */
1545 : /************************************************************************/
1546 :
1547 : /**
1548 : * Set a progress function.
1549 : *
1550 : * @param psOptions the options struct for GDALRasterize().
1551 : * @param pfnProgress the progress callback.
1552 : * @param pProgressData the user data for the progress callback.
1553 : *
1554 : * @since GDAL 2.1
1555 : */
1556 :
1557 8 : void GDALRasterizeOptionsSetProgress(GDALRasterizeOptions *psOptions,
1558 : GDALProgressFunc pfnProgress,
1559 : void *pProgressData)
1560 : {
1561 8 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
1562 8 : psOptions->pProgressData = pProgressData;
1563 8 : }
|
