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