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