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