Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: WAsP Translator
4 : * Purpose: Implements OGRWAsPLayer class.
5 : * Author: Vincent Mora, vincent dot mora at oslandia dot com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2014, Oslandia <info at oslandia dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "ogrwasp.h"
14 : #include "cpl_conv.h"
15 : #include "cpl_string.h"
16 : #include "ogrsf_frmts.h"
17 :
18 : #include <cassert>
19 : #include <map>
20 : #include <sstream>
21 :
22 : /************************************************************************/
23 : /* OGRWAsPLayer() */
24 : /************************************************************************/
25 :
26 9 : OGRWAsPLayer::OGRWAsPLayer(GDALDataset *poDS, const char *pszName,
27 : VSILFILE *hFileHandle,
28 9 : OGRSpatialReference *poSpatialRef)
29 : : m_poDS(poDS), bMerge(false), iFeatureCount(0), sName(pszName),
30 : hFile(hFileHandle), iFirstFieldIdx(0), iSecondFieldIdx(1),
31 9 : iGeomFieldIdx(0), poLayerDefn(new OGRFeatureDefn(pszName)),
32 18 : poSpatialReference(poSpatialRef), iOffsetFeatureBegin(VSIFTellL(hFile)),
33 18 : eMode(READ_ONLY)
34 : {
35 9 : SetDescription(poLayerDefn->GetName());
36 9 : poLayerDefn->Reference();
37 9 : poLayerDefn->GetGeomFieldDefn(0)->SetType(wkbLineString25D);
38 9 : poLayerDefn->GetGeomFieldDefn(0)->SetSpatialRef(poSpatialReference);
39 9 : if (poSpatialReference)
40 1 : poSpatialReference->Reference();
41 9 : }
42 :
43 24 : OGRWAsPLayer::OGRWAsPLayer(
44 : GDALDataset *poDS, const char *pszName, VSILFILE *hFileHandle,
45 : OGRSpatialReference *poSpatialRef, const CPLString &sFirstFieldParam,
46 : const CPLString &sSecondFieldParam, const CPLString &sGeomFieldParam,
47 : bool bMergeParam, double *pdfToleranceParam,
48 24 : double *pdfAdjacentPointToleranceParam, double *pdfPointToCircleRadiusParam)
49 : : m_poDS(poDS), bMerge(bMergeParam), iFeatureCount(0), sName(pszName),
50 : hFile(hFileHandle), sFirstField(sFirstFieldParam),
51 : sSecondField(sSecondFieldParam), sGeomField(sGeomFieldParam),
52 : iFirstFieldIdx(-1), iSecondFieldIdx(-1),
53 24 : iGeomFieldIdx(sGeomFieldParam.empty() ? 0 : -1),
54 24 : poLayerDefn(new OGRFeatureDefn(pszName)),
55 : poSpatialReference(poSpatialRef),
56 48 : iOffsetFeatureBegin(VSIFTellL(hFile)), // Avoids coverity warning.
57 : eMode(WRITE_ONLY), pdfTolerance(pdfToleranceParam),
58 : pdfAdjacentPointTolerance(pdfAdjacentPointToleranceParam),
59 72 : pdfPointToCircleRadius(pdfPointToCircleRadiusParam)
60 : {
61 24 : SetDescription(poLayerDefn->GetName());
62 24 : poLayerDefn->Reference();
63 24 : poLayerDefn->GetGeomFieldDefn(0)->SetType(wkbLineString25D);
64 24 : poLayerDefn->GetGeomFieldDefn(0)->SetSpatialRef(poSpatialReference);
65 24 : if (poSpatialReference)
66 3 : poSpatialReference->Reference();
67 24 : }
68 :
69 : /************************************************************************/
70 : /* ~OGRWAsPLayer() */
71 : /************************************************************************/
72 :
73 66 : OGRWAsPLayer::~OGRWAsPLayer()
74 :
75 : {
76 33 : if (bMerge)
77 : {
78 : /* If polygon where used, we have to merge lines before output */
79 : /* lines must be merged if they have the same left/right values */
80 : /* and touch at end points. */
81 : /* Those lines appear when polygon with the same roughness touch */
82 : /* since the boundary between them is not wanted */
83 : /* We do it here since we are sure we have all polygons */
84 : /* We first detect touching lines, then the kind of touching, */
85 : /* candidates for merging are pairs of neighbors with corresponding */
86 : /* left/right values. Finally we merge */
87 :
88 : typedef std::map<std::pair<double, double>, std::vector<int>> PointMap;
89 48 : PointMap oMap;
90 42 : for (int i = 0; i < static_cast<int>(oBoundaries.size()); i++)
91 : {
92 18 : const Boundary &p = oBoundaries[i];
93 36 : OGRPoint startP, endP;
94 18 : p.poLine->StartPoint(&startP);
95 18 : p.poLine->EndPoint(&endP);
96 18 : oMap[std::make_pair(startP.getX(), startP.getY())].push_back(i);
97 18 : oMap[std::make_pair(endP.getX(), endP.getY())].push_back(i);
98 : }
99 :
100 48 : std::vector<int> endNeighbors(oBoundaries.size(), -1);
101 48 : std::vector<int> startNeighbors(oBoundaries.size(), -1);
102 45 : for (PointMap::const_iterator it = oMap.begin(); it != oMap.end(); ++it)
103 : {
104 21 : if (it->second.size() != 2)
105 21 : continue;
106 0 : int i = it->second[0];
107 0 : int j = it->second[1];
108 :
109 0 : const Boundary &p = oBoundaries[i];
110 0 : OGRPoint startP, endP;
111 0 : p.poLine->StartPoint(&startP);
112 0 : p.poLine->EndPoint(&endP);
113 0 : const Boundary &q = oBoundaries[j];
114 0 : OGRPoint startQ, endQ;
115 0 : q.poLine->StartPoint(&startQ);
116 0 : q.poLine->EndPoint(&endQ);
117 0 : if (isEqual(p.dfRight, q.dfRight) && isEqual(p.dfLeft, q.dfLeft))
118 : {
119 0 : if (endP.Equals(&startQ))
120 : {
121 0 : endNeighbors[i] = j;
122 0 : startNeighbors[j] = i;
123 : }
124 0 : if (endQ.Equals(&startP))
125 : {
126 0 : endNeighbors[j] = i;
127 0 : startNeighbors[i] = j;
128 : }
129 : }
130 0 : if (isEqual(p.dfRight, q.dfLeft) && isEqual(p.dfLeft, q.dfRight))
131 : {
132 0 : if (startP.Equals(&startQ))
133 : {
134 0 : startNeighbors[i] = j;
135 0 : startNeighbors[j] = i;
136 : }
137 0 : if (endP.Equals(&endQ))
138 : {
139 0 : endNeighbors[j] = i;
140 0 : endNeighbors[i] = j;
141 : }
142 : }
143 : }
144 :
145 : /* output all end lines (one neighbor only) and all their neighbors*/
146 24 : if (!oBoundaries.empty())
147 : {
148 6 : std::vector<bool> oHasBeenMerged(oBoundaries.size(), false);
149 21 : for (size_t i = 0; i < oBoundaries.size(); i++)
150 : {
151 36 : if (!oHasBeenMerged[i] &&
152 18 : (startNeighbors[i] < 0 || endNeighbors[i] < 0))
153 : {
154 18 : oHasBeenMerged[i] = true;
155 18 : Boundary *p = &oBoundaries[i];
156 18 : int j = startNeighbors[i] < 0 ? endNeighbors[i]
157 0 : : startNeighbors[i];
158 18 : if (startNeighbors[i] >= 0)
159 : {
160 : /* reverse the line and left/right */
161 0 : p->poLine->reversePoints();
162 0 : std::swap(p->dfLeft, p->dfRight);
163 : }
164 18 : while (j >= 0)
165 : {
166 0 : assert(!oHasBeenMerged[j]);
167 0 : oHasBeenMerged[j] = true;
168 :
169 0 : OGRLineString *other = oBoundaries[j].poLine;
170 0 : OGRPoint endP, startOther;
171 0 : p->poLine->EndPoint(&endP);
172 0 : other->StartPoint(&startOther);
173 0 : if (!endP.Equals(&startOther))
174 0 : other->reversePoints();
175 0 : p->poLine->addSubLineString(other, 1);
176 :
177 : /* next neighbor */
178 0 : if (endNeighbors[j] >= 0 &&
179 0 : !oHasBeenMerged[endNeighbors[j]])
180 0 : j = endNeighbors[j];
181 0 : else if (startNeighbors[j] >= 0 &&
182 0 : !oHasBeenMerged[startNeighbors[j]])
183 0 : j = startNeighbors[j];
184 : else
185 0 : j = -1;
186 : }
187 18 : WriteRoughness(p->poLine, p->dfLeft, p->dfRight);
188 : }
189 : }
190 : /* output all rings */
191 21 : for (size_t i = 0; i < oBoundaries.size(); i++)
192 : {
193 18 : if (oHasBeenMerged[i])
194 18 : continue;
195 0 : oHasBeenMerged[i] = true;
196 0 : Boundary *p = &oBoundaries[i];
197 : int j =
198 0 : startNeighbors[i] < 0 ? endNeighbors[i] : startNeighbors[i];
199 0 : assert(j != -1);
200 0 : if (startNeighbors[i] >= 0)
201 : {
202 : /* reverse the line and left/right */
203 0 : p->poLine->reversePoints();
204 0 : std::swap(p->dfLeft, p->dfRight);
205 : }
206 0 : while (!oHasBeenMerged[j])
207 : {
208 0 : oHasBeenMerged[j] = true;
209 :
210 0 : OGRLineString *other = oBoundaries[j].poLine;
211 0 : OGRPoint endP, startOther;
212 0 : p->poLine->EndPoint(&endP);
213 0 : other->StartPoint(&startOther);
214 0 : if (!endP.Equals(&startOther))
215 0 : other->reversePoints();
216 0 : p->poLine->addSubLineString(other, 1);
217 :
218 : /* next neighbor */
219 0 : if (endNeighbors[j] >= 0)
220 0 : j = endNeighbors[j];
221 0 : else if (startNeighbors[j] >= 0)
222 0 : j = startNeighbors[j];
223 : else
224 0 : assert(false); /* there must be a neighbor since it is a
225 : ring */
226 : }
227 0 : WriteRoughness(p->poLine, p->dfLeft, p->dfRight);
228 : }
229 : }
230 : }
231 : else
232 : {
233 9 : for (size_t i = 0; i < oBoundaries.size(); i++)
234 : {
235 0 : Boundary *p = &oBoundaries[i];
236 0 : WriteRoughness(p->poLine, p->dfLeft, p->dfRight);
237 : }
238 : }
239 33 : poLayerDefn->Release();
240 33 : if (poSpatialReference)
241 4 : poSpatialReference->Release();
242 55 : for (size_t i = 0; i < oZones.size(); i++)
243 22 : delete oZones[i].poPolygon;
244 51 : for (size_t i = 0; i < oBoundaries.size(); i++)
245 18 : delete oBoundaries[i].poLine;
246 66 : }
247 :
248 : /************************************************************************/
249 : /* Simplify() */
250 : /************************************************************************/
251 70 : OGRLineString *OGRWAsPLayer::Simplify(const OGRLineString &line) const
252 : {
253 70 : if (!line.getNumPoints())
254 0 : return line.clone();
255 :
256 : std::unique_ptr<OGRLineString> poLine(
257 80 : (pdfTolerance.get() && *pdfTolerance > 0
258 10 : ? line.Simplify(*pdfTolerance)->toLineString()
259 150 : : line.clone()));
260 :
261 140 : OGRPoint startPt, endPt;
262 70 : poLine->StartPoint(&startPt);
263 70 : poLine->EndPoint(&endPt);
264 70 : const bool isRing = CPL_TO_BOOL(startPt.Equals(&endPt));
265 :
266 70 : if (pdfAdjacentPointTolerance.get() && *pdfAdjacentPointTolerance > 0)
267 : {
268 : /* remove consecutive points that are too close */
269 0 : auto newLine = std::make_unique<OGRLineString>();
270 0 : const double dist = *pdfAdjacentPointTolerance;
271 0 : OGRPoint pt;
272 0 : poLine->StartPoint(&pt);
273 0 : newLine->addPoint(&pt);
274 0 : const int iNumPoints = poLine->getNumPoints();
275 0 : for (int v = 1; v < iNumPoints; v++)
276 : {
277 0 : if (fabs(poLine->getX(v) - pt.getX()) > dist ||
278 0 : fabs(poLine->getY(v) - pt.getY()) > dist)
279 : {
280 0 : poLine->getPoint(v, &pt);
281 0 : newLine->addPoint(&pt);
282 : }
283 : }
284 :
285 : /* force closed loop if initially closed */
286 0 : if (isRing)
287 0 : newLine->setPoint(newLine->getNumPoints() - 1, &startPt);
288 :
289 0 : poLine.reset(newLine.release());
290 : }
291 :
292 70 : if (pdfPointToCircleRadius.get() && *pdfPointToCircleRadius > 0)
293 : {
294 0 : const double radius = *pdfPointToCircleRadius;
295 :
296 : #undef WASP_EXPERIMENTAL_CODE
297 : #ifdef WASP_EXPERIMENTAL_CODE
298 : if (3 == poLine->getNumPoints() && isRing)
299 : {
300 : OGRPoint p0, p1;
301 : poLine->getPoint(0, &p0);
302 : poLine->getPoint(1, &p1);
303 : const double dir[2] = {p1.getX() - p0.getX(),
304 : p1.getY() - p0.getY()};
305 : const double dirNrm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
306 : if (dirNrm > radius)
307 : {
308 : /* convert to rectangle by finding the direction */
309 : /* and offsetting */
310 : const double ortho[2] = {-radius * dir[1] / dirNrm,
311 : radius * dir[0] / dirNrm};
312 : poLine->setNumPoints(5);
313 : poLine->setPoint(0, p0.getX() - ortho[0], p0.getY() - ortho[1]);
314 : poLine->setPoint(1, p1.getX() - ortho[0], p1.getY() - ortho[1]);
315 : poLine->setPoint(2, p1.getX() + ortho[0], p1.getY() + ortho[1]);
316 : poLine->setPoint(3, p0.getX() + ortho[0], p0.getY() + ortho[1]);
317 : poLine->setPoint(4, p0.getX() - ortho[0], p0.getY() - ortho[1]);
318 : }
319 : else
320 : {
321 : /* reduce to a point to be dealt with just after*/
322 : poLine->setNumPoints(1);
323 : poLine->setPoint(0, 0.5 * (p0.getX() + p1.getX()),
324 : 0.5 * (p0.getY() + p1.getY()));
325 : }
326 : }
327 : #endif
328 :
329 0 : if (1 == poLine->getNumPoints())
330 : {
331 0 : const int nbPt = 8;
332 0 : const double cx = poLine->getX(0);
333 0 : const double cy = poLine->getY(0);
334 0 : poLine->setNumPoints(nbPt + 1);
335 0 : for (int v = 0; v <= nbPt; v++)
336 : {
337 : /* the % is necessary to make sure the ring */
338 : /* is really closed and not open due to */
339 : /* roundoff error of cos(2pi) and sin(2pi) */
340 0 : poLine->setPoint(
341 0 : v, cx + radius * cos((v % nbPt) * (2 * M_PI / nbPt)),
342 0 : cy + radius * sin((v % nbPt) * (2 * M_PI / nbPt)));
343 : }
344 : }
345 : }
346 :
347 70 : return poLine.release();
348 : }
349 :
350 : /************************************************************************/
351 : /* WriteElevation() */
352 : /************************************************************************/
353 :
354 42 : OGRErr OGRWAsPLayer::WriteElevation(OGRLineString *poGeom, const double &dfZ)
355 :
356 : {
357 84 : std::unique_ptr<OGRLineString> poLine(Simplify(*poGeom));
358 :
359 42 : const int iNumPoints = poLine->getNumPoints();
360 42 : if (!iNumPoints)
361 0 : return OGRERR_NONE; /* empty geom */
362 :
363 42 : VSIFPrintfL(hFile, "%11.3f %11d", dfZ, iNumPoints);
364 :
365 156 : for (int v = 0; v < iNumPoints; v++)
366 : {
367 114 : if (!(v % 3))
368 42 : VSIFPrintfL(hFile, "\n");
369 114 : VSIFPrintfL(hFile, "%11.1f %11.1f ", poLine->getX(v), poLine->getY(v));
370 : }
371 42 : VSIFPrintfL(hFile, "\n");
372 :
373 42 : return OGRERR_NONE;
374 : }
375 :
376 45 : OGRErr OGRWAsPLayer::WriteElevation(OGRGeometry *poGeom, const double &dfZ)
377 :
378 : {
379 45 : switch (poGeom->getGeometryType())
380 : {
381 42 : case wkbLineString:
382 : case wkbLineString25D:
383 42 : return WriteElevation(poGeom->toLineString(), dfZ);
384 1 : case wkbMultiLineString25D:
385 : case wkbMultiLineString:
386 : {
387 2 : for (auto &&poMember : poGeom->toGeometryCollection())
388 : {
389 1 : const OGRErr err = WriteElevation(poMember, dfZ);
390 1 : if (OGRERR_NONE != err)
391 0 : return err;
392 : }
393 1 : return OGRERR_NONE;
394 : }
395 2 : default:
396 : {
397 2 : CPLError(CE_Failure, CPLE_NotSupported,
398 : "Cannot handle geometry of type %s",
399 2 : OGRGeometryTypeToName(poGeom->getGeometryType()));
400 2 : break;
401 : }
402 : }
403 2 : return OGRERR_FAILURE; /* avoid visual warning */
404 : }
405 :
406 : /************************************************************************/
407 : /* WriteRoughness() */
408 : /************************************************************************/
409 :
410 22 : OGRErr OGRWAsPLayer::WriteRoughness(OGRPolygon *poGeom, const double &dfZ)
411 :
412 : {
413 : /* text intersection with polygons in the stack */
414 : /* for linestrings intersections, write linestring */
415 : /* for polygon intersection error */
416 : /* for point intersection do nothing */
417 :
418 22 : OGRErr err = OGRERR_NONE;
419 22 : OGREnvelope oEnvelope;
420 22 : poGeom->getEnvelope(&oEnvelope);
421 67 : for (size_t i = 0; i < oZones.size(); i++)
422 : {
423 : const bool bIntersects =
424 45 : CPL_TO_BOOL(oEnvelope.Intersects(oZones[i].oEnvelope));
425 90 : if (bIntersects &&
426 45 : (!bMerge || !isEqual(dfZ, oZones[i].dfZ))) /* boundary */
427 : {
428 : OGRGeometry *poIntersection =
429 39 : oZones[i].poPolygon->Intersection(poGeom);
430 39 : if (poIntersection)
431 : {
432 39 : switch (poIntersection->getGeometryType())
433 : {
434 18 : case wkbLineString:
435 : case wkbLineString25D:
436 : {
437 18 : Boundary oB = {poIntersection->toLineString()->clone(),
438 18 : dfZ, oZones[i].dfZ};
439 18 : oBoundaries.push_back(oB);
440 : }
441 18 : break;
442 0 : case wkbMultiLineString:
443 : case wkbMultiLineString25D:
444 : {
445 : /*TODO join the multilinestring into linestring*/
446 0 : OGRLineString *poLine = nullptr;
447 0 : OGRPoint *poStart = new OGRPoint;
448 0 : OGRPoint *poEnd = new OGRPoint;
449 0 : for (auto &&poSubLine :
450 0 : poIntersection->toMultiLineString())
451 : {
452 0 : poSubLine->StartPoint(poStart);
453 :
454 0 : if (poLine == nullptr)
455 : {
456 0 : poLine = poSubLine->clone();
457 : }
458 0 : else if (poLine->getNumPoints() == 0 ||
459 0 : poStart->Equals(poEnd))
460 : {
461 0 : poLine->addSubLineString(poSubLine, 1);
462 : }
463 : else
464 : {
465 0 : Boundary oB = {poLine, dfZ, oZones[i].dfZ};
466 0 : oBoundaries.push_back(oB);
467 0 : poLine = poSubLine->clone();
468 : }
469 0 : poLine->EndPoint(poEnd);
470 : }
471 0 : Boundary oB = {poLine, dfZ, oZones[i].dfZ};
472 0 : oBoundaries.push_back(oB);
473 0 : delete poStart;
474 0 : delete poEnd;
475 : }
476 0 : break;
477 0 : case wkbPolygon:
478 : case wkbPolygon25D:
479 : {
480 0 : OGREnvelope oErrorRegion = oZones[i].oEnvelope;
481 0 : oErrorRegion.Intersect(oEnvelope);
482 0 : CPLError(CE_Failure, CPLE_NotSupported,
483 : "Overlapping polygons in rectangle (%.16g "
484 : "%.16g, %.16g %.16g))",
485 : oErrorRegion.MinX, oErrorRegion.MinY,
486 : oErrorRegion.MaxX, oErrorRegion.MaxY);
487 0 : err = OGRERR_FAILURE;
488 : }
489 0 : break;
490 0 : case wkbGeometryCollection:
491 : case wkbGeometryCollection25D:
492 : {
493 0 : for (auto &&poMember :
494 0 : poIntersection->toGeometryCollection())
495 : {
496 : const OGRwkbGeometryType eType =
497 0 : poMember->getGeometryType();
498 0 : if (wkbFlatten(eType) == wkbPolygon)
499 : {
500 0 : OGREnvelope oErrorRegion = oZones[i].oEnvelope;
501 0 : oErrorRegion.Intersect(oEnvelope);
502 0 : CPLError(CE_Failure, CPLE_NotSupported,
503 : "Overlapping polygons in rectangle "
504 : "(%.16g %.16g, %.16g %.16g))",
505 : oErrorRegion.MinX, oErrorRegion.MinY,
506 : oErrorRegion.MaxX, oErrorRegion.MaxY);
507 0 : err = OGRERR_FAILURE;
508 : }
509 : }
510 : }
511 0 : break;
512 21 : case wkbPoint:
513 : case wkbPoint25D:
514 : /* do nothing */
515 21 : break;
516 0 : default:
517 0 : CPLError(CE_Failure, CPLE_NotSupported,
518 : "Unhandled polygon intersection of type %s",
519 : OGRGeometryTypeToName(
520 0 : poIntersection->getGeometryType()));
521 0 : err = OGRERR_FAILURE;
522 : }
523 : }
524 39 : delete poIntersection;
525 : }
526 : }
527 :
528 22 : Zone oZ = {oEnvelope, poGeom->clone(), dfZ};
529 22 : oZones.push_back(oZ);
530 22 : return err;
531 : }
532 :
533 28 : OGRErr OGRWAsPLayer::WriteRoughness(OGRLineString *poGeom,
534 : const double &dfZleft,
535 : const double &dfZright)
536 :
537 : {
538 56 : std::unique_ptr<OGRLineString> poLine(Simplify(*poGeom));
539 :
540 28 : const int iNumPoints = poLine->getNumPoints();
541 28 : if (!iNumPoints)
542 0 : return OGRERR_NONE; /* empty geom */
543 :
544 28 : VSIFPrintfL(hFile, "%11.3f %11.3f %11d", dfZleft, dfZright, iNumPoints);
545 :
546 94 : for (int v = 0; v < iNumPoints; v++)
547 : {
548 66 : if (!(v % 3))
549 28 : VSIFPrintfL(hFile, "\n ");
550 66 : VSIFPrintfL(hFile, "%11.1f %11.1f ", poLine->getX(v), poLine->getY(v));
551 : }
552 28 : VSIFPrintfL(hFile, "\n");
553 :
554 28 : return OGRERR_NONE;
555 : }
556 :
557 34 : OGRErr OGRWAsPLayer::WriteRoughness(OGRGeometry *poGeom, const double &dfZleft,
558 : const double &dfZright)
559 :
560 : {
561 34 : switch (poGeom->getGeometryType())
562 : {
563 10 : case wkbLineString:
564 : case wkbLineString25D:
565 10 : return WriteRoughness(poGeom->toLineString(), dfZleft, dfZright);
566 22 : case wkbPolygon:
567 : case wkbPolygon25D:
568 22 : return WriteRoughness(poGeom->toPolygon(), dfZleft);
569 2 : case wkbMultiPolygon:
570 : case wkbMultiPolygon25D:
571 : case wkbMultiLineString25D:
572 : case wkbMultiLineString:
573 : {
574 4 : for (auto &&poMember : poGeom->toGeometryCollection())
575 : {
576 2 : const OGRErr err = WriteRoughness(poMember, dfZleft, dfZright);
577 2 : if (OGRERR_NONE != err)
578 0 : return err;
579 : }
580 2 : return OGRERR_NONE;
581 : }
582 0 : default:
583 : {
584 0 : CPLError(CE_Failure, CPLE_NotSupported,
585 : "Cannot handle geometry of type %s",
586 0 : OGRGeometryTypeToName(poGeom->getGeometryType()));
587 0 : break;
588 : }
589 : }
590 0 : return OGRERR_FAILURE; /* avoid visual warning */
591 : }
592 :
593 : /************************************************************************/
594 : /* ICreateFeature() */
595 : /************************************************************************/
596 :
597 100 : OGRErr OGRWAsPLayer::ICreateFeature(OGRFeature *poFeature)
598 :
599 : {
600 100 : if (WRITE_ONLY != eMode)
601 : {
602 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Layer is open read only");
603 0 : return OGRERR_FAILURE;
604 : }
605 :
606 : /* This mainly checks for errors or inconsistencies */
607 : /* the real work is done by WriteElevation or WriteRoughness */
608 100 : if (-1 == iFirstFieldIdx && !sFirstField.empty())
609 : {
610 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Cannot find field %s",
611 : sFirstField.c_str());
612 0 : return OGRERR_FAILURE;
613 : }
614 100 : if (-1 == iSecondFieldIdx && !sSecondField.empty())
615 : {
616 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Cannot find field %s",
617 : sSecondField.c_str());
618 0 : return OGRERR_FAILURE;
619 : }
620 100 : if (-1 == iGeomFieldIdx && !sGeomField.empty())
621 : {
622 0 : CPLError(CE_Failure, CPLE_IllegalArg, "Cannot find field %s",
623 : sSecondField.c_str());
624 0 : return OGRERR_FAILURE;
625 : }
626 100 : OGRGeometry *geom = poFeature->GetGeomFieldRef(iGeomFieldIdx);
627 100 : if (!geom)
628 16 : return OGRERR_NONE; /* null geom, nothing to do */
629 :
630 84 : const OGRwkbGeometryType geomType = geom->getGeometryType();
631 84 : const bool bPolygon =
632 82 : (geomType == wkbPolygon) || (geomType == wkbPolygon25D) ||
633 166 : (geomType == wkbMultiPolygon) || (geomType == wkbMultiPolygon25D);
634 84 : const bool bRoughness = (-1 != iSecondFieldIdx) || bPolygon;
635 :
636 84 : double z1 = 0.0;
637 84 : if (-1 != iFirstFieldIdx)
638 : {
639 26 : if (!poFeature->IsFieldSetAndNotNull(iFirstFieldIdx))
640 : {
641 0 : CPLError(CE_Failure, CPLE_NotSupported, "Field %d %s is NULL",
642 : iFirstFieldIdx, sFirstField.c_str());
643 0 : return OGRERR_FAILURE;
644 : }
645 26 : z1 = poFeature->GetFieldAsDouble(iFirstFieldIdx);
646 : }
647 : else
648 : {
649 : /* Case of z value for elevation or roughness, so we compute it */
650 58 : OGRPoint centroid;
651 58 : if (geom->getCoordinateDimension() != 3)
652 : {
653 :
654 8 : CPLError(CE_Failure, CPLE_NotSupported,
655 : "No field defined and no Z coordinate");
656 8 : return OGRERR_FAILURE;
657 : }
658 50 : z1 = AvgZ(geom);
659 : }
660 :
661 76 : double z2 = 0.0;
662 76 : if (-1 != iSecondFieldIdx)
663 : {
664 10 : if (!poFeature->IsFieldSetAndNotNull(iSecondFieldIdx))
665 : {
666 0 : CPLError(CE_Failure, CPLE_NotSupported, "Field %d %s is NULL",
667 : iSecondFieldIdx, sSecondField.c_str());
668 0 : return OGRERR_FAILURE;
669 : }
670 10 : z2 = poFeature->GetFieldAsDouble(iSecondFieldIdx);
671 : }
672 66 : else if (bRoughness && !bPolygon)
673 : {
674 0 : CPLError(CE_Failure, CPLE_IllegalArg, "No right roughness field");
675 0 : return OGRERR_FAILURE;
676 : }
677 :
678 76 : return bRoughness ? WriteRoughness(geom, z1, z2) : WriteElevation(geom, z1);
679 : }
680 :
681 : /************************************************************************/
682 : /* CreateField() */
683 : /************************************************************************/
684 :
685 49 : OGRErr OGRWAsPLayer::CreateField(const OGRFieldDefn *poField,
686 : CPL_UNUSED int bApproxOK)
687 : {
688 49 : poLayerDefn->AddFieldDefn(poField);
689 :
690 : /* Update field indexes */
691 49 : if (-1 == iFirstFieldIdx && !sFirstField.empty())
692 4 : iFirstFieldIdx = poLayerDefn->GetFieldIndex(sFirstField.c_str());
693 49 : if (-1 == iSecondFieldIdx && !sSecondField.empty())
694 3 : iSecondFieldIdx = poLayerDefn->GetFieldIndex(sSecondField.c_str());
695 :
696 49 : return OGRERR_NONE;
697 : }
698 :
699 : /************************************************************************/
700 : /* CreateGeomField() */
701 : /************************************************************************/
702 :
703 0 : OGRErr OGRWAsPLayer::CreateGeomField(const OGRGeomFieldDefn *poGeomFieldIn,
704 : CPL_UNUSED int bApproxOK)
705 : {
706 0 : OGRGeomFieldDefn oFieldDefn(poGeomFieldIn);
707 0 : auto poSRSOri = poGeomFieldIn->GetSpatialRef();
708 0 : if (poSRSOri)
709 : {
710 0 : auto poSRS = poSRSOri->Clone();
711 0 : poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
712 0 : oFieldDefn.SetSpatialRef(poSRS);
713 0 : poSRS->Release();
714 : }
715 0 : poLayerDefn->AddGeomFieldDefn(&oFieldDefn);
716 :
717 : /* Update geom field index */
718 0 : if (-1 == iGeomFieldIdx)
719 : {
720 0 : iGeomFieldIdx = poLayerDefn->GetGeomFieldIndex(sGeomField.c_str());
721 : }
722 :
723 0 : return OGRERR_NONE;
724 : }
725 :
726 : /************************************************************************/
727 : /* GetNextRawFeature() */
728 : /************************************************************************/
729 :
730 19 : OGRFeature *OGRWAsPLayer::GetNextRawFeature()
731 :
732 : {
733 19 : if (READ_ONLY != eMode)
734 : {
735 8 : CPLError(CE_Failure, CPLE_IllegalArg, "Layer is open write only");
736 8 : return nullptr;
737 : }
738 :
739 11 : const char *pszLine = CPLReadLineL(hFile);
740 11 : if (!pszLine)
741 1 : return nullptr;
742 :
743 10 : double dfValues[4] = {0};
744 10 : int iNumValues = 0;
745 : {
746 20 : std::istringstream iss(pszLine);
747 30 : while (iNumValues < 4 && (iss >> dfValues[iNumValues]))
748 : {
749 20 : ++iNumValues;
750 : }
751 :
752 10 : if (iNumValues < 2)
753 : {
754 0 : if (iNumValues)
755 0 : CPLError(CE_Failure, CPLE_FileIO, "No enough values");
756 0 : return nullptr;
757 : }
758 : }
759 :
760 10 : if (poLayerDefn->GetFieldCount() != iNumValues - 1)
761 : {
762 0 : CPLError(CE_Failure, CPLE_FileIO,
763 : "looking for %d values and found %d on line: %s",
764 0 : poLayerDefn->GetFieldCount(), iNumValues - 1, pszLine);
765 0 : return nullptr;
766 : }
767 10 : const double dfNumPairToRead = dfValues[iNumValues - 1];
768 10 : if (!(dfNumPairToRead >= 0 && dfNumPairToRead < 1000000) ||
769 10 : static_cast<int>(dfNumPairToRead) != dfNumPairToRead)
770 : {
771 0 : CPLError(CE_Failure, CPLE_FileIO, "Invalid coordinate number: %f",
772 : dfNumPairToRead);
773 0 : return nullptr;
774 : }
775 :
776 20 : auto poFeature = std::make_unique<OGRFeature>(poLayerDefn);
777 10 : poFeature->SetFID(++iFeatureCount);
778 20 : for (int i = 0; i < iNumValues - 1; i++)
779 10 : poFeature->SetField(i, dfValues[i]);
780 :
781 10 : const int iNumValuesToRead = static_cast<int>(2 * dfNumPairToRead);
782 10 : int iReadValues = 0;
783 20 : std::vector<double> values(iNumValuesToRead);
784 20 : for (pszLine = CPLReadLineL(hFile); pszLine;
785 10 : pszLine = iNumValuesToRead > iReadValues ? CPLReadLineL(hFile)
786 : : nullptr)
787 : {
788 30 : std::istringstream iss(pszLine);
789 70 : while (iNumValuesToRead > iReadValues && (iss >> values[iReadValues]))
790 : {
791 60 : ++iReadValues;
792 : }
793 : }
794 10 : if (iNumValuesToRead != iReadValues)
795 : {
796 0 : CPLError(CE_Failure, CPLE_FileIO, "No enough values for linestring");
797 0 : return nullptr;
798 : }
799 20 : auto poLine = std::make_unique<OGRLineString>();
800 10 : poLine->setCoordinateDimension(3);
801 10 : poLine->assignSpatialReference(poSpatialReference);
802 40 : for (int i = 0; i < iNumValuesToRead; i += 2)
803 : {
804 30 : poLine->addPoint(values[i], values[i + 1], 0);
805 : }
806 10 : poFeature->SetGeomFieldDirectly(0, poLine.release());
807 :
808 10 : return poFeature.release();
809 : }
810 :
811 : /************************************************************************/
812 : /* TestCapability() */
813 : /************************************************************************/
814 :
815 48 : int OGRWAsPLayer::TestCapability(const char *pszCap)
816 :
817 : {
818 88 : return (WRITE_ONLY == eMode && (EQUAL(pszCap, OLCSequentialWrite) ||
819 40 : EQUAL(pszCap, OLCCreateField) ||
820 32 : EQUAL(pszCap, OLCCreateGeomField) ||
821 80 : EQUAL(pszCap, OLCZGeometries)));
822 : }
823 :
824 : /************************************************************************/
825 : /* TestCapability() */
826 : /************************************************************************/
827 :
828 8 : void OGRWAsPLayer::ResetReading()
829 : {
830 8 : iFeatureCount = 0;
831 8 : VSIFSeekL(hFile, iOffsetFeatureBegin, SEEK_SET);
832 8 : }
833 :
834 : /************************************************************************/
835 : /* AvgZ() */
836 : /************************************************************************/
837 :
838 48 : double OGRWAsPLayer::AvgZ(OGRLineString *poGeom)
839 :
840 : {
841 48 : const int iNumPoints = poGeom->getNumPoints();
842 48 : double sum = 0;
843 210 : for (int v = 0; v < iNumPoints; v++)
844 : {
845 162 : sum += poGeom->getZ(v);
846 : }
847 48 : return iNumPoints ? sum / iNumPoints : 0;
848 : }
849 :
850 16 : double OGRWAsPLayer::AvgZ(OGRPolygon *poGeom)
851 :
852 : {
853 16 : return AvgZ(poGeom->getExteriorRing());
854 : }
855 :
856 3 : double OGRWAsPLayer::AvgZ(OGRGeometryCollection *poGeom)
857 :
858 : {
859 3 : return poGeom->getNumGeometries() ? AvgZ(poGeom->getGeometryRef(0)) : 0;
860 : }
861 :
862 53 : double OGRWAsPLayer::AvgZ(OGRGeometry *poGeom)
863 :
864 : {
865 53 : switch (poGeom->getGeometryType())
866 : {
867 32 : case wkbLineString:
868 : case wkbLineString25D:
869 32 : return AvgZ(static_cast<OGRLineString *>(poGeom));
870 16 : case wkbPolygon:
871 : case wkbPolygon25D:
872 16 : return AvgZ(static_cast<OGRPolygon *>(poGeom));
873 3 : case wkbMultiLineString:
874 : case wkbMultiLineString25D:
875 :
876 : case wkbMultiPolygon:
877 : case wkbMultiPolygon25D:
878 3 : return AvgZ(static_cast<OGRGeometryCollection *>(poGeom));
879 2 : default:
880 2 : CPLError(CE_Warning, CPLE_NotSupported,
881 : "Unsupported geometry type in OGRWAsPLayer::AvgZ()");
882 2 : break;
883 : }
884 2 : return 0; /* avoid warning */
885 : }
886 :
887 : /************************************************************************/
888 : /* DouglasPeucker() */
889 : /************************************************************************/
890 :
891 : // void DouglasPeucker(PointList[], epsilon)
892 : //
893 : //{
894 : // // Find the point with the maximum distance
895 : // double dmax = 0;
896 : // int index = 0;
897 : // int end = length(PointList).
898 : // for (int i = 1; i<end; i++)
899 : // {
900 : // const double d = shortestDistanceToSegment(PointList[i],
901 : // Line(PointList[0], PointList[end])) if ( d > dmax )
902 : // {
903 : // index = i
904 : // dmax = d
905 : // }
906 : // }
907 : // // If max distance is greater than epsilon, recursively simplify
908 : // if ( dmax > epsilon )
909 : // {
910 : // // Recursive call
911 : // recResults1[] = DouglasPeucker(PointList[1...index], epsilon)
912 : // recResults2[] = DouglasPeucker(PointList[index...end], epsilon)
913 : //
914 : // // Build the result list
915 : // ResultList[] = {recResults1[1...end-1] recResults2[1...end]}
916 : // }
917 : // else
918 : // {
919 : // ResultList[] = {PointList[1], PointList[end]}
920 : // }
921 : // // Return the result
922 : // return ResultList[]
923 : // }
|
