Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Contour Generation
4 : * Purpose: Core algorithm implementation for contour line generation.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2003, Applied Coherent Technology Corporation, www.actgate.com
10 : * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
11 : * Copyright (c) 2018, Oslandia <infos at oslandia dot com>
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include "level_generator.h"
17 : #include "polygon_ring_appender.h"
18 : #include "utility.h"
19 : #include "contour_generator.h"
20 : #include "segment_merger.h"
21 : #include <algorithm>
22 :
23 : #include "gdal.h"
24 : #include "gdal_alg.h"
25 : #include "cpl_conv.h"
26 : #include "cpl_error_internal.h"
27 : #include "cpl_string.h"
28 : #include "ogr_api.h"
29 : #include "ogr_srs_api.h"
30 : #include "ogr_geometry.h"
31 :
32 : #include <climits>
33 : #include <limits>
34 :
35 96 : static CPLErr OGRPolygonContourWriter(double dfLevelMin, double dfLevelMax,
36 : const OGRMultiPolygon &multipoly,
37 : void *pInfo)
38 :
39 : {
40 96 : OGRContourWriterInfo *poInfo = static_cast<OGRContourWriterInfo *>(pInfo);
41 :
42 : OGRFeatureDefnH hFDefn =
43 96 : OGR_L_GetLayerDefn(static_cast<OGRLayerH>(poInfo->hLayer));
44 :
45 96 : OGRFeatureH hFeat = OGR_F_Create(hFDefn);
46 :
47 96 : if (poInfo->nIDField != -1)
48 96 : OGR_F_SetFieldInteger(hFeat, poInfo->nIDField, poInfo->nNextID++);
49 :
50 96 : if (poInfo->nElevFieldMin != -1)
51 96 : OGR_F_SetFieldDouble(hFeat, poInfo->nElevFieldMin, dfLevelMin);
52 :
53 96 : if (poInfo->nElevFieldMax != -1)
54 96 : OGR_F_SetFieldDouble(hFeat, poInfo->nElevFieldMax, dfLevelMax);
55 :
56 96 : const bool bHasZ = wkbHasZ(OGR_FD_GetGeomType(hFDefn));
57 : OGRGeometryH hGeom =
58 96 : OGR_G_CreateGeometry(bHasZ ? wkbMultiPolygon25D : wkbMultiPolygon);
59 :
60 444 : for (int iPart = 0; iPart < multipoly.getNumGeometries(); iPart++)
61 : {
62 348 : OGRPolygon *poNewPoly = new OGRPolygon();
63 : const OGRPolygon *poPolygon =
64 348 : static_cast<const OGRPolygon *>(multipoly.getGeometryRef(iPart));
65 :
66 974 : for (int iRing = 0; iRing < poPolygon->getNumInteriorRings() + 1;
67 : iRing++)
68 : {
69 : const OGRLinearRing *poRing =
70 626 : iRing == 0 ? poPolygon->getExteriorRing()
71 278 : : poPolygon->getInteriorRing(iRing - 1);
72 :
73 626 : OGRLinearRing *poNewRing = new OGRLinearRing();
74 53403 : for (int iPoint = 0; iPoint < poRing->getNumPoints(); iPoint++)
75 : {
76 : const double dfX =
77 105554 : poInfo->adfGeoTransform[0] +
78 52777 : poInfo->adfGeoTransform[1] * poRing->getX(iPoint) +
79 52777 : poInfo->adfGeoTransform[2] * poRing->getY(iPoint);
80 : const double dfY =
81 105554 : poInfo->adfGeoTransform[3] +
82 52777 : poInfo->adfGeoTransform[4] * poRing->getX(iPoint) +
83 52777 : poInfo->adfGeoTransform[5] * poRing->getY(iPoint);
84 52777 : if (bHasZ)
85 0 : OGR_G_SetPoint(OGRGeometry::ToHandle(poNewRing), iPoint,
86 : dfX, dfY, dfLevelMax);
87 : else
88 52777 : OGR_G_SetPoint_2D(OGRGeometry::ToHandle(poNewRing), iPoint,
89 : dfX, dfY);
90 : }
91 626 : poNewPoly->addRingDirectly(poNewRing);
92 : }
93 348 : OGR_G_AddGeometryDirectly(hGeom, OGRGeometry::ToHandle(poNewPoly));
94 : }
95 :
96 96 : OGR_F_SetGeometryDirectly(hFeat, hGeom);
97 :
98 : OGRErr eErr =
99 96 : OGR_L_CreateFeature(static_cast<OGRLayerH>(poInfo->hLayer), hFeat);
100 96 : OGR_F_Destroy(hFeat);
101 :
102 96 : if (eErr == OGRERR_NONE && poInfo->nTransactionCommitInterval > 0)
103 : {
104 12 : if (++poInfo->nWrittenFeatureCountSinceLastCommit ==
105 12 : poInfo->nTransactionCommitInterval)
106 : {
107 6 : poInfo->nWrittenFeatureCountSinceLastCommit = 0;
108 : // CPLDebug("CONTOUR", "Flush transaction");
109 : eErr =
110 6 : OGR_L_CommitTransaction(static_cast<OGRLayerH>(poInfo->hLayer));
111 6 : if (eErr == OGRERR_NONE)
112 : {
113 6 : eErr = OGR_L_StartTransaction(
114 6 : static_cast<OGRLayerH>(poInfo->hLayer));
115 : }
116 : }
117 : }
118 :
119 96 : return eErr == OGRERR_NONE ? CE_None : CE_Failure;
120 : }
121 :
122 : struct PolygonContourWriter
123 : {
124 : CPL_DISALLOW_COPY_ASSIGN(PolygonContourWriter)
125 :
126 29 : explicit PolygonContourWriter(OGRContourWriterInfo *poInfo, double minLevel)
127 29 : : poInfo_(poInfo), currentLevel_(minLevel)
128 : {
129 29 : }
130 :
131 112 : void startPolygon(double level)
132 : {
133 112 : previousLevel_ = currentLevel_;
134 112 : currentGeometry_.reset(new OGRMultiPolygon());
135 112 : currentLevel_ = level;
136 112 : }
137 :
138 112 : void endPolygon()
139 : {
140 112 : if (currentPart_)
141 : {
142 96 : currentGeometry_->addGeometryDirectly(currentPart_);
143 : }
144 :
145 112 : if (currentGeometry_->getNumGeometries() > 0)
146 : {
147 96 : OGRPolygonContourWriter(previousLevel_, currentLevel_,
148 96 : *currentGeometry_, poInfo_);
149 : }
150 :
151 112 : currentGeometry_.reset(nullptr);
152 112 : currentPart_ = nullptr;
153 112 : }
154 :
155 348 : void addPart(const marching_squares::LineString &ring)
156 : {
157 348 : if (currentPart_)
158 : {
159 252 : currentGeometry_->addGeometryDirectly(currentPart_);
160 : }
161 :
162 348 : OGRLinearRing *poNewRing = new OGRLinearRing();
163 348 : poNewRing->setNumPoints(int(ring.size()));
164 348 : int iPoint = 0;
165 43631 : for (const auto &p : ring)
166 : {
167 43283 : poNewRing->setPoint(iPoint, p.x, p.y);
168 43283 : iPoint++;
169 : }
170 348 : currentPart_ = new OGRPolygon();
171 348 : currentPart_->addRingDirectly(poNewRing);
172 348 : }
173 :
174 278 : void addInteriorRing(const marching_squares::LineString &ring)
175 : {
176 278 : OGRLinearRing *poNewRing = new OGRLinearRing();
177 9772 : for (const auto &p : ring)
178 : {
179 9494 : poNewRing->addPoint(p.x, p.y);
180 : }
181 278 : currentPart_->addRingDirectly(poNewRing);
182 278 : }
183 :
184 : std::unique_ptr<OGRMultiPolygon> currentGeometry_ = {};
185 : OGRPolygon *currentPart_ = nullptr;
186 : OGRContourWriterInfo *poInfo_ = nullptr;
187 : double currentLevel_ = 0;
188 : double previousLevel_ = 0;
189 : };
190 :
191 : struct GDALRingAppender
192 : {
193 : CPL_DISALLOW_COPY_ASSIGN(GDALRingAppender)
194 :
195 46 : GDALRingAppender(GDALContourWriter write, void *data)
196 46 : : write_(write), data_(data)
197 : {
198 46 : }
199 :
200 4717 : void addLine(double level, marching_squares::LineString &ls,
201 : bool /*closed*/)
202 : {
203 4717 : const size_t sz = ls.size();
204 14151 : std::vector<double> xs(sz), ys(sz);
205 4717 : size_t i = 0;
206 85713 : for (const auto &pt : ls)
207 : {
208 80996 : xs[i] = pt.x;
209 80996 : ys[i] = pt.y;
210 80996 : i++;
211 : }
212 :
213 4717 : if (write_(level, int(sz), &xs[0], &ys[0], data_) != CE_None)
214 0 : CPLError(CE_Failure, CPLE_AppDefined, "cannot write linestring");
215 4717 : }
216 :
217 : private:
218 : GDALContourWriter write_;
219 : void *data_;
220 : };
221 :
222 : /************************************************************************/
223 : /* ==================================================================== */
224 : /* Additional C Callable Functions */
225 : /* ==================================================================== */
226 : /************************************************************************/
227 :
228 : /************************************************************************/
229 : /* OGRContourWriter() */
230 : /************************************************************************/
231 :
232 4717 : CPLErr OGRContourWriter(double dfLevel, int nPoints, double *padfX,
233 : double *padfY, void *pInfo)
234 :
235 : {
236 4717 : OGRContourWriterInfo *poInfo = static_cast<OGRContourWriterInfo *>(pInfo);
237 :
238 : OGRFeatureDefnH hFDefn =
239 4717 : OGR_L_GetLayerDefn(static_cast<OGRLayerH>(poInfo->hLayer));
240 :
241 4717 : OGRFeatureH hFeat = OGR_F_Create(hFDefn);
242 :
243 4717 : if (poInfo->nIDField != -1)
244 4717 : OGR_F_SetFieldInteger(hFeat, poInfo->nIDField, poInfo->nNextID++);
245 :
246 4717 : if (poInfo->nElevField != -1)
247 301 : OGR_F_SetFieldDouble(hFeat, poInfo->nElevField, dfLevel);
248 :
249 4717 : const bool bHasZ = wkbHasZ(OGR_FD_GetGeomType(hFDefn));
250 : OGRGeometryH hGeom =
251 4717 : OGR_G_CreateGeometry(bHasZ ? wkbLineString25D : wkbLineString);
252 :
253 85713 : for (int iPoint = nPoints - 1; iPoint >= 0; iPoint--)
254 : {
255 80996 : const double dfX = poInfo->adfGeoTransform[0] +
256 80996 : poInfo->adfGeoTransform[1] * padfX[iPoint] +
257 80996 : poInfo->adfGeoTransform[2] * padfY[iPoint];
258 80996 : const double dfY = poInfo->adfGeoTransform[3] +
259 80996 : poInfo->adfGeoTransform[4] * padfX[iPoint] +
260 80996 : poInfo->adfGeoTransform[5] * padfY[iPoint];
261 80996 : if (bHasZ)
262 49759 : OGR_G_SetPoint(hGeom, iPoint, dfX, dfY, dfLevel);
263 : else
264 31237 : OGR_G_SetPoint_2D(hGeom, iPoint, dfX, dfY);
265 : }
266 :
267 4717 : OGR_F_SetGeometryDirectly(hFeat, hGeom);
268 :
269 : const OGRErr eErr =
270 4717 : OGR_L_CreateFeature(static_cast<OGRLayerH>(poInfo->hLayer), hFeat);
271 4717 : OGR_F_Destroy(hFeat);
272 :
273 4717 : return eErr == OGRERR_NONE ? CE_None : CE_Failure;
274 : }
275 :
276 : /************************************************************************/
277 : /* GDALContourGenerate() */
278 : /************************************************************************/
279 :
280 : /**
281 : * Create vector contours from raster DEM.
282 : *
283 : * This function is kept for compatibility reason and will call the new
284 : * variant GDALContourGenerateEx that is more extensible and provide more
285 : * options.
286 : *
287 : * Details about the algorithm are also given in the documentation of the
288 : * new GDALContourenerateEx function.
289 : *
290 : * @param hBand The band to read raster data from. The whole band will be
291 : * processed.
292 : *
293 : * @param dfContourInterval The elevation interval between contours generated.
294 : *
295 : * @param dfContourBase The "base" relative to which contour intervals are
296 : * applied. This is normally zero, but could be different. To generate 10m
297 : * contours at 5, 15, 25, ... the ContourBase would be 5.
298 : *
299 : * @param nFixedLevelCount The number of fixed levels. If this is greater than
300 : * zero, then fixed levels will be used, and ContourInterval and ContourBase
301 : * are ignored.
302 : *
303 : * @param padfFixedLevels The list of fixed contour levels at which contours
304 : * should be generated. It will contain FixedLevelCount entries, and may be
305 : * NULL if fixed levels are disabled (FixedLevelCount = 0).
306 : *
307 : * @param bUseNoData If TRUE the dfNoDataValue will be used.
308 : *
309 : * @param dfNoDataValue The value to use as a "nodata" value. That is, a
310 : * pixel value which should be ignored in generating contours as if the value
311 : * of the pixel were not known.
312 : *
313 : * @param hLayer The layer to which new contour vectors will be written.
314 : * Each contour will have a LINESTRING geometry attached to it. This
315 : * is really of type OGRLayerH, but void * is used to avoid pulling the
316 : * ogr_api.h file in here.
317 : *
318 : * @param iIDField If not -1 this will be used as a field index to indicate
319 : * where a unique id should be written for each feature (contour) written.
320 : *
321 : * @param iElevField If not -1 this will be used as a field index to indicate
322 : * where the elevation value of the contour should be written.
323 : *
324 : * @param pfnProgress A GDALProgressFunc that may be used to report progress
325 : * to the user, or to interrupt the algorithm. May be NULL if not required.
326 : *
327 : * @param pProgressArg The callback data for the pfnProgress function.
328 : *
329 : * @return CE_None on success or CE_Failure if an error occurs.
330 : */
331 :
332 4 : CPLErr GDALContourGenerate(GDALRasterBandH hBand, double dfContourInterval,
333 : double dfContourBase, int nFixedLevelCount,
334 : double *padfFixedLevels, int bUseNoData,
335 : double dfNoDataValue, void *hLayer, int iIDField,
336 : int iElevField, GDALProgressFunc pfnProgress,
337 : void *pProgressArg)
338 : {
339 4 : char **options = nullptr;
340 4 : if (nFixedLevelCount > 0)
341 : {
342 1 : std::string values = "FIXED_LEVELS=";
343 4 : for (int i = 0; i < nFixedLevelCount; i++)
344 : {
345 3 : const int sz = 32;
346 3 : char *newValue = new char[sz + 1];
347 3 : if (i == nFixedLevelCount - 1)
348 : {
349 1 : CPLsnprintf(newValue, sz + 1, "%f", padfFixedLevels[i]);
350 : }
351 : else
352 : {
353 2 : CPLsnprintf(newValue, sz + 1, "%f,", padfFixedLevels[i]);
354 : }
355 3 : values = values + std::string(newValue);
356 3 : delete[] newValue;
357 : }
358 1 : options = CSLAddString(options, values.c_str());
359 : }
360 3 : else if (dfContourInterval != 0.0)
361 : {
362 : options =
363 3 : CSLAppendPrintf(options, "LEVEL_INTERVAL=%f", dfContourInterval);
364 : }
365 :
366 4 : if (dfContourBase != 0.0)
367 : {
368 0 : options = CSLAppendPrintf(options, "LEVEL_BASE=%f", dfContourBase);
369 : }
370 :
371 4 : if (bUseNoData)
372 : {
373 0 : options = CSLAppendPrintf(options, "NODATA=%.19g", dfNoDataValue);
374 : }
375 4 : if (iIDField != -1)
376 : {
377 4 : options = CSLAppendPrintf(options, "ID_FIELD=%d", iIDField);
378 : }
379 4 : if (iElevField != -1)
380 : {
381 3 : options = CSLAppendPrintf(options, "ELEV_FIELD=%d", iElevField);
382 : }
383 :
384 4 : CPLErr err = GDALContourGenerateEx(hBand, hLayer, options, pfnProgress,
385 : pProgressArg);
386 4 : CSLDestroy(options);
387 :
388 4 : return err;
389 : }
390 :
391 : /**
392 : * Create vector contours from raster DEM.
393 : *
394 : * This algorithm is an implementation of "Marching squares" [1] that will
395 : * generate contour vectors for the input raster band on the requested set
396 : * of contour levels. The vector contours are written to the passed in OGR
397 : * vector layer. Also, a NODATA value may be specified to identify pixels
398 : * that should not be considered in contour line generation.
399 : *
400 : * The gdal/apps/gdal_contour_bin.cpp mainline can be used as an example of
401 : * how to use this function.
402 : *
403 : * [1] see https://en.wikipedia.org/wiki/Marching_squares
404 : *
405 : * ALGORITHM RULES
406 :
407 : For contouring purposes raster pixel values are assumed to represent a point
408 : value at the center of the corresponding pixel region. For the purpose of
409 : contour generation we virtually connect each pixel center to the values to
410 : the left, right, top and bottom. We assume that the pixel value is linearly
411 : interpolated between the pixel centers along each line, and determine where
412 : (if any) contour lines will appear along these line segments. Then the
413 : contour crossings are connected.
414 :
415 : This means that contour lines' nodes will not actually be on pixel edges, but
416 : rather along vertical and horizontal lines connecting the pixel centers.
417 :
418 : \verbatim
419 : General Case:
420 :
421 : 5 | | 3
422 : -- + ---------------- + --
423 : | |
424 : | |
425 : | |
426 : | |
427 : 10 + |
428 : |\ |
429 : | \ |
430 : -- + -+-------------- + --
431 : 12 | 10 | 1
432 :
433 : Saddle Point:
434 :
435 : 5 | | 12
436 : -- + -------------+-- + --
437 : | \ |
438 : | \|
439 : | +
440 : | |
441 : + |
442 : |\ |
443 : | \ |
444 : -- + -+-------------- + --
445 : 12 | | 1
446 :
447 : or:
448 :
449 : 5 | | 12
450 : -- + -------------+-- + --
451 : | __/ |
452 : | ___/ |
453 : | ___/ __+
454 : | / __/ |
455 : +' __/ |
456 : | __/ |
457 : | ,__/ |
458 : -- + -+-------------- + --
459 : 12 | | 1
460 : \endverbatim
461 :
462 : Nodata:
463 :
464 : In the "nodata" case we treat the whole nodata pixel as a no-mans land.
465 : We extend the corner pixels near the nodata out to half way and then
466 : construct extra lines from those points to the center which is assigned
467 : an averaged value from the two nearby points (in this case (12+3+5)/3).
468 :
469 : \verbatim
470 : 5 | | 3
471 : -- + ---------------- + --
472 : | |
473 : | |
474 : | 6.7 |
475 : | +---------+ 3
476 : 10 +___ |
477 : | \____+ 10
478 : | |
479 : -- + -------+ +
480 : 12 | 12 (nodata)
481 :
482 : \endverbatim
483 :
484 : *
485 : * @param hBand The band to read raster data from. The whole band will be
486 : * processed.
487 : *
488 : * @param hLayer The layer to which new contour vectors will be written.
489 : * Each contour will have a LINESTRING geometry attached to it
490 : * (or POLYGON if POLYGONIZE=YES). This is really of type OGRLayerH, but
491 : * void * is used to avoid pulling the ogr_api.h file in here.
492 : *
493 : * @param pfnProgress A GDALProgressFunc that may be used to report progress
494 : * to the user, or to interrupt the algorithm. May be NULL if not required.
495 : *
496 : * @param pProgressArg The callback data for the pfnProgress function.
497 : *
498 : * @param options List of options
499 : *
500 : * Options:
501 : *
502 : * LEVEL_INTERVAL=f
503 : *
504 : * The elevation interval between contours generated.
505 : *
506 : * LEVEL_BASE=f
507 : *
508 : * The "base" relative to which contour intervals are
509 : * applied. This is normally zero, but could be different. To generate 10m
510 : * contours at 5, 15, 25, ... the LEVEL_BASE would be 5.
511 : *
512 : * LEVEL_EXP_BASE=f
513 : *
514 : * If greater than 0, contour levels are generated on an
515 : * exponential scale. Levels will then be generated by LEVEL_EXP_BASE^k
516 : * where k is a positive integer.
517 : *
518 : * FIXED_LEVELS=f[,f]*
519 : *
520 : * The list of fixed contour levels at which contours should be generated.
521 : * This option has precedence on LEVEL_INTERVAL. MIN and MAX can be used
522 : * as special values to represent the minimum and maximum values of the
523 : * raster.
524 : *
525 : * NODATA=f
526 : *
527 : * The value to use as a "nodata" value. That is, a pixel value which
528 : * should be ignored in generating contours as if the value of the pixel
529 : * were not known.
530 : *
531 : * ID_FIELD=d
532 : *
533 : * This will be used as a field index to indicate where a unique id should
534 : * be written for each feature (contour) written.
535 : *
536 : * ELEV_FIELD=d
537 : *
538 : * This will be used as a field index to indicate where the elevation value
539 : * of the contour should be written. Only used in line contouring mode.
540 : *
541 : * ELEV_FIELD_MIN=d
542 : *
543 : * This will be used as a field index to indicate where the minimum elevation
544 : value
545 : * of the polygon contour should be written. Only used in polygonal contouring
546 : mode.
547 : *
548 : * ELEV_FIELD_MAX=d
549 : *
550 : * This will be used as a field index to indicate where the maximum elevation
551 : value
552 : * of the polygon contour should be written. Only used in polygonal contouring
553 : mode.
554 : *
555 : * POLYGONIZE=YES|NO
556 : *
557 : * If YES, contour polygons will be created, rather than polygon lines.
558 : *
559 : *
560 : * COMMIT_INTERVAL=num
561 : *
562 : * (GDAL >= 3.10) Interval in number of features at which transactions must be
563 : * flushed. The default value of 0 means that no transactions are opened.
564 : * A negative value means a single transaction. The function takes care of
565 : * issuing the starting transaction and committing the final one.
566 : *
567 : * @return CE_None on success or CE_Failure if an error occurs.
568 : */
569 77 : CPLErr GDALContourGenerateEx(GDALRasterBandH hBand, void *hLayer,
570 : CSLConstList options, GDALProgressFunc pfnProgress,
571 : void *pProgressArg)
572 : {
573 77 : VALIDATE_POINTER1(hBand, "GDALContourGenerateEx", CE_Failure);
574 :
575 77 : if (pfnProgress == nullptr)
576 57 : pfnProgress = GDALDummyProgress;
577 :
578 77 : double contourInterval = 0.0;
579 77 : const char *opt = CSLFetchNameValue(options, "LEVEL_INTERVAL");
580 77 : if (opt)
581 : {
582 48 : contourInterval = CPLAtof(opt);
583 : // Written this way to catch NaN as well.
584 48 : if (!(contourInterval > 0))
585 : {
586 1 : CPLError(CE_Failure, CPLE_AppDefined,
587 : "Invalid value for LEVEL_INTERVAL. Should be strictly "
588 : "positive.");
589 1 : return CE_Failure;
590 : }
591 : }
592 :
593 76 : double contourBase = 0.0;
594 76 : opt = CSLFetchNameValue(options, "LEVEL_BASE");
595 76 : if (opt)
596 : {
597 4 : contourBase = CPLAtof(opt);
598 : }
599 :
600 76 : double expBase = 0.0;
601 76 : opt = CSLFetchNameValue(options, "LEVEL_EXP_BASE");
602 76 : if (opt)
603 : {
604 8 : expBase = CPLAtof(opt);
605 : }
606 :
607 152 : std::vector<double> fixedLevels;
608 76 : opt = CSLFetchNameValue(options, "FIXED_LEVELS");
609 76 : if (opt)
610 : {
611 : const CPLStringList aosLevels(
612 33 : CSLTokenizeStringComplex(opt, ",", FALSE, FALSE));
613 33 : fixedLevels.resize(aosLevels.size());
614 130 : for (size_t i = 0; i < fixedLevels.size(); i++)
615 : {
616 : // Handle min/max values
617 97 : if (EQUAL(aosLevels[i], "MIN"))
618 : {
619 6 : fixedLevels[i] = std::numeric_limits<double>::lowest();
620 : }
621 91 : else if (EQUAL(aosLevels[i], "MAX"))
622 : {
623 4 : fixedLevels[i] = std::numeric_limits<double>::max();
624 : }
625 : else
626 : {
627 87 : fixedLevels[i] = CPLAtof(aosLevels[i]);
628 : }
629 97 : if (i > 0 && !(fixedLevels[i] >= fixedLevels[i - 1]))
630 : {
631 0 : CPLError(CE_Failure, CPLE_AppDefined,
632 : "FIXED_LEVELS should be strictly increasing");
633 0 : return CE_Failure;
634 : }
635 : }
636 : }
637 :
638 76 : bool useNoData = false;
639 76 : double noDataValue = 0.0;
640 76 : opt = CSLFetchNameValue(options, "NODATA");
641 76 : if (opt)
642 : {
643 13 : useNoData = true;
644 13 : noDataValue = CPLAtof(opt);
645 13 : if (GDALGetRasterDataType(hBand) == GDT_Float32)
646 : {
647 2 : int bClamped = FALSE;
648 2 : int bRounded = FALSE;
649 2 : noDataValue = GDALAdjustValueToDataType(GDT_Float32, noDataValue,
650 : &bClamped, &bRounded);
651 : }
652 : }
653 :
654 76 : int idField = -1;
655 76 : opt = CSLFetchNameValue(options, "ID_FIELD");
656 76 : if (opt)
657 : {
658 76 : idField = atoi(opt);
659 : }
660 :
661 76 : int elevField = -1;
662 76 : opt = CSLFetchNameValue(options, "ELEV_FIELD");
663 76 : if (opt)
664 : {
665 32 : elevField = atoi(opt);
666 : }
667 :
668 76 : int elevFieldMin = -1;
669 76 : opt = CSLFetchNameValue(options, "ELEV_FIELD_MIN");
670 76 : if (opt)
671 : {
672 29 : elevFieldMin = atoi(opt);
673 : }
674 :
675 76 : int elevFieldMax = -1;
676 76 : opt = CSLFetchNameValue(options, "ELEV_FIELD_MAX");
677 76 : if (opt)
678 : {
679 29 : elevFieldMax = atoi(opt);
680 : }
681 :
682 76 : bool polygonize = CPLFetchBool(options, "POLYGONIZE", false);
683 :
684 76 : int bSuccessMin = FALSE;
685 76 : double dfMinimum = GDALGetRasterMinimum(hBand, &bSuccessMin);
686 76 : int bSuccessMax = FALSE;
687 76 : double dfMaximum = GDALGetRasterMaximum(hBand, &bSuccessMax);
688 76 : if ((!bSuccessMin || !bSuccessMax))
689 : {
690 : double adfMinMax[2];
691 76 : CPLErrorAccumulator oAccumulator;
692 : bool bGotMinMax;
693 : {
694 152 : CPLErrorStateBackuper oBackuper;
695 76 : auto oAccumulatorContext = oAccumulator.InstallForCurrentScope();
696 76 : CPL_IGNORE_RET_VAL(oAccumulatorContext);
697 76 : bGotMinMax =
698 76 : (GDALComputeRasterMinMax(hBand, false, adfMinMax) == CE_None);
699 : }
700 76 : if (bGotMinMax)
701 : {
702 74 : dfMinimum = adfMinMax[0];
703 74 : dfMaximum = adfMinMax[1];
704 : }
705 : else
706 : {
707 2 : bool bNoValidPixels = false;
708 6 : for (const auto &sError : oAccumulator.GetErrors())
709 : {
710 4 : if (sError.msg.find("no valid pixels found in sampling") !=
711 : std::string::npos)
712 : {
713 1 : bNoValidPixels = true;
714 : }
715 3 : else if (sError.type == CE_Debug)
716 : {
717 0 : CPLDebug("GDAL", "%s", sError.msg.c_str());
718 : }
719 : else
720 : {
721 3 : CPLError(sError.type, sError.no, "%s", sError.msg.c_str());
722 : }
723 : }
724 2 : if (bNoValidPixels && oAccumulator.GetErrors().size() == 1)
725 : {
726 1 : CPLDebug("GDAL",
727 : "contour: no valid pixels found in input band");
728 1 : return CE_None;
729 : }
730 1 : return CE_Failure;
731 : }
732 : }
733 :
734 : using namespace marching_squares;
735 :
736 : OGRContourWriterInfo oCWI;
737 74 : oCWI.hLayer = static_cast<OGRLayerH>(hLayer);
738 74 : oCWI.nElevField = elevField;
739 74 : oCWI.nElevFieldMin = elevFieldMin;
740 74 : oCWI.nElevFieldMax = elevFieldMax;
741 74 : oCWI.nIDField = idField;
742 74 : oCWI.adfGeoTransform[0] = 0.0;
743 74 : oCWI.adfGeoTransform[1] = 1.0;
744 74 : oCWI.adfGeoTransform[2] = 0.0;
745 74 : oCWI.adfGeoTransform[3] = 0.0;
746 74 : oCWI.adfGeoTransform[4] = 0.0;
747 74 : oCWI.adfGeoTransform[5] = 1.0;
748 74 : GDALDatasetH hSrcDS = GDALGetBandDataset(hBand);
749 74 : if (hSrcDS != nullptr)
750 74 : GDALGetGeoTransform(hSrcDS, oCWI.adfGeoTransform);
751 74 : oCWI.nNextID = 0;
752 74 : oCWI.nWrittenFeatureCountSinceLastCommit = 0;
753 74 : oCWI.nTransactionCommitInterval =
754 74 : CPLAtoGIntBig(CSLFetchNameValueDef(options, "COMMIT_INTERVAL", "0"));
755 :
756 74 : if (oCWI.nTransactionCommitInterval < 0)
757 1 : oCWI.nTransactionCommitInterval = std::numeric_limits<GIntBig>::max();
758 74 : if (oCWI.nTransactionCommitInterval > 0)
759 : {
760 4 : if (OGR_L_StartTransaction(static_cast<OGRLayerH>(hLayer)) !=
761 : OGRERR_NONE)
762 : {
763 0 : oCWI.nTransactionCommitInterval = 0;
764 : }
765 : }
766 :
767 74 : bool ok = true;
768 :
769 : // Replace fixed levels min/max values with raster min/max values
770 74 : if (!fixedLevels.empty())
771 : {
772 33 : if (fixedLevels[0] == std::numeric_limits<double>::lowest())
773 : {
774 6 : fixedLevels[0] = dfMinimum;
775 : }
776 33 : if (fixedLevels.back() == std::numeric_limits<double>::max())
777 : {
778 4 : fixedLevels.back() = dfMaximum;
779 : }
780 : }
781 :
782 : try
783 : {
784 74 : if (polygonize)
785 : {
786 :
787 29 : if (!fixedLevels.empty())
788 : {
789 :
790 : // If the minimum raster value is larger than the first requested
791 : // level, select the requested level that is just below the
792 : // minimum raster value
793 23 : if (fixedLevels[0] < dfMinimum)
794 : {
795 8 : for (size_t i = 1; i < fixedLevels.size(); ++i)
796 : {
797 8 : if (fixedLevels[i] >= dfMinimum)
798 : {
799 6 : dfMinimum = fixedLevels[i - 1];
800 6 : break;
801 : }
802 : }
803 : }
804 :
805 : // Largest requested level (levels are sorted)
806 23 : const double maxLevel{fixedLevels.back()};
807 :
808 : // If the maximum raster value is smaller than the last requested
809 : // level, select the requested level that is just above the
810 : // maximum raster value
811 23 : if (maxLevel > dfMaximum)
812 : {
813 11 : for (size_t i = fixedLevels.size() - 1; i > 0; --i)
814 : {
815 10 : if (fixedLevels[i] <= dfMaximum)
816 : {
817 4 : dfMaximum = fixedLevels[i + 1];
818 4 : break;
819 : }
820 : }
821 : }
822 :
823 : // If the maximum raster value is equal to the last requested
824 : // level, add a small value to the maximum to avoid skipping the
825 : // last level polygons
826 23 : if (maxLevel == dfMaximum)
827 : {
828 10 : dfMaximum = std::nextafter(
829 : dfMaximum, std::numeric_limits<double>::infinity());
830 : }
831 : }
832 :
833 58 : PolygonContourWriter w(&oCWI, dfMinimum);
834 : typedef PolygonRingAppender<PolygonContourWriter> RingAppender;
835 58 : RingAppender appender(w);
836 :
837 29 : if (expBase > 0.0)
838 : {
839 : // Do not provide the actual minimum value to level iterator
840 : // in polygonal case, otherwise it can result in a polygon
841 : // with a degenerate min=max range.
842 : ExponentialLevelRangeIterator generator(
843 5 : expBase, -std::numeric_limits<double>::infinity());
844 5 : auto levelIt{generator.range(dfMinimum, dfMaximum)};
845 14 : for (auto i = levelIt.begin(); i != levelIt.end(); ++i)
846 : {
847 9 : const double level = (*i).second;
848 9 : fixedLevels.push_back(level);
849 : }
850 : // Append minimum value to fixed levels
851 5 : fixedLevels.push_back(dfMinimum);
852 : // Append maximum value to fixed levels
853 5 : fixedLevels.push_back(dfMaximum);
854 : }
855 24 : else if (contourInterval != 0)
856 : {
857 : // Do not provide the actual minimum value to level iterator
858 : // in polygonal case, otherwise it can result in a polygon
859 : // with a degenerate min=max range.
860 : IntervalLevelRangeIterator generator(
861 : contourBase, contourInterval,
862 7 : -std::numeric_limits<double>::infinity());
863 7 : auto levelIt{generator.range(dfMinimum, dfMaximum)};
864 32 : for (auto i = levelIt.begin(); i != levelIt.end(); ++i)
865 : {
866 25 : const double level = (*i).second;
867 25 : fixedLevels.push_back(level);
868 : }
869 : // Append minimum value to fixed levels
870 7 : fixedLevels.push_back(dfMinimum);
871 : // Append maximum value to fixed levels
872 7 : fixedLevels.push_back(dfMaximum);
873 : }
874 :
875 29 : if (!fixedLevels.empty())
876 : {
877 29 : std::sort(fixedLevels.begin(), fixedLevels.end());
878 : auto uniqueIt =
879 29 : std::unique(fixedLevels.begin(), fixedLevels.end());
880 29 : fixedLevels.erase(uniqueIt, fixedLevels.end());
881 : // Do not provide the actual minimum value to level iterator
882 : // in polygonal case, otherwise it can result in a polygon
883 : // with a degenerate min=max range.
884 : FixedLevelRangeIterator levels(
885 29 : &fixedLevels[0], fixedLevels.size(),
886 58 : -std::numeric_limits<double>::infinity(), dfMaximum);
887 : SegmentMerger<RingAppender, FixedLevelRangeIterator> writer(
888 58 : appender, levels, /* polygonize */ true);
889 58 : std::vector<int> aoiSkipLevels;
890 : // Skip first and last levels (min/max) in polygonal case
891 29 : aoiSkipLevels.push_back(0);
892 29 : aoiSkipLevels.push_back(static_cast<int>(levels.levelsCount()));
893 29 : writer.setSkipLevels(aoiSkipLevels);
894 : ContourGeneratorFromRaster<decltype(writer),
895 : FixedLevelRangeIterator>
896 29 : cg(hBand, useNoData, noDataValue, writer, levels);
897 29 : ok = cg.process(pfnProgress, pProgressArg);
898 : }
899 : }
900 : else
901 : {
902 45 : GDALRingAppender appender(OGRContourWriter, &oCWI);
903 :
904 : // Append all exp levels to fixed levels
905 45 : if (expBase > 0.0)
906 : {
907 3 : ExponentialLevelRangeIterator generator(expBase, dfMinimum);
908 3 : auto levelIt{generator.range(dfMinimum, dfMaximum)};
909 3 : for (auto i = levelIt.begin(); i != levelIt.end(); ++i)
910 : {
911 2 : const double level = (*i).second;
912 2 : fixedLevels.push_back(level);
913 : }
914 : }
915 42 : else if (contourInterval != 0)
916 : {
917 : IntervalLevelRangeIterator levels(contourBase, contourInterval,
918 36 : dfMinimum);
919 36 : auto levelIt{levels.range(dfMinimum, dfMaximum)};
920 458 : for (auto i = levelIt.begin(); i != levelIt.end(); ++i)
921 : {
922 424 : const double level = (*i).second;
923 424 : fixedLevels.push_back(level);
924 : }
925 : }
926 :
927 41 : if (!fixedLevels.empty())
928 : {
929 39 : std::sort(fixedLevels.begin(), fixedLevels.end());
930 : auto uniqueIt =
931 39 : std::unique(fixedLevels.begin(), fixedLevels.end());
932 39 : fixedLevels.erase(uniqueIt, fixedLevels.end());
933 : FixedLevelRangeIterator levels(
934 39 : &fixedLevels[0], fixedLevels.size(), dfMinimum, dfMaximum);
935 : SegmentMerger<GDALRingAppender, FixedLevelRangeIterator> writer(
936 78 : appender, levels, /* polygonize */ false);
937 : ContourGeneratorFromRaster<decltype(writer),
938 : FixedLevelRangeIterator>
939 39 : cg(hBand, useNoData, noDataValue, writer, levels);
940 39 : ok = cg.process(pfnProgress, pProgressArg);
941 : }
942 : }
943 : }
944 4 : catch (const std::exception &e)
945 : {
946 4 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
947 4 : return CE_Failure;
948 : }
949 :
950 70 : if (oCWI.nTransactionCommitInterval > 0)
951 : {
952 : // CPLDebug("CONTOUR", "Flush transaction");
953 4 : if (OGR_L_CommitTransaction(static_cast<OGRLayerH>(hLayer)) !=
954 : OGRERR_NONE)
955 : {
956 0 : ok = false;
957 : }
958 : }
959 :
960 70 : if (ok)
961 : {
962 70 : pfnProgress(1.0, "", pProgressArg);
963 : }
964 :
965 70 : return ok ? CE_None : CE_Failure;
966 : }
967 :
968 : /************************************************************************/
969 : /* GDAL_CG_Create() */
970 : /************************************************************************/
971 :
972 : namespace marching_squares
973 : {
974 :
975 : // Opaque type used by the C API
976 : struct ContourGeneratorOpaque
977 : {
978 : typedef SegmentMerger<GDALRingAppender, IntervalLevelRangeIterator>
979 : SegmentMergerT;
980 : typedef ContourGenerator<SegmentMergerT, IntervalLevelRangeIterator>
981 : ContourGeneratorT;
982 :
983 1 : ContourGeneratorOpaque(int nWidth, int nHeight, int bNoDataSet,
984 : double dfNoDataValue, double dfContourInterval,
985 : double dfContourBase, GDALContourWriter pfnWriter,
986 : void *pCBData)
987 1 : : levels(dfContourBase, dfContourInterval,
988 1 : -std::numeric_limits<double>::infinity()),
989 : writer(pfnWriter, pCBData),
990 1 : merger(writer, levels, /* polygonize */ false),
991 : contourGenerator(nWidth, nHeight, bNoDataSet != 0, dfNoDataValue,
992 1 : merger, levels)
993 : {
994 1 : }
995 :
996 : IntervalLevelRangeIterator levels;
997 : GDALRingAppender writer;
998 : SegmentMergerT merger;
999 : ContourGeneratorT contourGenerator;
1000 : };
1001 :
1002 : } // namespace marching_squares
1003 :
1004 : /** Create contour generator */
1005 1 : GDALContourGeneratorH GDAL_CG_Create(int nWidth, int nHeight, int bNoDataSet,
1006 : double dfNoDataValue,
1007 : double dfContourInterval,
1008 : double dfContourBase,
1009 : GDALContourWriter pfnWriter, void *pCBData)
1010 :
1011 : {
1012 : auto cg = new marching_squares::ContourGeneratorOpaque(
1013 : nWidth, nHeight, bNoDataSet, dfNoDataValue, dfContourInterval,
1014 1 : dfContourBase, pfnWriter, pCBData);
1015 :
1016 1 : return reinterpret_cast<GDALContourGeneratorH>(cg);
1017 : }
1018 :
1019 : /************************************************************************/
1020 : /* GDAL_CG_FeedLine() */
1021 : /************************************************************************/
1022 :
1023 : /** Feed a line to the contour generator */
1024 1 : CPLErr GDAL_CG_FeedLine(GDALContourGeneratorH hCG, double *padfScanline)
1025 :
1026 : {
1027 1 : VALIDATE_POINTER1(hCG, "GDAL_CG_FeedLine", CE_Failure);
1028 : return reinterpret_cast<marching_squares::ContourGeneratorOpaque *>(hCG)
1029 1 : ->contourGenerator.feedLine(padfScanline);
1030 : }
1031 :
1032 : /************************************************************************/
1033 : /* GDAL_CG_Destroy() */
1034 : /************************************************************************/
1035 :
1036 : /** Destroy contour generator */
1037 1 : void GDAL_CG_Destroy(GDALContourGeneratorH hCG)
1038 :
1039 : {
1040 1 : delete reinterpret_cast<marching_squares::ContourGeneratorOpaque *>(hCG);
1041 1 : }
|
