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