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