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