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 62 : GDALRasterizeOptionsGetParser(GDALRasterizeOptions *psOptions,
80 : GDALRasterizeOptionsForBinary *psOptionsForBinary)
81 : {
82 : auto argParser = std::make_unique<GDALArgumentParser>(
83 62 : "gdal_rasterize", /* bForBinary=*/psOptionsForBinary != nullptr);
84 :
85 62 : argParser->add_description(_("Burns vector geometries into a raster."));
86 :
87 62 : argParser->add_epilog(
88 : _("This program burns vector geometries (points, lines, and polygons) "
89 62 : "into the raster band(s) of a raster image."));
90 :
91 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
92 62 : argParser->add_argument("-b")
93 124 : .metavar("<band>")
94 62 : .append()
95 62 : .scan<'i', int>()
96 : //.nargs(argparse::nargs_pattern::at_least_one)
97 62 : .help(_("The band(s) to burn values into."));
98 :
99 62 : argParser->add_argument("-i")
100 62 : .flag()
101 62 : .store_into(psOptions->bInverse)
102 62 : .help(_("Invert rasterization."));
103 :
104 62 : argParser->add_argument("-at")
105 62 : .flag()
106 : .action(
107 28 : [psOptions](const std::string &) {
108 : psOptions->aosRasterizeOptions.SetNameValue("ALL_TOUCHED",
109 28 : "TRUE");
110 62 : })
111 62 : .help(_("Enables the ALL_TOUCHED rasterization option."));
112 :
113 : // Mutually exclusive options: -burn, -3d, -a
114 : {
115 : // Required if options for binary
116 62 : auto &group = argParser->add_mutually_exclusive_group(
117 62 : psOptionsForBinary != nullptr);
118 :
119 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
120 62 : group.add_argument("-burn")
121 124 : .metavar("<value>")
122 62 : .scan<'g', double>()
123 62 : .append()
124 : //.nargs(argparse::nargs_pattern::at_least_one)
125 62 : .help(_("A fixed value to burn into the raster band(s)."));
126 :
127 62 : group.add_argument("-a")
128 124 : .metavar("<attribute_name>")
129 62 : .store_into(psOptions->osBurnAttribute)
130 : .help(_("Name of the field in the input layer to get the burn "
131 62 : "values from."));
132 :
133 62 : group.add_argument("-3d")
134 62 : .flag()
135 62 : .store_into(psOptions->b3D)
136 : .action(
137 5 : [psOptions](const std::string &) {
138 : psOptions->aosRasterizeOptions.SetNameValue(
139 5 : "BURN_VALUE_FROM", "Z");
140 62 : })
141 : .help(_("Indicates that a burn value should be extracted from the "
142 62 : "\"Z\" values of the feature."));
143 : }
144 :
145 62 : argParser->add_argument("-add")
146 62 : .flag()
147 : .action(
148 1 : [psOptions](const std::string &) {
149 1 : psOptions->aosRasterizeOptions.SetNameValue("MERGE_ALG", "ADD");
150 62 : })
151 : .help(_("Instead of burning a new value, this adds the new value to "
152 62 : "the existing raster."));
153 :
154 : // Undocumented
155 62 : argParser->add_argument("-chunkysize")
156 62 : .flag()
157 62 : .hidden()
158 : .action(
159 0 : [psOptions](const std::string &s) {
160 : psOptions->aosRasterizeOptions.SetNameValue("CHUNKYSIZE",
161 0 : s.c_str());
162 62 : });
163 :
164 : // Mutually exclusive -l, -sql
165 : {
166 62 : auto &group = argParser->add_mutually_exclusive_group(false);
167 :
168 62 : group.add_argument("-l")
169 124 : .metavar("<layer_name>")
170 62 : .append()
171 62 : .store_into(psOptions->aosLayers)
172 62 : .help(_("Name of the layer(s) to process."));
173 :
174 62 : group.add_argument("-sql")
175 124 : .metavar("<sql_statement>")
176 62 : .store_into(psOptions->osSQL)
177 : .action(
178 7 : [psOptions](const std::string &sql)
179 : {
180 6 : GByte *pabyRet = nullptr;
181 7 : if (!sql.empty() && sql.at(0) == '@' &&
182 7 : 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 : psOptions->osSQL =
189 1 : CPLRemoveSQLComments(pszSQLStatement);
190 1 : VSIFree(pszSQLStatement);
191 : }
192 68 : })
193 : .help(
194 : _("An SQL statement to be evaluated against the datasource to "
195 62 : "produce a virtual layer of features to be burned in."));
196 : }
197 :
198 62 : argParser->add_argument("-where")
199 124 : .metavar("<expression>")
200 62 : .store_into(psOptions->osWHERE)
201 : .help(_("An optional SQL WHERE style query expression to be applied to "
202 : "select features "
203 62 : "to burn in from the input layer(s)."));
204 :
205 62 : argParser->add_argument("-dialect")
206 124 : .metavar("<sql_dialect>")
207 62 : .store_into(psOptions->osDialect)
208 62 : .help(_("The SQL dialect to use for the SQL expression."));
209 :
210 : // Store later
211 62 : argParser->add_argument("-a_nodata")
212 124 : .metavar("<value>")
213 62 : .help(_("Assign a specified nodata value to output bands."));
214 :
215 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
216 62 : argParser->add_argument("-init")
217 124 : .metavar("<value>")
218 62 : .append()
219 : //.nargs(argparse::nargs_pattern::at_least_one)
220 62 : .scan<'g', double>()
221 62 : .help(_("Initialize the output bands to the specified value."));
222 :
223 62 : argParser->add_argument("-a_srs")
224 124 : .metavar("<srs_def>")
225 : .action(
226 4 : [psOptions](const std::string &osOutputSRSDef)
227 : {
228 2 : if (psOptions->oOutputSRS.SetFromUserInput(
229 2 : 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 2 : psOptions->bCreateOutput = true;
236 64 : })
237 62 : .help(_("The spatial reference system to use for the output raster."));
238 :
239 62 : argParser->add_argument("-to")
240 124 : .metavar("<NAME>=<VALUE>")
241 62 : .append()
242 1 : .action([psOptions](const std::string &s)
243 63 : { psOptions->aosTO.AddString(s.c_str()); })
244 62 : .help(_("Set a transformer option."));
245 :
246 : // Store later
247 62 : argParser->add_argument("-te")
248 124 : .metavar("<xmin> <ymin> <xmax> <ymax>")
249 62 : .nargs(4)
250 62 : .scan<'g', double>()
251 62 : .help(_("Set georeferenced extents of output file to be created."));
252 :
253 : // Mutex with tr
254 : {
255 62 : auto &group = argParser->add_mutually_exclusive_group(false);
256 :
257 : // Store later
258 62 : group.add_argument("-tr")
259 124 : .metavar("<xres> <yres>")
260 62 : .nargs(2)
261 62 : .scan<'g', double>()
262 : .help(
263 62 : _("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 62 : auto &arg = group.add_argument("-ts")
269 124 : .metavar("<width> <height>")
270 62 : .nargs(2)
271 62 : .scan<'g', double>()
272 62 : .help(_("Set output file size in pixels and lines."));
273 :
274 62 : argParser->add_hidden_alias_for(arg, "-outsize");
275 : }
276 :
277 62 : argParser->add_argument("-tap")
278 62 : .flag()
279 62 : .store_into(psOptions->bTargetAlignedPixels)
280 1 : .action([psOptions](const std::string &)
281 62 : { psOptions->bCreateOutput = true; })
282 : .help(_("Align the coordinates of the extent to the values of the "
283 62 : "output raster."));
284 :
285 62 : argParser->add_argument("-optim")
286 124 : .metavar("AUTO|VECTOR|RASTER")
287 : .action(
288 31 : [psOptions](const std::string &s) {
289 31 : psOptions->aosRasterizeOptions.SetNameValue("OPTIM", s.c_str());
290 62 : })
291 62 : .help(_("Force the algorithm used."));
292 :
293 62 : argParser->add_creation_options_argument(psOptions->aosCreationOptions)
294 2 : .action([psOptions](const std::string &)
295 62 : { psOptions->bCreateOutput = true; });
296 :
297 62 : argParser->add_output_type_argument(psOptions->eOutputType)
298 2 : .action([psOptions](const std::string &)
299 62 : { psOptions->bCreateOutput = true; });
300 :
301 62 : argParser->add_output_format_argument(psOptions->osFormat)
302 11 : .action([psOptions](const std::string &)
303 62 : { 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 62 : psOptionsForBinary ? &(psOptionsForBinary->bQuiet) : nullptr);
309 :
310 62 : 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 62 : 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 3 : static void InvertGeometries(GDALDatasetH hDstDS,
357 : std::vector<OGRGeometryH> &ahGeometries)
358 :
359 : {
360 3 : OGRMultiPolygon *poInvertMP = new OGRMultiPolygon();
361 :
362 : /* -------------------------------------------------------------------- */
363 : /* Create a ring that is a bit outside the raster dataset. */
364 : /* -------------------------------------------------------------------- */
365 3 : const int brx = GDALGetRasterXSize(hDstDS) + 2;
366 3 : const int bry = GDALGetRasterYSize(hDstDS) + 2;
367 :
368 3 : double adfGeoTransform[6] = {};
369 3 : GDALGetGeoTransform(hDstDS, adfGeoTransform);
370 :
371 3 : auto poUniverseRing = std::make_unique<OGRLinearRing>();
372 :
373 3 : poUniverseRing->addPoint(
374 3 : adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
375 3 : adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);
376 :
377 3 : poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
378 3 : -2 * adfGeoTransform[2],
379 3 : adfGeoTransform[3] + brx * adfGeoTransform[4] +
380 3 : -2 * adfGeoTransform[5]);
381 :
382 3 : poUniverseRing->addPoint(adfGeoTransform[0] + brx * adfGeoTransform[1] +
383 3 : bry * adfGeoTransform[2],
384 3 : adfGeoTransform[3] + brx * adfGeoTransform[4] +
385 3 : bry * adfGeoTransform[5]);
386 :
387 3 : poUniverseRing->addPoint(adfGeoTransform[0] + -2 * adfGeoTransform[1] +
388 3 : bry * adfGeoTransform[2],
389 3 : adfGeoTransform[3] + -2 * adfGeoTransform[4] +
390 3 : bry * adfGeoTransform[5]);
391 :
392 3 : poUniverseRing->addPoint(
393 3 : adfGeoTransform[0] + -2 * adfGeoTransform[1] + -2 * adfGeoTransform[2],
394 3 : adfGeoTransform[3] + -2 * adfGeoTransform[4] + -2 * adfGeoTransform[5]);
395 :
396 3 : auto poUniversePoly = std::make_unique<OGRPolygon>();
397 3 : poUniversePoly->addRing(std::move(poUniverseRing));
398 3 : poInvertMP->addGeometry(std::move(poUniversePoly));
399 :
400 3 : bool bFoundNonPoly = false;
401 : // If we have GEOS, use it to "subtract" each polygon from the universe
402 : // multipolygon
403 3 : if (OGRGeometryFactory::haveGEOS())
404 : {
405 3 : OGRGeometry *poInvertMPAsGeom = poInvertMP;
406 3 : poInvertMP = nullptr;
407 3 : CPL_IGNORE_RET_VAL(poInvertMP);
408 10 : for (unsigned int iGeom = 0; iGeom < ahGeometries.size(); iGeom++)
409 : {
410 7 : auto poGeom = OGRGeometry::FromHandle(ahGeometries[iGeom]);
411 7 : const auto eGType = OGR_GT_Flatten(poGeom->getGeometryType());
412 7 : 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 6 : auto poNewGeom = poInvertMPAsGeom->Difference(poGeom);
424 6 : if (poNewGeom)
425 : {
426 6 : delete poInvertMPAsGeom;
427 6 : poInvertMPAsGeom = poNewGeom;
428 : }
429 : }
430 :
431 7 : delete poGeom;
432 : }
433 :
434 3 : ahGeometries.resize(1);
435 3 : ahGeometries[0] = OGRGeometry::ToHandle(poInvertMPAsGeom);
436 3 : 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 49 : 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 49 : OGRCoordinateTransformationH hCT = nullptr;
527 49 : if (!bSRSIsSet)
528 : {
529 47 : OGRSpatialReferenceH hDstSRS = GDALGetSpatialRef(hDstDS);
530 :
531 47 : if (hDstSRS)
532 13 : hDstSRS = OSRClone(hDstSRS);
533 34 : 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 47 : OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
542 47 : if (hDstSRS != nullptr && hSrcSRS != nullptr)
543 : {
544 15 : 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 32 : 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 32 : 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 47 : if (hDstSRS != nullptr)
575 : {
576 15 : OSRDestroySpatialReference(hDstSRS);
577 : }
578 : }
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* Get field index, and check. */
582 : /* -------------------------------------------------------------------- */
583 49 : int iBurnField = -1;
584 49 : bool bUseInt64 = false;
585 49 : if (!osBurnAttribute.empty())
586 : {
587 5 : OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hSrcLayer);
588 5 : iBurnField = OGR_FD_GetFieldIndex(hLayerDefn, osBurnAttribute.c_str());
589 5 : if (iBurnField == -1)
590 : {
591 1 : 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 1 : if (hCT != nullptr)
596 0 : OCTDestroyCoordinateTransformation(hCT);
597 1 : return CE_Failure;
598 : }
599 4 : 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 48 : OGRFeatureH hFeat = nullptr;
615 96 : std::vector<OGRGeometryH> ahGeometries;
616 96 : std::vector<double> adfFullBurnValues;
617 48 : std::vector<int64_t> anFullBurnValues;
618 :
619 48 : OGR_L_ResetReading(hSrcLayer);
620 :
621 2026 : while ((hFeat = OGR_L_GetNextFeature(hSrcLayer)) != nullptr)
622 : {
623 1978 : OGRGeometryH hGeom = OGR_F_StealGeometry(hFeat);
624 1978 : if (hGeom == nullptr)
625 : {
626 5 : OGR_F_Destroy(hFeat);
627 5 : continue;
628 : }
629 :
630 1973 : 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 1973 : ahGeometries.push_back(hGeom);
640 :
641 4072 : for (unsigned int iBand = 0; iBand < anBandList.size(); iBand++)
642 : {
643 2099 : if (!adfBurnValues.empty())
644 220 : adfFullBurnValues.push_back(adfBurnValues[std::min(
645 : iBand,
646 440 : static_cast<unsigned int>(adfBurnValues.size()) - 1)]);
647 1879 : else if (!osBurnAttribute.empty())
648 : {
649 16 : if (bUseInt64)
650 1 : anFullBurnValues.push_back(
651 1 : OGR_F_GetFieldAsInteger64(hFeat, iBurnField));
652 : else
653 15 : adfFullBurnValues.push_back(
654 15 : OGR_F_GetFieldAsDouble(hFeat, iBurnField));
655 : }
656 1863 : 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 1863 : adfFullBurnValues.push_back(0.0);
662 : }
663 : }
664 :
665 1973 : OGR_F_Destroy(hFeat);
666 : }
667 :
668 48 : 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 48 : if (bInverse)
677 : {
678 3 : 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 3 : 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 48 : void *pTransformArg = nullptr;
708 48 : GDALTransformerFunc pfnTransformer = nullptr;
709 48 : CPLErr eErr = CE_None;
710 48 : 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 48 : if (eErr == CE_None)
737 : {
738 48 : 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 94 : eErr = GDALRasterizeGeometries(
749 47 : hDstDS, static_cast<int>(anBandList.size()), anBandList.data(),
750 47 : static_cast<int>(ahGeometries.size()), ahGeometries.data(),
751 47 : pfnTransformer, pTransformArg, adfFullBurnValues.data(),
752 : papszRasterizeOptions, pfnProgress, pProgressData);
753 : }
754 : }
755 :
756 : /* -------------------------------------------------------------------- */
757 : /* Cleanup */
758 : /* -------------------------------------------------------------------- */
759 :
760 48 : if (pTransformArg)
761 1 : GDALDestroyTransformer(pTransformArg);
762 :
763 2017 : for (int iGeom = static_cast<int>(ahGeometries.size()) - 1; iGeom >= 0;
764 : iGeom--)
765 1969 : OGR_G_DestroyGeometry(ahGeometries[iGeom]);
766 :
767 48 : return eErr;
768 : }
769 :
770 : /************************************************************************/
771 : /* CreateOutputDataset() */
772 : /************************************************************************/
773 :
774 25 : 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 25 : bool bFirstLayer = true;
782 25 : char *pszWKT = nullptr;
783 25 : const bool bBoundsSpecifiedByUser = sEnvelop.IsInit();
784 :
785 49 : for (unsigned int i = 0; i < ahLayers.size(); i++)
786 : {
787 25 : OGRLayerH hLayer = ahLayers[i];
788 :
789 25 : if (!bBoundsSpecifiedByUser)
790 : {
791 22 : OGREnvelope sLayerEnvelop;
792 :
793 22 : 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 21 : if (!bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)
803 : {
804 9 : sLayerEnvelop.MinX -= dfXRes / 2;
805 9 : sLayerEnvelop.MaxX += dfXRes / 2;
806 9 : sLayerEnvelop.MinY -= dfYRes / 2;
807 9 : sLayerEnvelop.MaxY += dfYRes / 2;
808 : }
809 :
810 21 : sEnvelop.Merge(sLayerEnvelop);
811 : }
812 :
813 24 : if (bFirstLayer)
814 : {
815 24 : if (hSRS == nullptr)
816 22 : hSRS = OGR_L_GetSpatialRef(hLayer);
817 :
818 24 : bFirstLayer = false;
819 : }
820 : }
821 :
822 24 : if (!sEnvelop.IsInit())
823 : {
824 0 : CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
825 0 : return nullptr;
826 : }
827 :
828 24 : if (dfXRes == 0 && dfYRes == 0)
829 : {
830 11 : if (nXSize == 0 || nYSize == 0)
831 : {
832 2 : CPLError(CE_Failure, CPLE_AppDefined,
833 : "Size and resolution are missing");
834 2 : return nullptr;
835 : }
836 9 : dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;
837 9 : dfYRes = (sEnvelop.MaxY - sEnvelop.MinY) / nYSize;
838 : }
839 13 : else if (bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)
840 : {
841 1 : sEnvelop.MinX = floor(sEnvelop.MinX / dfXRes) * dfXRes;
842 1 : sEnvelop.MaxX = ceil(sEnvelop.MaxX / dfXRes) * dfXRes;
843 1 : sEnvelop.MinY = floor(sEnvelop.MinY / dfYRes) * dfYRes;
844 1 : sEnvelop.MaxY = ceil(sEnvelop.MaxY / dfYRes) * dfYRes;
845 : }
846 :
847 22 : if (dfXRes == 0 || dfYRes == 0)
848 : {
849 0 : CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
850 0 : return nullptr;
851 : }
852 :
853 22 : double adfProjection[6] = {sEnvelop.MinX, dfXRes, 0.0,
854 22 : sEnvelop.MaxY, 0.0, -dfYRes};
855 :
856 22 : if (nXSize == 0 && nYSize == 0)
857 : {
858 : // coverity[divide_by_zero]
859 13 : const double dfXSize = 0.5 + (sEnvelop.MaxX - sEnvelop.MinX) / dfXRes;
860 : // coverity[divide_by_zero]
861 13 : const double dfYSize = 0.5 + (sEnvelop.MaxY - sEnvelop.MinY) / dfYRes;
862 13 : if (dfXSize > std::numeric_limits<int>::max() ||
863 12 : dfXSize < std::numeric_limits<int>::min() ||
864 36 : dfYSize > std::numeric_limits<int>::max() ||
865 11 : 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 11 : nXSize = static_cast<int>(dfXSize);
873 11 : nYSize = static_cast<int>(dfYSize);
874 : }
875 :
876 : GDALDatasetH hDstDS =
877 20 : GDALCreate(hDriver, pszDest, nXSize, nYSize, nBandCount, eOutputType,
878 : papszCreationOptions);
879 20 : if (hDstDS == nullptr)
880 : {
881 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszDest);
882 0 : return nullptr;
883 : }
884 :
885 20 : GDALSetGeoTransform(hDstDS, adfProjection);
886 :
887 20 : if (hSRS)
888 6 : OSRExportToWkt(hSRS, &pszWKT);
889 20 : if (pszWKT)
890 6 : GDALSetProjection(hDstDS, pszWKT);
891 20 : 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 20 : if (pszNoData)
904 : {
905 56 : for (int iBand = 0; iBand < nBandCount; iBand++)
906 : {
907 36 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
908 36 : if (GDALGetRasterDataType(hBand) == GDT_Int64)
909 1 : GDALSetRasterNoDataValueAsInt64(hBand,
910 1 : CPLAtoGIntBig(pszNoData));
911 : else
912 35 : GDALSetRasterNoDataValue(hBand, CPLAtof(pszNoData));
913 : }
914 : }
915 :
916 20 : if (!adfInitVals.empty())
917 : {
918 10 : for (int iBand = 0;
919 10 : iBand < std::min(nBandCount, static_cast<int>(adfInitVals.size()));
920 : iBand++)
921 : {
922 6 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
923 6 : GDALFillRaster(hBand, adfInitVals[iBand], 0);
924 : }
925 : }
926 :
927 20 : 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 58 : GDALDatasetH GDALRasterize(const char *pszDest, GDALDatasetH hDstDS,
960 : GDALDatasetH hSrcDataset,
961 : const GDALRasterizeOptions *psOptionsIn,
962 : int *pbUsageError)
963 : {
964 58 : 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 58 : 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 58 : 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 116 : psOptionsToFree(nullptr, GDALRasterizeOptionsFree);
994 58 : const GDALRasterizeOptions *psOptions = psOptionsIn;
995 58 : if (psOptions == nullptr)
996 : {
997 0 : psOptionsToFree.reset(GDALRasterizeOptionsNew(nullptr, nullptr));
998 0 : psOptions = psOptionsToFree.get();
999 : }
1000 :
1001 58 : const bool bCloseOutDSOnError = hDstDS == nullptr;
1002 58 : if (pszDest == nullptr)
1003 13 : pszDest = GDALGetDescription(hDstDS);
1004 :
1005 83 : if (psOptions->osSQL.empty() && psOptions->aosLayers.empty() &&
1006 25 : 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 58 : const bool bCreateOutput = psOptions->bCreateOutput || hDstDS == nullptr;
1021 :
1022 58 : GDALDriverH hDriver = nullptr;
1023 58 : if (bCreateOutput)
1024 : {
1025 29 : CPLString osFormat;
1026 29 : if (psOptions->osFormat.empty())
1027 : {
1028 18 : osFormat = GetOutputDriverForRaster(pszDest);
1029 18 : if (osFormat.empty())
1030 : {
1031 0 : return nullptr;
1032 : }
1033 : }
1034 : else
1035 : {
1036 11 : osFormat = psOptions->osFormat;
1037 : }
1038 :
1039 : /* --------------------------------------------------------------------
1040 : */
1041 : /* Find the output driver. */
1042 : /* --------------------------------------------------------------------
1043 : */
1044 29 : hDriver = GDALGetDriverByName(osFormat);
1045 : char **papszDriverMD =
1046 29 : hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;
1047 29 : if (hDriver == nullptr)
1048 : {
1049 1 : CPLError(CE_Failure, CPLE_NotSupported,
1050 : "Output driver `%s' not recognised.", osFormat.c_str());
1051 1 : return nullptr;
1052 : }
1053 28 : if (!CPLTestBool(
1054 : CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_RASTER, "FALSE")))
1055 : {
1056 1 : CPLError(CE_Failure, CPLE_NotSupported,
1057 : "Output driver `%s' is not a raster driver.",
1058 : osFormat.c_str());
1059 1 : return nullptr;
1060 : }
1061 27 : if (!CPLTestBool(
1062 : CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))
1063 : {
1064 1 : CPLError(CE_Failure, CPLE_NotSupported,
1065 : "Output driver `%s' does not support direct output file "
1066 : "creation. "
1067 : "To write a file to this format, first write to a "
1068 : "different format such as "
1069 : "GeoTIFF and then convert the output.",
1070 : osFormat.c_str());
1071 1 : return nullptr;
1072 : }
1073 : }
1074 :
1075 23 : const auto GetOutputDataType = [&](OGRLayerH hLayer)
1076 : {
1077 23 : CPLAssert(bCreateOutput);
1078 23 : CPLAssert(hDriver);
1079 47 : GDALDataType eOutputType = psOptions->eOutputType;
1080 23 : if (eOutputType == GDT_Unknown && !psOptions->osBurnAttribute.empty())
1081 : {
1082 1 : OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hLayer);
1083 1 : const int iBurnField = OGR_FD_GetFieldIndex(
1084 1 : hLayerDefn, psOptions->osBurnAttribute.c_str());
1085 1 : if (iBurnField >= 0 && OGR_Fld_GetType(OGR_FD_GetFieldDefn(
1086 : hLayerDefn, iBurnField)) == OFTInteger64)
1087 : {
1088 1 : const char *pszMD = GDALGetMetadataItem(
1089 : hDriver, GDAL_DMD_CREATIONDATATYPES, nullptr);
1090 2 : if (pszMD && CPLStringList(CSLTokenizeString2(pszMD, " ", 0))
1091 1 : .FindString("Int64") >= 0)
1092 : {
1093 1 : eOutputType = GDT_Int64;
1094 : }
1095 : }
1096 : }
1097 23 : if (eOutputType == GDT_Unknown)
1098 : {
1099 22 : eOutputType = GDT_Float64;
1100 : }
1101 23 : return eOutputType;
1102 55 : };
1103 :
1104 : // Store SRS handle
1105 : OGRSpatialReferenceH hSRS =
1106 55 : psOptions->oOutputSRS.IsEmpty()
1107 55 : ? nullptr
1108 2 : : OGRSpatialReference::ToHandle(
1109 2 : const_cast<OGRSpatialReference *>(&psOptions->oOutputSRS));
1110 :
1111 : /* -------------------------------------------------------------------- */
1112 : /* Process SQL request. */
1113 : /* -------------------------------------------------------------------- */
1114 55 : CPLErr eErr = CE_Failure;
1115 :
1116 55 : if (!psOptions->osSQL.empty())
1117 : {
1118 : OGRLayerH hLayer =
1119 6 : GDALDatasetExecuteSQL(hSrcDataset, psOptions->osSQL.c_str(),
1120 6 : nullptr, psOptions->osDialect.c_str());
1121 6 : if (hLayer != nullptr)
1122 : {
1123 6 : if (bCreateOutput)
1124 : {
1125 2 : std::vector<OGRLayerH> ahLayers;
1126 2 : ahLayers.push_back(hLayer);
1127 :
1128 2 : const GDALDataType eOutputType = GetOutputDataType(hLayer);
1129 2 : hDstDS = CreateOutputDataset(
1130 2 : ahLayers, hSRS, psOptions->sEnvelop, hDriver, pszDest,
1131 2 : psOptions->nXSize, psOptions->nYSize, psOptions->dfXRes,
1132 2 : psOptions->dfYRes, psOptions->bTargetAlignedPixels,
1133 2 : static_cast<int>(psOptions->anBandList.size()), eOutputType,
1134 2 : psOptions->aosCreationOptions, psOptions->adfInitVals,
1135 2 : psOptions->osNoData.c_str());
1136 2 : if (hDstDS == nullptr)
1137 : {
1138 0 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1139 0 : return nullptr;
1140 : }
1141 : }
1142 :
1143 18 : eErr = ProcessLayer(
1144 6 : hLayer, hSRS != nullptr, hDstDS, psOptions->anBandList,
1145 6 : psOptions->adfBurnValues, psOptions->b3D, psOptions->bInverse,
1146 6 : psOptions->osBurnAttribute.c_str(),
1147 6 : psOptions->aosRasterizeOptions, psOptions->aosTO,
1148 6 : psOptions->pfnProgress, psOptions->pProgressData);
1149 :
1150 6 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1151 : }
1152 : }
1153 :
1154 : /* -------------------------------------------------------------------- */
1155 : /* Create output file if necessary. */
1156 : /* -------------------------------------------------------------------- */
1157 : const int nLayerCount =
1158 104 : (psOptions->osSQL.empty() && psOptions->aosLayers.empty())
1159 110 : ? 1
1160 33 : : static_cast<int>(psOptions->aosLayers.size());
1161 :
1162 55 : if (bCreateOutput && hDstDS == nullptr)
1163 : {
1164 24 : std::vector<OGRLayerH> ahLayers;
1165 :
1166 24 : GDALDataType eOutputType = psOptions->eOutputType;
1167 :
1168 47 : for (int i = 0; i < nLayerCount; i++)
1169 : {
1170 : OGRLayerH hLayer;
1171 24 : if (psOptions->aosLayers.size() > static_cast<size_t>(i))
1172 13 : hLayer = GDALDatasetGetLayerByName(
1173 13 : hSrcDataset, psOptions->aosLayers[i].c_str());
1174 : else
1175 11 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1176 24 : if (hLayer == nullptr)
1177 : {
1178 1 : CPLError(CE_Failure, CPLE_AppDefined,
1179 : "Unable to find layer \"%s\".",
1180 1 : psOptions->aosLayers.size() > static_cast<size_t>(i)
1181 1 : ? psOptions->aosLayers[i].c_str()
1182 : : "0");
1183 1 : return nullptr;
1184 : }
1185 23 : if (eOutputType == GDT_Unknown)
1186 : {
1187 21 : if (GetOutputDataType(hLayer) == GDT_Int64)
1188 1 : eOutputType = GDT_Int64;
1189 : }
1190 :
1191 23 : ahLayers.push_back(hLayer);
1192 : }
1193 :
1194 23 : if (eOutputType == GDT_Unknown)
1195 : {
1196 20 : eOutputType = GDT_Float64;
1197 : }
1198 :
1199 23 : hDstDS = CreateOutputDataset(
1200 23 : ahLayers, hSRS, psOptions->sEnvelop, hDriver, pszDest,
1201 23 : psOptions->nXSize, psOptions->nYSize, psOptions->dfXRes,
1202 23 : psOptions->dfYRes, psOptions->bTargetAlignedPixels,
1203 23 : static_cast<int>(psOptions->anBandList.size()), eOutputType,
1204 23 : psOptions->aosCreationOptions, psOptions->adfInitVals,
1205 23 : psOptions->osNoData.c_str());
1206 23 : if (hDstDS == nullptr)
1207 : {
1208 5 : return nullptr;
1209 : }
1210 : }
1211 :
1212 : /* -------------------------------------------------------------------- */
1213 : /* Process each layer. */
1214 : /* -------------------------------------------------------------------- */
1215 :
1216 91 : for (int i = 0; i < nLayerCount; i++)
1217 : {
1218 : OGRLayerH hLayer;
1219 43 : if (psOptions->aosLayers.size() > static_cast<size_t>(i))
1220 25 : hLayer = GDALDatasetGetLayerByName(hSrcDataset,
1221 25 : psOptions->aosLayers[i].c_str());
1222 : else
1223 18 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1224 43 : if (hLayer == nullptr)
1225 : {
1226 0 : CPLError(CE_Failure, CPLE_AppDefined,
1227 : "Unable to find layer \"%s\".",
1228 0 : psOptions->aosLayers.size() > static_cast<size_t>(i)
1229 0 : ? psOptions->aosLayers[i].c_str()
1230 : : "0");
1231 0 : eErr = CE_Failure;
1232 0 : break;
1233 : }
1234 :
1235 43 : if (!psOptions->osWHERE.empty())
1236 : {
1237 1 : if (OGR_L_SetAttributeFilter(hLayer, psOptions->osWHERE.c_str()) !=
1238 : OGRERR_NONE)
1239 : {
1240 0 : eErr = CE_Failure;
1241 0 : break;
1242 : }
1243 : }
1244 :
1245 86 : void *pScaledProgress = GDALCreateScaledProgress(
1246 43 : 0.0, 1.0 * (i + 1) / nLayerCount, psOptions->pfnProgress,
1247 43 : psOptions->pProgressData);
1248 :
1249 129 : eErr = ProcessLayer(hLayer, !psOptions->oOutputSRS.IsEmpty(), hDstDS,
1250 43 : psOptions->anBandList, psOptions->adfBurnValues,
1251 43 : psOptions->b3D, psOptions->bInverse,
1252 43 : psOptions->osBurnAttribute.c_str(),
1253 43 : psOptions->aosRasterizeOptions, psOptions->aosTO,
1254 : GDALScaledProgress, pScaledProgress);
1255 :
1256 43 : GDALDestroyScaledProgress(pScaledProgress);
1257 43 : if (eErr != CE_None)
1258 1 : break;
1259 : }
1260 :
1261 49 : if (eErr != CE_None)
1262 : {
1263 1 : if (bCloseOutDSOnError)
1264 0 : GDALClose(hDstDS);
1265 1 : return nullptr;
1266 : }
1267 :
1268 48 : return hDstDS;
1269 : }
1270 :
1271 : /************************************************************************/
1272 : /* ArgIsNumericRasterize() */
1273 : /************************************************************************/
1274 :
1275 359 : static bool ArgIsNumericRasterize(const char *pszArg)
1276 :
1277 : {
1278 359 : char *pszEnd = nullptr;
1279 359 : CPLStrtod(pszArg, &pszEnd);
1280 359 : return pszEnd != nullptr && pszEnd[0] == '\0';
1281 : }
1282 :
1283 : /************************************************************************/
1284 : /* GDALRasterizeOptionsNew() */
1285 : /************************************************************************/
1286 :
1287 : /**
1288 : * Allocates a GDALRasterizeOptions struct.
1289 : *
1290 : * @param papszArgv NULL terminated list of options (potentially including
1291 : * filename and open options too), or NULL. The accepted options are the ones of
1292 : * the <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
1293 : * @param psOptionsForBinary (output) may be NULL (and should generally be
1294 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
1295 : * GDALRasterizeOptionsForBinaryNew() prior to this
1296 : * function. Will be filled with potentially present filename, open options,...
1297 : * @return pointer to the allocated GDALRasterizeOptions struct. Must be freed
1298 : * with GDALRasterizeOptionsFree().
1299 : *
1300 : * @since GDAL 2.1
1301 : */
1302 :
1303 : GDALRasterizeOptions *
1304 62 : GDALRasterizeOptionsNew(char **papszArgv,
1305 : GDALRasterizeOptionsForBinary *psOptionsForBinary)
1306 : {
1307 :
1308 124 : auto psOptions = std::make_unique<GDALRasterizeOptions>();
1309 :
1310 : /*-------------------------------------------------------------------- */
1311 : /* Parse arguments. */
1312 : /*-------------------------------------------------------------------- */
1313 :
1314 124 : CPLStringList aosArgv;
1315 :
1316 : /* -------------------------------------------------------------------- */
1317 : /* Pre-processing for custom syntax that ArgumentParser does not */
1318 : /* support. */
1319 : /* -------------------------------------------------------------------- */
1320 62 : const int argc = CSLCount(papszArgv);
1321 611 : for (int i = 0; i < argc && papszArgv != nullptr && papszArgv[i] != nullptr;
1322 : i++)
1323 : {
1324 : // argparser will be confused if the value of a string argument
1325 : // starts with a negative sign.
1326 549 : if (EQUAL(papszArgv[i], "-a_nodata") && papszArgv[i + 1])
1327 : {
1328 5 : ++i;
1329 5 : psOptions->osNoData = papszArgv[i];
1330 5 : psOptions->bCreateOutput = true;
1331 : }
1332 :
1333 : // argparser is confused by arguments that have at_least_one
1334 : // cardinality, if they immediately precede positional arguments.
1335 544 : else if (EQUAL(papszArgv[i], "-burn") && papszArgv[i + 1])
1336 : {
1337 100 : if (strchr(papszArgv[i + 1], ' '))
1338 : {
1339 : const CPLStringList aosTokens(
1340 0 : CSLTokenizeString(papszArgv[i + 1]));
1341 0 : for (const char *pszToken : aosTokens)
1342 : {
1343 0 : psOptions->adfBurnValues.push_back(CPLAtof(pszToken));
1344 : }
1345 0 : i += 1;
1346 : }
1347 : else
1348 : {
1349 200 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1350 : {
1351 100 : psOptions->adfBurnValues.push_back(
1352 100 : CPLAtof(papszArgv[i + 1]));
1353 100 : i += 1;
1354 : }
1355 : }
1356 :
1357 : // Dummy value to make argparse happy, as at least one of
1358 : // -burn, -a or -3d is required
1359 100 : aosArgv.AddString("-burn");
1360 100 : aosArgv.AddString("0");
1361 : }
1362 444 : else if (EQUAL(papszArgv[i], "-init") && papszArgv[i + 1])
1363 : {
1364 12 : if (strchr(papszArgv[i + 1], ' '))
1365 : {
1366 : const CPLStringList aosTokens(
1367 0 : CSLTokenizeString(papszArgv[i + 1]));
1368 0 : for (const char *pszToken : aosTokens)
1369 : {
1370 0 : psOptions->adfInitVals.push_back(CPLAtof(pszToken));
1371 : }
1372 0 : i += 1;
1373 : }
1374 : else
1375 : {
1376 24 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1377 : {
1378 12 : psOptions->adfInitVals.push_back(CPLAtof(papszArgv[i + 1]));
1379 12 : i += 1;
1380 : }
1381 : }
1382 12 : psOptions->bCreateOutput = true;
1383 : }
1384 432 : else if (EQUAL(papszArgv[i], "-b") && papszArgv[i + 1])
1385 : {
1386 72 : if (strchr(papszArgv[i + 1], ' '))
1387 : {
1388 : const CPLStringList aosTokens(
1389 0 : CSLTokenizeString(papszArgv[i + 1]));
1390 0 : for (const char *pszToken : aosTokens)
1391 : {
1392 0 : psOptions->anBandList.push_back(atoi(pszToken));
1393 : }
1394 0 : i += 1;
1395 : }
1396 : else
1397 : {
1398 144 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1399 : {
1400 72 : psOptions->anBandList.push_back(atoi(papszArgv[i + 1]));
1401 72 : i += 1;
1402 : }
1403 72 : }
1404 : }
1405 : else
1406 : {
1407 360 : aosArgv.AddString(papszArgv[i]);
1408 : }
1409 : }
1410 :
1411 : try
1412 : {
1413 : auto argParser =
1414 64 : GDALRasterizeOptionsGetParser(psOptions.get(), psOptionsForBinary);
1415 62 : argParser->parse_args_without_binary_name(aosArgv.List());
1416 :
1417 : // Check all no store_into args
1418 63 : if (auto oTe = argParser->present<std::vector<double>>("-te"))
1419 : {
1420 3 : psOptions->sEnvelop.MinX = oTe.value()[0];
1421 3 : psOptions->sEnvelop.MinY = oTe.value()[1];
1422 3 : psOptions->sEnvelop.MaxX = oTe.value()[2];
1423 3 : psOptions->sEnvelop.MaxY = oTe.value()[3];
1424 3 : psOptions->bCreateOutput = true;
1425 : }
1426 :
1427 60 : if (auto oTr = argParser->present<std::vector<double>>("-tr"))
1428 : {
1429 13 : psOptions->dfXRes = oTr.value()[0];
1430 13 : psOptions->dfYRes = oTr.value()[1];
1431 :
1432 13 : if (psOptions->dfXRes <= 0 || psOptions->dfYRes <= 0)
1433 : {
1434 0 : CPLError(CE_Failure, CPLE_AppDefined,
1435 : "Wrong value for -tr parameter.");
1436 0 : return nullptr;
1437 : }
1438 :
1439 13 : psOptions->bCreateOutput = true;
1440 : }
1441 :
1442 60 : if (auto oTs = argParser->present<std::vector<double>>("-ts"))
1443 : {
1444 15 : const int nXSize = static_cast<int>(oTs.value()[0]);
1445 15 : const int nYSize = static_cast<int>(oTs.value()[1]);
1446 :
1447 : // Warn the user if the conversion to int looses precision
1448 15 : if (nXSize != oTs.value()[0] || nYSize != oTs.value()[1])
1449 : {
1450 1 : CPLError(CE_Warning, CPLE_AppDefined,
1451 : "-ts values parsed as %d %d.", nXSize, nYSize);
1452 : }
1453 :
1454 15 : psOptions->nXSize = nXSize;
1455 15 : psOptions->nYSize = nYSize;
1456 :
1457 15 : if (psOptions->nXSize <= 0 || psOptions->nYSize <= 0)
1458 : {
1459 0 : CPLError(CE_Failure, CPLE_AppDefined,
1460 : "Wrong value for -ts parameter.");
1461 0 : return nullptr;
1462 : }
1463 :
1464 15 : psOptions->bCreateOutput = true;
1465 : }
1466 :
1467 60 : if (psOptions->bCreateOutput)
1468 : {
1469 45 : if (psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
1470 45 : psOptions->nXSize == 0 && psOptions->nYSize == 0)
1471 : {
1472 1 : CPLError(CE_Failure, CPLE_NotSupported,
1473 : "'-tr xres yres' or '-ts xsize ysize' is required.");
1474 1 : return nullptr;
1475 : }
1476 :
1477 28 : if (psOptions->bTargetAlignedPixels && psOptions->dfXRes == 0 &&
1478 0 : psOptions->dfYRes == 0)
1479 : {
1480 0 : CPLError(CE_Failure, CPLE_NotSupported,
1481 : "-tap option cannot be used without using -tr.");
1482 0 : return nullptr;
1483 : }
1484 :
1485 28 : if (!psOptions->anBandList.empty())
1486 : {
1487 1 : CPLError(
1488 : CE_Failure, CPLE_NotSupported,
1489 : "-b option cannot be used when creating a GDAL dataset.");
1490 1 : return nullptr;
1491 : }
1492 :
1493 27 : int nBandCount = 1;
1494 :
1495 27 : if (!psOptions->adfBurnValues.empty())
1496 13 : nBandCount = static_cast<int>(psOptions->adfBurnValues.size());
1497 :
1498 27 : if (static_cast<int>(psOptions->adfInitVals.size()) > nBandCount)
1499 0 : nBandCount = static_cast<int>(psOptions->adfInitVals.size());
1500 :
1501 27 : if (psOptions->adfInitVals.size() == 1)
1502 : {
1503 3 : for (int i = 1; i <= nBandCount - 1; i++)
1504 0 : psOptions->adfInitVals.push_back(psOptions->adfInitVals[0]);
1505 : }
1506 :
1507 70 : for (int i = 1; i <= nBandCount; i++)
1508 43 : psOptions->anBandList.push_back(i);
1509 : }
1510 : else
1511 : {
1512 31 : if (psOptions->anBandList.empty())
1513 11 : psOptions->anBandList.push_back(1);
1514 : }
1515 :
1516 58 : if (!psOptions->osDialect.empty() && !psOptions->osWHERE.empty() &&
1517 0 : !psOptions->osSQL.empty())
1518 : {
1519 0 : CPLError(CE_Warning, CPLE_AppDefined,
1520 : "-dialect is ignored with -where. Use -sql instead");
1521 : }
1522 :
1523 58 : if (psOptionsForBinary)
1524 : {
1525 11 : psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
1526 11 : if (!psOptions->osFormat.empty())
1527 0 : psOptionsForBinary->osFormat = psOptions->osFormat;
1528 : }
1529 63 : else if (psOptions->adfBurnValues.empty() &&
1530 63 : psOptions->osBurnAttribute.empty() && !psOptions->b3D)
1531 : {
1532 11 : psOptions->adfBurnValues.push_back(255);
1533 : }
1534 : }
1535 2 : catch (const std::exception &e)
1536 : {
1537 2 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
1538 2 : return nullptr;
1539 : }
1540 :
1541 58 : return psOptions.release();
1542 : }
1543 :
1544 : /************************************************************************/
1545 : /* GDALRasterizeOptionsFree() */
1546 : /************************************************************************/
1547 :
1548 : /**
1549 : * Frees the GDALRasterizeOptions struct.
1550 : *
1551 : * @param psOptions the options struct for GDALRasterize().
1552 : *
1553 : * @since GDAL 2.1
1554 : */
1555 :
1556 58 : void GDALRasterizeOptionsFree(GDALRasterizeOptions *psOptions)
1557 : {
1558 58 : delete psOptions;
1559 58 : }
1560 :
1561 : /************************************************************************/
1562 : /* GDALRasterizeOptionsSetProgress() */
1563 : /************************************************************************/
1564 :
1565 : /**
1566 : * Set a progress function.
1567 : *
1568 : * @param psOptions the options struct for GDALRasterize().
1569 : * @param pfnProgress the progress callback.
1570 : * @param pProgressData the user data for the progress callback.
1571 : *
1572 : * @since GDAL 2.1
1573 : */
1574 :
1575 34 : void GDALRasterizeOptionsSetProgress(GDALRasterizeOptions *psOptions,
1576 : GDALProgressFunc pfnProgress,
1577 : void *pProgressData)
1578 : {
1579 34 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
1580 34 : psOptions->pProgressData = pProgressData;
1581 34 : }
|
