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 66 : GDALRasterizeOptionsGetParser(GDALRasterizeOptions *psOptions,
80 : GDALRasterizeOptionsForBinary *psOptionsForBinary)
81 : {
82 : auto argParser = std::make_unique<GDALArgumentParser>(
83 66 : "gdal_rasterize", /* bForBinary=*/psOptionsForBinary != nullptr);
84 :
85 66 : argParser->add_description(_("Burns vector geometries into a raster."));
86 :
87 66 : argParser->add_epilog(
88 : _("This program burns vector geometries (points, lines, and polygons) "
89 66 : "into the raster band(s) of a raster image."));
90 :
91 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
92 66 : argParser->add_argument("-b")
93 132 : .metavar("<band>")
94 66 : .append()
95 66 : .scan<'i', int>()
96 : //.nargs(argparse::nargs_pattern::at_least_one)
97 66 : .help(_("The band(s) to burn values into."));
98 :
99 66 : argParser->add_argument("-i")
100 66 : .flag()
101 66 : .store_into(psOptions->bInverse)
102 66 : .help(_("Invert rasterization."));
103 :
104 66 : argParser->add_argument("-at")
105 66 : .flag()
106 : .action(
107 32 : [psOptions](const std::string &) {
108 : psOptions->aosRasterizeOptions.SetNameValue("ALL_TOUCHED",
109 32 : "TRUE");
110 66 : })
111 66 : .help(_("Enables the ALL_TOUCHED rasterization option."));
112 :
113 : // Mutually exclusive options: -burn, -3d, -a
114 : {
115 : // Required if options for binary
116 66 : auto &group = argParser->add_mutually_exclusive_group(
117 66 : psOptionsForBinary != nullptr);
118 :
119 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
120 66 : group.add_argument("-burn")
121 132 : .metavar("<value>")
122 66 : .scan<'g', double>()
123 66 : .append()
124 : //.nargs(argparse::nargs_pattern::at_least_one)
125 66 : .help(_("A fixed value to burn into the raster band(s)."));
126 :
127 66 : group.add_argument("-a")
128 132 : .metavar("<attribute_name>")
129 66 : .store_into(psOptions->osBurnAttribute)
130 : .help(_("Name of the field in the input layer to get the burn "
131 66 : "values from."));
132 :
133 66 : group.add_argument("-3d")
134 66 : .flag()
135 66 : .store_into(psOptions->b3D)
136 : .action(
137 5 : [psOptions](const std::string &) {
138 : psOptions->aosRasterizeOptions.SetNameValue(
139 5 : "BURN_VALUE_FROM", "Z");
140 66 : })
141 : .help(_("Indicates that a burn value should be extracted from the "
142 66 : "\"Z\" values of the feature."));
143 : }
144 :
145 66 : argParser->add_argument("-add")
146 66 : .flag()
147 : .action(
148 1 : [psOptions](const std::string &) {
149 1 : psOptions->aosRasterizeOptions.SetNameValue("MERGE_ALG", "ADD");
150 66 : })
151 : .help(_("Instead of burning a new value, this adds the new value to "
152 66 : "the existing raster."));
153 :
154 : // Undocumented
155 66 : argParser->add_argument("-chunkysize")
156 66 : .flag()
157 66 : .hidden()
158 : .action(
159 0 : [psOptions](const std::string &s) {
160 : psOptions->aosRasterizeOptions.SetNameValue("CHUNKYSIZE",
161 0 : s.c_str());
162 66 : });
163 :
164 : // Mutually exclusive -l, -sql
165 : {
166 66 : auto &group = argParser->add_mutually_exclusive_group(false);
167 :
168 66 : group.add_argument("-l")
169 132 : .metavar("<layer_name>")
170 66 : .append()
171 66 : .store_into(psOptions->aosLayers)
172 66 : .help(_("Name of the layer(s) to process."));
173 :
174 66 : group.add_argument("-sql")
175 132 : .metavar("<sql_statement>")
176 66 : .store_into(psOptions->osSQL)
177 : .action(
178 10 : [psOptions](const std::string &sql)
179 : {
180 9 : GByte *pabyRet = nullptr;
181 10 : if (!sql.empty() && sql.at(0) == '@' &&
182 10 : VSIIngestFile(nullptr, sql.substr(1).c_str(), &pabyRet,
183 : nullptr, 10 * 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 75 : })
193 : .help(
194 : _("An SQL statement to be evaluated against the datasource to "
195 66 : "produce a virtual layer of features to be burned in."));
196 : }
197 :
198 66 : argParser->add_argument("-where")
199 132 : .metavar("<expression>")
200 66 : .store_into(psOptions->osWHERE)
201 : .help(_("An optional SQL WHERE style query expression to be applied to "
202 : "select features "
203 66 : "to burn in from the input layer(s)."));
204 :
205 66 : argParser->add_argument("-dialect")
206 132 : .metavar("<sql_dialect>")
207 66 : .store_into(psOptions->osDialect)
208 66 : .help(_("The SQL dialect to use for the SQL expression."));
209 :
210 : // Store later
211 66 : argParser->add_argument("-a_nodata")
212 132 : .metavar("<value>")
213 66 : .help(_("Assign a specified nodata value to output bands."));
214 :
215 : // Dealt manually as argparse::nargs_pattern::at_least_one is problematic
216 66 : argParser->add_argument("-init")
217 132 : .metavar("<value>")
218 66 : .append()
219 : //.nargs(argparse::nargs_pattern::at_least_one)
220 66 : .scan<'g', double>()
221 66 : .help(_("Initialize the output bands to the specified value."));
222 :
223 66 : argParser->add_argument("-a_srs")
224 132 : .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 68 : })
237 66 : .help(_("The spatial reference system to use for the output raster."));
238 :
239 66 : argParser->add_argument("-to")
240 132 : .metavar("<NAME>=<VALUE>")
241 66 : .append()
242 1 : .action([psOptions](const std::string &s)
243 67 : { psOptions->aosTO.AddString(s.c_str()); })
244 66 : .help(_("Set a transformer option."));
245 :
246 : // Store later
247 66 : argParser->add_argument("-te")
248 132 : .metavar("<xmin> <ymin> <xmax> <ymax>")
249 66 : .nargs(4)
250 66 : .scan<'g', double>()
251 66 : .help(_("Set georeferenced extents of output file to be created."));
252 :
253 : // Mutex with tr
254 : {
255 66 : auto &group = argParser->add_mutually_exclusive_group(false);
256 :
257 : // Store later
258 66 : group.add_argument("-tr")
259 132 : .metavar("<xres> <yres>")
260 66 : .nargs(2)
261 66 : .scan<'g', double>()
262 : .help(
263 66 : _("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 66 : auto &arg = group.add_argument("-ts")
269 132 : .metavar("<width> <height>")
270 66 : .nargs(2)
271 66 : .scan<'g', double>()
272 66 : .help(_("Set output file size in pixels and lines."));
273 :
274 66 : argParser->add_hidden_alias_for(arg, "-outsize");
275 : }
276 :
277 66 : argParser->add_argument("-tap")
278 66 : .flag()
279 66 : .store_into(psOptions->bTargetAlignedPixels)
280 1 : .action([psOptions](const std::string &)
281 66 : { psOptions->bCreateOutput = true; })
282 : .help(_("Align the coordinates of the extent to the values of the "
283 66 : "output raster."));
284 :
285 66 : argParser->add_argument("-optim")
286 132 : .metavar("AUTO|VECTOR|RASTER")
287 : .action(
288 35 : [psOptions](const std::string &s) {
289 35 : psOptions->aosRasterizeOptions.SetNameValue("OPTIM", s.c_str());
290 66 : })
291 66 : .help(_("Force the algorithm used."));
292 :
293 66 : argParser->add_creation_options_argument(psOptions->aosCreationOptions)
294 2 : .action([psOptions](const std::string &)
295 66 : { psOptions->bCreateOutput = true; });
296 :
297 66 : argParser->add_output_type_argument(psOptions->eOutputType)
298 2 : .action([psOptions](const std::string &)
299 66 : { psOptions->bCreateOutput = true; });
300 :
301 66 : argParser->add_output_format_argument(psOptions->osFormat)
302 11 : .action([psOptions](const std::string &)
303 66 : { 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 66 : psOptionsForBinary ? &(psOptionsForBinary->bQuiet) : nullptr);
309 :
310 66 : 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 66 : 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 54 : 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 54 : OGRCoordinateTransformationH hCT = nullptr;
527 54 : if (!bSRSIsSet)
528 : {
529 52 : OGRSpatialReferenceH hDstSRS = GDALGetSpatialRef(hDstDS);
530 :
531 52 : if (hDstSRS)
532 13 : hDstSRS = OSRClone(hDstSRS);
533 39 : 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 52 : OGRSpatialReferenceH hSrcSRS = OGR_L_GetSpatialRef(hSrcLayer);
542 52 : 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 37 : 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 37 : 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 52 : if (hDstSRS != nullptr)
575 : {
576 15 : OSRDestroySpatialReference(hDstSRS);
577 : }
578 : }
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* Get field index, and check. */
582 : /* -------------------------------------------------------------------- */
583 54 : int iBurnField = -1;
584 54 : bool bUseInt64 = false;
585 54 : 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 53 : OGRFeatureH hFeat = nullptr;
615 106 : std::vector<OGRGeometryH> ahGeometries;
616 106 : std::vector<double> adfFullBurnValues;
617 53 : std::vector<int64_t> anFullBurnValues;
618 :
619 53 : OGR_L_ResetReading(hSrcLayer);
620 :
621 2046 : while ((hFeat = OGR_L_GetNextFeature(hSrcLayer)) != nullptr)
622 : {
623 1993 : OGRGeometryH hGeom = OGR_F_StealGeometry(hFeat);
624 1993 : if (hGeom == nullptr)
625 : {
626 5 : OGR_F_Destroy(hFeat);
627 5 : continue;
628 : }
629 :
630 1988 : 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 1988 : ahGeometries.push_back(hGeom);
640 :
641 4132 : for (unsigned int iBand = 0; iBand < anBandList.size(); iBand++)
642 : {
643 2144 : if (!adfBurnValues.empty())
644 265 : adfFullBurnValues.push_back(adfBurnValues[std::min(
645 : iBand,
646 530 : 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 1988 : OGR_F_Destroy(hFeat);
666 : }
667 :
668 53 : 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 53 : 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 53 : void *pTransformArg = nullptr;
708 53 : GDALTransformerFunc pfnTransformer = nullptr;
709 53 : CPLErr eErr = CE_None;
710 53 : if (papszTO != nullptr)
711 : {
712 1 : GDALDataset *poDS = GDALDataset::FromHandle(hDstDS);
713 1 : char **papszTransformerOptions = CSLDuplicate(papszTO);
714 1 : GDALGeoTransform gt;
715 2 : if (poDS->GetGeoTransform(gt) != CE_None && poDS->GetGCPCount() == 0 &&
716 1 : 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 53 : if (eErr == CE_None)
737 : {
738 53 : 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 104 : eErr = GDALRasterizeGeometries(
749 52 : hDstDS, static_cast<int>(anBandList.size()), anBandList.data(),
750 52 : static_cast<int>(ahGeometries.size()), ahGeometries.data(),
751 52 : pfnTransformer, pTransformArg, adfFullBurnValues.data(),
752 : papszRasterizeOptions, pfnProgress, pProgressData);
753 : }
754 : }
755 :
756 : /* -------------------------------------------------------------------- */
757 : /* Cleanup */
758 : /* -------------------------------------------------------------------- */
759 :
760 53 : if (pTransformArg)
761 1 : GDALDestroyTransformer(pTransformArg);
762 :
763 2037 : for (int iGeom = static_cast<int>(ahGeometries.size()) - 1; iGeom >= 0;
764 : iGeom--)
765 1984 : OGR_G_DestroyGeometry(ahGeometries[iGeom]);
766 :
767 53 : return eErr;
768 : }
769 :
770 : /************************************************************************/
771 : /* CreateOutputDataset() */
772 : /************************************************************************/
773 :
774 29 : 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 29 : bool bFirstLayer = true;
782 29 : char *pszWKT = nullptr;
783 29 : const bool bBoundsSpecifiedByUser = sEnvelop.IsInit();
784 :
785 57 : for (unsigned int i = 0; i < ahLayers.size(); i++)
786 : {
787 29 : OGRLayerH hLayer = ahLayers[i];
788 :
789 29 : if (!bBoundsSpecifiedByUser)
790 : {
791 26 : OGREnvelope sLayerEnvelop;
792 :
793 26 : 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 25 : 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 25 : sEnvelop.Merge(sLayerEnvelop);
811 : }
812 :
813 28 : if (bFirstLayer)
814 : {
815 28 : if (hSRS == nullptr)
816 26 : hSRS = OGR_L_GetSpatialRef(hLayer);
817 :
818 28 : bFirstLayer = false;
819 : }
820 : }
821 :
822 28 : if (!sEnvelop.IsInit())
823 : {
824 0 : CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
825 0 : return nullptr;
826 : }
827 :
828 28 : if (dfXRes == 0 && dfYRes == 0)
829 : {
830 15 : if (nXSize == 0 || nYSize == 0)
831 : {
832 1 : CPLError(CE_Failure, CPLE_AppDefined,
833 : "Size and resolution are missing");
834 1 : return nullptr;
835 : }
836 14 : dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;
837 14 : 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 27 : if (dfXRes == 0 || dfYRes == 0)
848 : {
849 0 : CPLError(CE_Failure, CPLE_AppDefined, "Could not determine bounds");
850 0 : return nullptr;
851 : }
852 :
853 27 : double adfProjection[6] = {sEnvelop.MinX, dfXRes, 0.0,
854 27 : sEnvelop.MaxY, 0.0, -dfYRes};
855 :
856 27 : 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 25 : GDALCreate(hDriver, pszDest, nXSize, nYSize, nBandCount, eOutputType,
878 : papszCreationOptions);
879 25 : if (hDstDS == nullptr)
880 : {
881 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszDest);
882 0 : return nullptr;
883 : }
884 :
885 25 : GDALSetGeoTransform(hDstDS, adfProjection);
886 :
887 25 : if (hSRS)
888 6 : OSRExportToWkt(hSRS, &pszWKT);
889 25 : if (pszWKT)
890 6 : GDALSetProjection(hDstDS, pszWKT);
891 25 : 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 25 : if (pszNoData)
904 : {
905 76 : for (int iBand = 0; iBand < nBandCount; iBand++)
906 : {
907 51 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
908 51 : if (GDALGetRasterDataType(hBand) == GDT_Int64)
909 1 : GDALSetRasterNoDataValueAsInt64(hBand,
910 1 : CPLAtoGIntBig(pszNoData));
911 : else
912 50 : GDALSetRasterNoDataValue(hBand, CPLAtof(pszNoData));
913 : }
914 : }
915 :
916 25 : if (!adfInitVals.empty())
917 : {
918 30 : for (int iBand = 0;
919 30 : iBand < std::min(nBandCount, static_cast<int>(adfInitVals.size()));
920 : iBand++)
921 : {
922 21 : GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);
923 21 : GDALFillRaster(hBand, adfInitVals[iBand], 0);
924 : }
925 : }
926 :
927 25 : 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 63 : GDALDatasetH GDALRasterize(const char *pszDest, GDALDatasetH hDstDS,
960 : GDALDatasetH hSrcDataset,
961 : const GDALRasterizeOptions *psOptionsIn,
962 : int *pbUsageError)
963 : {
964 63 : 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 63 : 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 63 : 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 126 : psOptionsToFree(nullptr, GDALRasterizeOptionsFree);
994 126 : GDALRasterizeOptions sOptions;
995 63 : if (psOptionsIn)
996 63 : sOptions = *psOptionsIn;
997 :
998 63 : const bool bCloseOutDSOnError = hDstDS == nullptr;
999 63 : if (pszDest == nullptr)
1000 13 : pszDest = GDALGetDescription(hDstDS);
1001 :
1002 88 : if (sOptions.osSQL.empty() && sOptions.aosLayers.empty() &&
1003 25 : GDALDatasetGetLayerCount(hSrcDataset) != 1)
1004 : {
1005 0 : CPLError(CE_Failure, CPLE_NotSupported,
1006 : "Neither -sql nor -l are specified, but the source dataset "
1007 : "has not one single layer.");
1008 0 : if (pbUsageError)
1009 0 : *pbUsageError = TRUE;
1010 0 : return nullptr;
1011 : }
1012 :
1013 : /* -------------------------------------------------------------------- */
1014 : /* Open target raster file. Eventually we will add optional */
1015 : /* creation. */
1016 : /* -------------------------------------------------------------------- */
1017 63 : const bool bCreateOutput = sOptions.bCreateOutput || hDstDS == nullptr;
1018 :
1019 63 : GDALDriverH hDriver = nullptr;
1020 63 : if (bCreateOutput)
1021 : {
1022 34 : CPLString osFormat;
1023 34 : if (sOptions.osFormat.empty())
1024 : {
1025 23 : osFormat = GetOutputDriverForRaster(pszDest);
1026 23 : if (osFormat.empty())
1027 : {
1028 0 : return nullptr;
1029 : }
1030 : }
1031 : else
1032 : {
1033 11 : osFormat = sOptions.osFormat;
1034 : }
1035 :
1036 : /* --------------------------------------------------------------------
1037 : */
1038 : /* Find the output driver. */
1039 : /* --------------------------------------------------------------------
1040 : */
1041 34 : hDriver = GDALGetDriverByName(osFormat);
1042 : char **papszDriverMD =
1043 34 : hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;
1044 34 : if (hDriver == nullptr)
1045 : {
1046 1 : CPLError(CE_Failure, CPLE_NotSupported,
1047 : "Output driver `%s' not recognised.", osFormat.c_str());
1048 1 : return nullptr;
1049 : }
1050 33 : if (!CPLTestBool(
1051 : CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_RASTER, "FALSE")))
1052 : {
1053 1 : CPLError(CE_Failure, CPLE_NotSupported,
1054 : "Output driver `%s' is not a raster driver.",
1055 : osFormat.c_str());
1056 1 : return nullptr;
1057 : }
1058 32 : if (!CPLTestBool(
1059 : CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))
1060 : {
1061 1 : CPLError(CE_Failure, CPLE_NotSupported,
1062 : "Output driver `%s' does not support direct output file "
1063 : "creation. "
1064 : "To write a file to this format, first write to a "
1065 : "different format such as "
1066 : "GeoTIFF and then convert the output.",
1067 : osFormat.c_str());
1068 1 : return nullptr;
1069 : }
1070 : }
1071 :
1072 4 : auto calculateSize = [&](const OGREnvelope &sEnvelope) -> bool
1073 : {
1074 4 : const double width{sEnvelope.MaxX - sEnvelope.MinX};
1075 4 : if (std::isnan(width))
1076 : {
1077 0 : return false;
1078 : }
1079 :
1080 4 : const double height{sEnvelope.MaxY - sEnvelope.MinY};
1081 4 : if (std::isnan(height))
1082 : {
1083 0 : return false;
1084 : }
1085 :
1086 4 : if (height == 0 || width == 0)
1087 : {
1088 0 : return false;
1089 : }
1090 :
1091 4 : if (sOptions.nXSize == 0)
1092 : {
1093 2 : const double xSize{
1094 2 : (sEnvelope.MaxX - sEnvelope.MinX) /
1095 2 : ((sEnvelope.MaxY - sEnvelope.MinY) / sOptions.nYSize)};
1096 4 : if (std::isnan(xSize) || xSize > std::numeric_limits<int>::max() ||
1097 2 : xSize < std::numeric_limits<int>::min())
1098 : {
1099 0 : return false;
1100 : }
1101 2 : sOptions.nXSize = static_cast<int>(xSize);
1102 : }
1103 : else
1104 : {
1105 2 : const double ySize{
1106 2 : (sEnvelope.MaxY - sEnvelope.MinY) /
1107 2 : ((sEnvelope.MaxX - sEnvelope.MinX) / sOptions.nXSize)};
1108 4 : if (std::isnan(ySize) || ySize > std::numeric_limits<int>::max() ||
1109 2 : ySize < std::numeric_limits<int>::min())
1110 : {
1111 0 : return false;
1112 : }
1113 2 : sOptions.nYSize = static_cast<int>(ySize);
1114 : }
1115 4 : return sOptions.nXSize > 0 && sOptions.nYSize > 0;
1116 60 : };
1117 :
1118 : const int nLayerCount =
1119 111 : (sOptions.osSQL.empty() && sOptions.aosLayers.empty())
1120 120 : ? 1
1121 38 : : static_cast<int>(sOptions.aosLayers.size());
1122 :
1123 60 : const bool bOneSizeNeedsCalculation{
1124 60 : static_cast<bool>((sOptions.nXSize == 0) ^ (sOptions.nYSize == 0))};
1125 :
1126 : // Calculate the size if either nXSize or nYSize is 0
1127 60 : if (sOptions.osSQL.empty() && bOneSizeNeedsCalculation)
1128 : {
1129 2 : CPLErr eErr = CE_None;
1130 : // Get the extent of the source dataset
1131 2 : OGREnvelope sEnvelope;
1132 2 : bool bFirstLayer = true;
1133 4 : for (int i = 0; i < nLayerCount; i++)
1134 : {
1135 : OGRLayerH hLayer;
1136 2 : if (sOptions.aosLayers.size() > static_cast<size_t>(i))
1137 2 : hLayer = GDALDatasetGetLayerByName(
1138 2 : hSrcDataset, sOptions.aosLayers[i].c_str());
1139 : else
1140 0 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1141 2 : if (hLayer == nullptr)
1142 : {
1143 0 : CPLError(CE_Failure, CPLE_AppDefined,
1144 : "Unable to find layer \"%s\".",
1145 0 : sOptions.aosLayers.size() > static_cast<size_t>(i)
1146 0 : ? sOptions.aosLayers[i].c_str()
1147 : : "0");
1148 0 : eErr = CE_Failure;
1149 0 : break;
1150 : }
1151 2 : OGREnvelope sLayerEnvelop;
1152 2 : if (OGR_L_GetExtent(hLayer, &sLayerEnvelop, TRUE) != OGRERR_NONE)
1153 : {
1154 0 : CPLError(CE_Failure, CPLE_AppDefined,
1155 : "Cannot get layer extent");
1156 0 : eErr = CE_Failure;
1157 0 : break;
1158 : }
1159 2 : if (bFirstLayer)
1160 : {
1161 2 : sEnvelope = sLayerEnvelop;
1162 2 : bFirstLayer = false;
1163 : }
1164 : else
1165 : {
1166 0 : sEnvelope.Merge(sLayerEnvelop);
1167 : }
1168 : }
1169 :
1170 2 : if (!calculateSize(sEnvelope))
1171 : {
1172 0 : CPLError(CE_Failure, CPLE_AppDefined,
1173 : "Cannot calculate size from layer extent");
1174 0 : eErr = CE_Failure;
1175 : }
1176 :
1177 2 : if (eErr == CE_Failure)
1178 : {
1179 0 : if (bCloseOutDSOnError && hDstDS)
1180 : {
1181 0 : GDALClose(hDstDS);
1182 : }
1183 0 : return nullptr;
1184 : }
1185 : }
1186 :
1187 27 : const auto GetOutputDataType = [&](OGRLayerH hLayer)
1188 : {
1189 27 : CPLAssert(bCreateOutput);
1190 27 : CPLAssert(hDriver);
1191 55 : GDALDataType eOutputType = sOptions.eOutputType;
1192 27 : if (eOutputType == GDT_Unknown && !sOptions.osBurnAttribute.empty())
1193 : {
1194 1 : OGRFeatureDefnH hLayerDefn = OGR_L_GetLayerDefn(hLayer);
1195 1 : const int iBurnField = OGR_FD_GetFieldIndex(
1196 : hLayerDefn, sOptions.osBurnAttribute.c_str());
1197 1 : if (iBurnField >= 0 && OGR_Fld_GetType(OGR_FD_GetFieldDefn(
1198 : hLayerDefn, iBurnField)) == OFTInteger64)
1199 : {
1200 1 : const char *pszMD = GDALGetMetadataItem(
1201 : hDriver, GDAL_DMD_CREATIONDATATYPES, nullptr);
1202 2 : if (pszMD && CPLStringList(CSLTokenizeString2(pszMD, " ", 0))
1203 1 : .FindString("Int64") >= 0)
1204 : {
1205 1 : eOutputType = GDT_Int64;
1206 : }
1207 : }
1208 : }
1209 27 : if (eOutputType == GDT_Unknown)
1210 : {
1211 26 : eOutputType = GDT_Float64;
1212 : }
1213 27 : return eOutputType;
1214 60 : };
1215 :
1216 : // Store SRS handle
1217 : OGRSpatialReferenceH hSRS =
1218 60 : sOptions.oOutputSRS.IsEmpty()
1219 60 : ? nullptr
1220 2 : : OGRSpatialReference::ToHandle(
1221 60 : const_cast<OGRSpatialReference *>(&sOptions.oOutputSRS));
1222 :
1223 : /* -------------------------------------------------------------------- */
1224 : /* Process SQL request. */
1225 : /* -------------------------------------------------------------------- */
1226 60 : CPLErr eErr = CE_Failure;
1227 :
1228 60 : if (!sOptions.osSQL.empty())
1229 : {
1230 : OGRLayerH hLayer =
1231 9 : GDALDatasetExecuteSQL(hSrcDataset, sOptions.osSQL.c_str(), nullptr,
1232 9 : sOptions.osDialect.c_str());
1233 9 : if (hLayer != nullptr)
1234 : {
1235 :
1236 9 : if (bOneSizeNeedsCalculation)
1237 : {
1238 3 : OGREnvelope sEnvelope;
1239 : bool bSizeCalculationError{
1240 3 : OGR_L_GetExtent(hLayer, &sEnvelope, TRUE) != OGRERR_NONE};
1241 3 : if (!bSizeCalculationError)
1242 : {
1243 2 : bSizeCalculationError = !calculateSize(sEnvelope);
1244 : }
1245 :
1246 3 : if (bSizeCalculationError)
1247 : {
1248 1 : CPLError(CE_Failure, CPLE_AppDefined,
1249 : "Cannot get layer extent");
1250 1 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1251 1 : return nullptr;
1252 : }
1253 : }
1254 :
1255 8 : if (bCreateOutput)
1256 : {
1257 4 : std::vector<OGRLayerH> ahLayers;
1258 4 : ahLayers.push_back(hLayer);
1259 :
1260 4 : const GDALDataType eOutputType = GetOutputDataType(hLayer);
1261 8 : hDstDS = CreateOutputDataset(
1262 : ahLayers, hSRS, sOptions.sEnvelop, hDriver, pszDest,
1263 : sOptions.nXSize, sOptions.nYSize, sOptions.dfXRes,
1264 4 : sOptions.dfYRes, sOptions.bTargetAlignedPixels,
1265 4 : static_cast<int>(sOptions.anBandList.size()), eOutputType,
1266 : sOptions.aosCreationOptions, sOptions.adfInitVals,
1267 4 : sOptions.osNoData.c_str());
1268 4 : if (hDstDS == nullptr)
1269 : {
1270 0 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1271 0 : return nullptr;
1272 : }
1273 : }
1274 :
1275 16 : eErr = ProcessLayer(
1276 : hLayer, hSRS != nullptr, hDstDS, sOptions.anBandList,
1277 8 : sOptions.adfBurnValues, sOptions.b3D, sOptions.bInverse,
1278 : sOptions.osBurnAttribute.c_str(), sOptions.aosRasterizeOptions,
1279 8 : sOptions.aosTO, sOptions.pfnProgress, sOptions.pProgressData);
1280 :
1281 8 : GDALDatasetReleaseResultSet(hSrcDataset, hLayer);
1282 : }
1283 : }
1284 :
1285 : /* -------------------------------------------------------------------- */
1286 : /* Create output file if necessary. */
1287 : /* -------------------------------------------------------------------- */
1288 :
1289 59 : if (bCreateOutput && hDstDS == nullptr)
1290 : {
1291 26 : std::vector<OGRLayerH> ahLayers;
1292 :
1293 26 : GDALDataType eOutputType = sOptions.eOutputType;
1294 :
1295 51 : for (int i = 0; i < nLayerCount; i++)
1296 : {
1297 : OGRLayerH hLayer;
1298 26 : if (sOptions.aosLayers.size() > static_cast<size_t>(i))
1299 15 : hLayer = GDALDatasetGetLayerByName(
1300 15 : hSrcDataset, sOptions.aosLayers[i].c_str());
1301 : else
1302 11 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1303 26 : if (hLayer == nullptr)
1304 : {
1305 1 : CPLError(CE_Failure, CPLE_AppDefined,
1306 : "Unable to find layer \"%s\".",
1307 1 : sOptions.aosLayers.size() > static_cast<size_t>(i)
1308 1 : ? sOptions.aosLayers[i].c_str()
1309 : : "0");
1310 1 : return nullptr;
1311 : }
1312 25 : if (eOutputType == GDT_Unknown)
1313 : {
1314 23 : if (GetOutputDataType(hLayer) == GDT_Int64)
1315 1 : eOutputType = GDT_Int64;
1316 : }
1317 :
1318 25 : ahLayers.push_back(hLayer);
1319 : }
1320 :
1321 25 : if (eOutputType == GDT_Unknown)
1322 : {
1323 22 : eOutputType = GDT_Float64;
1324 : }
1325 :
1326 50 : hDstDS = CreateOutputDataset(
1327 : ahLayers, hSRS, sOptions.sEnvelop, hDriver, pszDest,
1328 : sOptions.nXSize, sOptions.nYSize, sOptions.dfXRes, sOptions.dfYRes,
1329 25 : sOptions.bTargetAlignedPixels,
1330 25 : static_cast<int>(sOptions.anBandList.size()), eOutputType,
1331 : sOptions.aosCreationOptions, sOptions.adfInitVals,
1332 25 : sOptions.osNoData.c_str());
1333 25 : if (hDstDS == nullptr)
1334 : {
1335 4 : return nullptr;
1336 : }
1337 : }
1338 :
1339 : /* -------------------------------------------------------------------- */
1340 : /* Process each layer. */
1341 : /* -------------------------------------------------------------------- */
1342 :
1343 99 : for (int i = 0; i < nLayerCount; i++)
1344 : {
1345 : OGRLayerH hLayer;
1346 46 : if (sOptions.aosLayers.size() > static_cast<size_t>(i))
1347 28 : hLayer = GDALDatasetGetLayerByName(hSrcDataset,
1348 28 : sOptions.aosLayers[i].c_str());
1349 : else
1350 18 : hLayer = GDALDatasetGetLayer(hSrcDataset, 0);
1351 46 : if (hLayer == nullptr)
1352 : {
1353 0 : CPLError(CE_Failure, CPLE_AppDefined,
1354 : "Unable to find layer \"%s\".",
1355 0 : sOptions.aosLayers.size() > static_cast<size_t>(i)
1356 0 : ? sOptions.aosLayers[i].c_str()
1357 : : "0");
1358 0 : eErr = CE_Failure;
1359 0 : break;
1360 : }
1361 :
1362 46 : if (!sOptions.osWHERE.empty())
1363 : {
1364 1 : if (OGR_L_SetAttributeFilter(hLayer, sOptions.osWHERE.c_str()) !=
1365 : OGRERR_NONE)
1366 : {
1367 0 : eErr = CE_Failure;
1368 0 : break;
1369 : }
1370 : }
1371 :
1372 92 : void *pScaledProgress = GDALCreateScaledProgress(
1373 46 : 0.0, 1.0 * (i + 1) / nLayerCount, sOptions.pfnProgress,
1374 : sOptions.pProgressData);
1375 :
1376 92 : eErr = ProcessLayer(
1377 46 : hLayer, !sOptions.oOutputSRS.IsEmpty(), hDstDS, sOptions.anBandList,
1378 46 : sOptions.adfBurnValues, sOptions.b3D, sOptions.bInverse,
1379 : sOptions.osBurnAttribute.c_str(), sOptions.aosRasterizeOptions,
1380 46 : sOptions.aosTO, GDALScaledProgress, pScaledProgress);
1381 :
1382 46 : GDALDestroyScaledProgress(pScaledProgress);
1383 46 : if (eErr != CE_None)
1384 1 : break;
1385 : }
1386 :
1387 54 : if (eErr != CE_None)
1388 : {
1389 1 : if (bCloseOutDSOnError)
1390 0 : GDALClose(hDstDS);
1391 1 : return nullptr;
1392 : }
1393 :
1394 53 : return hDstDS;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* ArgIsNumericRasterize() */
1399 : /************************************************************************/
1400 :
1401 401 : static bool ArgIsNumericRasterize(const char *pszArg)
1402 :
1403 : {
1404 401 : char *pszEnd = nullptr;
1405 401 : CPLStrtod(pszArg, &pszEnd);
1406 401 : return pszEnd != nullptr && pszEnd[0] == '\0';
1407 : }
1408 :
1409 : /************************************************************************/
1410 : /* GDALRasterizeOptionsNew() */
1411 : /************************************************************************/
1412 :
1413 : /**
1414 : * Allocates a GDALRasterizeOptions struct.
1415 : *
1416 : * @param papszArgv NULL terminated list of options (potentially including
1417 : * filename and open options too), or NULL. The accepted options are the ones of
1418 : * the <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
1419 : * @param psOptionsForBinary (output) may be NULL (and should generally be
1420 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
1421 : * GDALRasterizeOptionsForBinaryNew() prior to this
1422 : * function. Will be filled with potentially present filename, open options,...
1423 : * @return pointer to the allocated GDALRasterizeOptions struct. Must be freed
1424 : * with GDALRasterizeOptionsFree().
1425 : *
1426 : * @since GDAL 2.1
1427 : */
1428 :
1429 : GDALRasterizeOptions *
1430 66 : GDALRasterizeOptionsNew(char **papszArgv,
1431 : GDALRasterizeOptionsForBinary *psOptionsForBinary)
1432 : {
1433 :
1434 132 : auto psOptions = std::make_unique<GDALRasterizeOptions>();
1435 :
1436 : /*-------------------------------------------------------------------- */
1437 : /* Parse arguments. */
1438 : /*-------------------------------------------------------------------- */
1439 :
1440 132 : CPLStringList aosArgv;
1441 :
1442 : /* -------------------------------------------------------------------- */
1443 : /* Pre-processing for custom syntax that ArgumentParser does not */
1444 : /* support. */
1445 : /* -------------------------------------------------------------------- */
1446 66 : const int argc = CSLCount(papszArgv);
1447 674 : for (int i = 0; i < argc && papszArgv != nullptr && papszArgv[i] != nullptr;
1448 : i++)
1449 : {
1450 : // argparser will be confused if the value of a string argument
1451 : // starts with a negative sign.
1452 608 : if (EQUAL(papszArgv[i], "-a_nodata") && papszArgv[i + 1])
1453 : {
1454 5 : ++i;
1455 5 : psOptions->osNoData = papszArgv[i];
1456 5 : psOptions->bCreateOutput = true;
1457 : }
1458 :
1459 : // argparser is confused by arguments that have at_least_one
1460 : // cardinality, if they immediately precede positional arguments.
1461 603 : else if (EQUAL(papszArgv[i], "-burn") && papszArgv[i + 1])
1462 : {
1463 112 : if (strchr(papszArgv[i + 1], ' '))
1464 : {
1465 : const CPLStringList aosTokens(
1466 0 : CSLTokenizeString(papszArgv[i + 1]));
1467 0 : for (const char *pszToken : aosTokens)
1468 : {
1469 0 : psOptions->adfBurnValues.push_back(CPLAtof(pszToken));
1470 : }
1471 0 : i += 1;
1472 : }
1473 : else
1474 : {
1475 224 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1476 : {
1477 112 : psOptions->adfBurnValues.push_back(
1478 112 : CPLAtof(papszArgv[i + 1]));
1479 112 : i += 1;
1480 : }
1481 : }
1482 :
1483 : // Dummy value to make argparse happy, as at least one of
1484 : // -burn, -a or -3d is required
1485 112 : aosArgv.AddString("-burn");
1486 112 : aosArgv.AddString("0");
1487 : }
1488 491 : else if (EQUAL(papszArgv[i], "-init") && papszArgv[i + 1])
1489 : {
1490 27 : if (strchr(papszArgv[i + 1], ' '))
1491 : {
1492 : const CPLStringList aosTokens(
1493 0 : CSLTokenizeString(papszArgv[i + 1]));
1494 0 : for (const char *pszToken : aosTokens)
1495 : {
1496 0 : psOptions->adfInitVals.push_back(CPLAtof(pszToken));
1497 : }
1498 0 : i += 1;
1499 : }
1500 : else
1501 : {
1502 54 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1503 : {
1504 27 : psOptions->adfInitVals.push_back(CPLAtof(papszArgv[i + 1]));
1505 27 : i += 1;
1506 : }
1507 : }
1508 27 : psOptions->bCreateOutput = true;
1509 : }
1510 464 : else if (EQUAL(papszArgv[i], "-b") && papszArgv[i + 1])
1511 : {
1512 66 : if (strchr(papszArgv[i + 1], ' '))
1513 : {
1514 : const CPLStringList aosTokens(
1515 0 : CSLTokenizeString(papszArgv[i + 1]));
1516 0 : for (const char *pszToken : aosTokens)
1517 : {
1518 0 : psOptions->anBandList.push_back(atoi(pszToken));
1519 : }
1520 0 : i += 1;
1521 : }
1522 : else
1523 : {
1524 132 : while (i < argc - 1 && ArgIsNumericRasterize(papszArgv[i + 1]))
1525 : {
1526 66 : psOptions->anBandList.push_back(atoi(papszArgv[i + 1]));
1527 66 : i += 1;
1528 : }
1529 66 : }
1530 : }
1531 : else
1532 : {
1533 398 : aosArgv.AddString(papszArgv[i]);
1534 : }
1535 : }
1536 :
1537 : try
1538 : {
1539 : auto argParser =
1540 68 : GDALRasterizeOptionsGetParser(psOptions.get(), psOptionsForBinary);
1541 66 : argParser->parse_args_without_binary_name(aosArgv.List());
1542 :
1543 : // Check all no store_into args
1544 67 : if (auto oTe = argParser->present<std::vector<double>>("-te"))
1545 : {
1546 3 : psOptions->sEnvelop.MinX = oTe.value()[0];
1547 3 : psOptions->sEnvelop.MinY = oTe.value()[1];
1548 3 : psOptions->sEnvelop.MaxX = oTe.value()[2];
1549 3 : psOptions->sEnvelop.MaxY = oTe.value()[3];
1550 3 : psOptions->bCreateOutput = true;
1551 : }
1552 :
1553 64 : if (auto oTr = argParser->present<std::vector<double>>("-tr"))
1554 : {
1555 13 : psOptions->dfXRes = oTr.value()[0];
1556 13 : psOptions->dfYRes = oTr.value()[1];
1557 :
1558 13 : if (psOptions->dfXRes <= 0 || psOptions->dfYRes <= 0)
1559 : {
1560 0 : CPLError(CE_Failure, CPLE_AppDefined,
1561 : "Wrong value for -tr parameter.");
1562 0 : return nullptr;
1563 : }
1564 :
1565 13 : psOptions->bCreateOutput = true;
1566 : }
1567 :
1568 64 : if (auto oTs = argParser->present<std::vector<double>>("-ts"))
1569 : {
1570 21 : const int nXSize = static_cast<int>(oTs.value()[0]);
1571 21 : const int nYSize = static_cast<int>(oTs.value()[1]);
1572 :
1573 : // Warn the user if the conversion to int looses precision
1574 21 : if (nXSize != oTs.value()[0] || nYSize != oTs.value()[1])
1575 : {
1576 1 : CPLError(CE_Warning, CPLE_AppDefined,
1577 : "-ts values parsed as %d %d.", nXSize, nYSize);
1578 : }
1579 :
1580 21 : psOptions->nXSize = nXSize;
1581 21 : psOptions->nYSize = nYSize;
1582 :
1583 21 : if (!(psOptions->nXSize > 0 || psOptions->nYSize > 0))
1584 : {
1585 0 : CPLError(CE_Failure, CPLE_AppDefined,
1586 : "Wrong value for -ts parameter: at least one of the "
1587 : "arguments must be greater than zero.");
1588 0 : return nullptr;
1589 : }
1590 :
1591 21 : psOptions->bCreateOutput = true;
1592 : }
1593 :
1594 64 : if (psOptions->bCreateOutput)
1595 : {
1596 55 : if (psOptions->dfXRes == 0 && psOptions->dfYRes == 0 &&
1597 55 : psOptions->nXSize == 0 && psOptions->nYSize == 0)
1598 : {
1599 0 : CPLError(CE_Failure, CPLE_NotSupported,
1600 : "'-tr xres yres' or '-ts xsize ysize' is required.");
1601 0 : return nullptr;
1602 : }
1603 :
1604 34 : if (psOptions->bTargetAlignedPixels && psOptions->dfXRes == 0 &&
1605 0 : psOptions->dfYRes == 0)
1606 : {
1607 0 : CPLError(CE_Failure, CPLE_NotSupported,
1608 : "-tap option cannot be used without using -tr.");
1609 0 : return nullptr;
1610 : }
1611 :
1612 34 : if (!psOptions->anBandList.empty())
1613 : {
1614 1 : CPLError(
1615 : CE_Failure, CPLE_NotSupported,
1616 : "-b option cannot be used when creating a GDAL dataset.");
1617 1 : return nullptr;
1618 : }
1619 :
1620 33 : int nBandCount = 1;
1621 :
1622 33 : if (!psOptions->adfBurnValues.empty())
1623 19 : nBandCount = static_cast<int>(psOptions->adfBurnValues.size());
1624 :
1625 33 : if (static_cast<int>(psOptions->adfInitVals.size()) > nBandCount)
1626 0 : nBandCount = static_cast<int>(psOptions->adfInitVals.size());
1627 :
1628 33 : if (psOptions->adfInitVals.size() == 1)
1629 : {
1630 3 : for (int i = 1; i <= nBandCount - 1; i++)
1631 0 : psOptions->adfInitVals.push_back(psOptions->adfInitVals[0]);
1632 : }
1633 :
1634 94 : for (int i = 1; i <= nBandCount; i++)
1635 61 : psOptions->anBandList.push_back(i);
1636 : }
1637 : else
1638 : {
1639 30 : if (psOptions->anBandList.empty())
1640 11 : psOptions->anBandList.push_back(1);
1641 : }
1642 :
1643 63 : if (!psOptions->osDialect.empty() && !psOptions->osWHERE.empty() &&
1644 0 : !psOptions->osSQL.empty())
1645 : {
1646 0 : CPLError(CE_Warning, CPLE_AppDefined,
1647 : "-dialect is ignored with -where. Use -sql instead");
1648 : }
1649 :
1650 63 : if (psOptionsForBinary)
1651 : {
1652 11 : psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
1653 11 : if (!psOptions->osFormat.empty())
1654 0 : psOptionsForBinary->osFormat = psOptions->osFormat;
1655 : }
1656 68 : else if (psOptions->adfBurnValues.empty() &&
1657 68 : psOptions->osBurnAttribute.empty() && !psOptions->b3D)
1658 : {
1659 11 : psOptions->adfBurnValues.push_back(255);
1660 : }
1661 : }
1662 2 : catch (const std::exception &e)
1663 : {
1664 2 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
1665 2 : return nullptr;
1666 : }
1667 :
1668 63 : return psOptions.release();
1669 : }
1670 :
1671 : /************************************************************************/
1672 : /* GDALRasterizeOptionsFree() */
1673 : /************************************************************************/
1674 :
1675 : /**
1676 : * Frees the GDALRasterizeOptions struct.
1677 : *
1678 : * @param psOptions the options struct for GDALRasterize().
1679 : *
1680 : * @since GDAL 2.1
1681 : */
1682 :
1683 63 : void GDALRasterizeOptionsFree(GDALRasterizeOptions *psOptions)
1684 : {
1685 63 : delete psOptions;
1686 63 : }
1687 :
1688 : /************************************************************************/
1689 : /* GDALRasterizeOptionsSetProgress() */
1690 : /************************************************************************/
1691 :
1692 : /**
1693 : * Set a progress function.
1694 : *
1695 : * @param psOptions the options struct for GDALRasterize().
1696 : * @param pfnProgress the progress callback.
1697 : * @param pProgressData the user data for the progress callback.
1698 : *
1699 : * @since GDAL 2.1
1700 : */
1701 :
1702 39 : void GDALRasterizeOptionsSetProgress(GDALRasterizeOptions *psOptions,
1703 : GDALProgressFunc pfnProgress,
1704 : void *pProgressData)
1705 : {
1706 39 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
1707 39 : psOptions->pProgressData = pProgressData;
1708 39 : }
|
