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