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