Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Interlis 1 Translator
4 : * Purpose: Implements OGRILI1Layer class.
5 : * Author: Pirmin Kalberer, Sourcepole AG
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, Pirmin Kalberer, Sourcepole AG
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_conv.h"
31 : #include "cpl_string.h"
32 : #include "ogr_geos.h"
33 : #include "ogr_ili1.h"
34 :
35 : #include <map>
36 : #include <vector>
37 :
38 : /************************************************************************/
39 : /* OGRILI1Layer() */
40 : /************************************************************************/
41 :
42 152 : OGRILI1Layer::OGRILI1Layer(OGRFeatureDefn *poFeatureDefnIn,
43 : const GeomFieldInfos &oGeomFieldInfosIn,
44 152 : OGRILI1DataSource *poDSIn)
45 : : poFeatureDefn(poFeatureDefnIn), oGeomFieldInfos(oGeomFieldInfosIn),
46 : nFeatures(0), papoFeatures(nullptr), nFeatureIdx(0), bGeomsJoined(FALSE),
47 152 : poDS(poDSIn)
48 : {
49 152 : SetDescription(poFeatureDefn->GetName());
50 152 : poFeatureDefn->Reference();
51 152 : }
52 :
53 : /************************************************************************/
54 : /* ~OGRILI1Layer() */
55 : /************************************************************************/
56 :
57 304 : OGRILI1Layer::~OGRILI1Layer()
58 : {
59 418 : for (int i = 0; i < nFeatures; i++)
60 : {
61 266 : delete papoFeatures[i];
62 : }
63 152 : CPLFree(papoFeatures);
64 :
65 152 : if (poFeatureDefn)
66 152 : poFeatureDefn->Release();
67 304 : }
68 :
69 266 : OGRErr OGRILI1Layer::AddFeature(OGRFeature *poFeature)
70 : {
71 266 : nFeatures++;
72 :
73 266 : papoFeatures = static_cast<OGRFeature **>(
74 266 : CPLRealloc(papoFeatures, sizeof(void *) * nFeatures));
75 :
76 266 : papoFeatures[nFeatures - 1] = poFeature;
77 :
78 266 : return OGRERR_NONE;
79 : }
80 :
81 : /************************************************************************/
82 : /* ResetReading() */
83 : /************************************************************************/
84 :
85 49 : void OGRILI1Layer::ResetReading()
86 : {
87 49 : nFeatureIdx = 0;
88 49 : }
89 :
90 : /************************************************************************/
91 : /* GetNextFeature() */
92 : /************************************************************************/
93 :
94 51 : OGRFeature *OGRILI1Layer::GetNextFeature()
95 : {
96 51 : if (!bGeomsJoined)
97 36 : JoinGeomLayers();
98 :
99 51 : while (nFeatureIdx < nFeatures)
100 : {
101 34 : OGRFeature *poFeature = GetNextFeatureRef();
102 34 : if (poFeature)
103 34 : return poFeature->Clone();
104 : }
105 17 : return nullptr;
106 : }
107 :
108 109 : OGRFeature *OGRILI1Layer::GetNextFeatureRef()
109 : {
110 109 : OGRFeature *poFeature = nullptr;
111 109 : if (nFeatureIdx < nFeatures)
112 : {
113 102 : poFeature = papoFeatures[nFeatureIdx++];
114 : // apply filters
115 204 : if ((m_poFilterGeom == nullptr ||
116 204 : FilterGeometry(poFeature->GetGeometryRef())) &&
117 102 : (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
118 102 : return poFeature;
119 : }
120 7 : return nullptr;
121 : }
122 :
123 : /************************************************************************/
124 : /* GetFeatureRef() */
125 : /************************************************************************/
126 :
127 0 : OGRFeature *OGRILI1Layer::GetFeatureRef(GIntBig nFID)
128 :
129 : {
130 0 : ResetReading();
131 :
132 0 : OGRFeature *poFeature = nullptr;
133 0 : while ((poFeature = GetNextFeatureRef()) != nullptr)
134 : {
135 0 : if (poFeature->GetFID() == nFID)
136 0 : return poFeature;
137 : }
138 :
139 0 : return nullptr;
140 : }
141 :
142 19 : OGRFeature *OGRILI1Layer::GetFeatureRef(const char *fid)
143 :
144 : {
145 19 : ResetReading();
146 :
147 19 : OGRFeature *poFeature = nullptr;
148 41 : while ((poFeature = GetNextFeatureRef()) != nullptr)
149 : {
150 41 : if (!strcmp(poFeature->GetFieldAsString(0), fid))
151 19 : return poFeature;
152 : }
153 :
154 0 : return nullptr;
155 : }
156 :
157 : /************************************************************************/
158 : /* GetFeatureCount() */
159 : /************************************************************************/
160 :
161 189 : GIntBig OGRILI1Layer::GetFeatureCount(int bForce)
162 : {
163 189 : if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr
164 : /* && poAreaLineLayer == NULL*/)
165 : {
166 189 : return nFeatures;
167 : }
168 : else
169 : {
170 0 : return OGRLayer::GetFeatureCount(bForce);
171 : }
172 : }
173 :
174 : #ifndef d2str_defined
175 : #define d2str_defined
176 :
177 346 : static const char *d2str(double val)
178 : {
179 346 : if (val == (int)val)
180 342 : return CPLSPrintf("%d", (int)val);
181 4 : if (fabs(val) < 370)
182 4 : return CPLSPrintf("%.16g", val);
183 0 : if (fabs(val) > 100000000.0)
184 0 : return CPLSPrintf("%.16g", val);
185 :
186 0 : return CPLSPrintf("%.3f", val);
187 : }
188 : #endif
189 :
190 37 : static void AppendCoordinateList(OGRLineString *poLine, OGRILI1DataSource *poDS)
191 : {
192 37 : const bool b3D = CPL_TO_BOOL(wkbHasZ(poLine->getGeometryType()));
193 :
194 165 : for (int iPoint = 0; iPoint < poLine->getNumPoints(); iPoint++)
195 : {
196 128 : if (iPoint == 0)
197 37 : VSIFPrintfL(poDS->GetTransferFile(), "STPT");
198 : else
199 91 : VSIFPrintfL(poDS->GetTransferFile(), "LIPT");
200 128 : VSIFPrintfL(poDS->GetTransferFile(), " %s",
201 : d2str(poLine->getX(iPoint)));
202 128 : VSIFPrintfL(poDS->GetTransferFile(), " %s",
203 : d2str(poLine->getY(iPoint)));
204 128 : if (b3D)
205 63 : VSIFPrintfL(poDS->GetTransferFile(), " %s",
206 : d2str(poLine->getZ(iPoint)));
207 128 : VSIFPrintfL(poDS->GetTransferFile(), "\n");
208 : }
209 37 : VSIFPrintfL(poDS->GetTransferFile(), "ELIN\n");
210 37 : }
211 :
212 1 : static void AppendCompoundCurve(OGRCompoundCurve *poCC, OGRILI1DataSource *poDS)
213 : {
214 3 : for (int iMember = 0; iMember < poCC->getNumCurves(); iMember++)
215 : {
216 2 : OGRCurve *poGeometry = poCC->getCurve(iMember);
217 2 : int b3D = wkbHasZ(poGeometry->getGeometryType());
218 3 : int bIsArc = (poGeometry->getGeometryType() == wkbCircularString ||
219 1 : poGeometry->getGeometryType() == wkbCircularStringZ);
220 2 : OGRSimpleCurve *poLine = (OGRSimpleCurve *)poGeometry;
221 7 : for (int iPoint = 0; iPoint < poLine->getNumPoints(); iPoint++)
222 : {
223 : // Skip last point in curve member
224 7 : if (iPoint == poLine->getNumPoints() - 1 &&
225 2 : iMember < poCC->getNumCurves() - 1)
226 1 : continue;
227 4 : if (iMember == 0 && iPoint == 0)
228 1 : VSIFPrintfL(poDS->GetTransferFile(), "STPT");
229 3 : else if (bIsArc && iPoint == 1)
230 1 : VSIFPrintfL(poDS->GetTransferFile(), "ARCP");
231 : else
232 2 : VSIFPrintfL(poDS->GetTransferFile(), "LIPT");
233 4 : VSIFPrintfL(poDS->GetTransferFile(), " %s",
234 : d2str(poLine->getX(iPoint)));
235 4 : VSIFPrintfL(poDS->GetTransferFile(), " %s",
236 : d2str(poLine->getY(iPoint)));
237 4 : if (b3D)
238 0 : VSIFPrintfL(poDS->GetTransferFile(), " %s",
239 : d2str(poLine->getZ(iPoint)));
240 4 : VSIFPrintfL(poDS->GetTransferFile(), "\n");
241 : }
242 : }
243 1 : VSIFPrintfL(poDS->GetTransferFile(), "ELIN\n");
244 1 : }
245 :
246 101 : int OGRILI1Layer::GeometryAppend(OGRGeometry *poGeometry)
247 : {
248 : #ifdef DEBUG_VERBOSE
249 : CPLDebug("OGR_ILI", "OGRILI1Layer::GeometryAppend OGRGeometryType: %s",
250 : OGRGeometryTypeToName(poGeometry->getGeometryType()));
251 : #endif
252 : /* -------------------------------------------------------------------- */
253 : /* 2D Point */
254 : /* -------------------------------------------------------------------- */
255 101 : if (poGeometry->getGeometryType() == wkbPoint)
256 : {
257 : /* embedded in from non-geometry fields */
258 : }
259 : /* -------------------------------------------------------------------- */
260 : /* 3D Point */
261 : /* -------------------------------------------------------------------- */
262 90 : else if (poGeometry->getGeometryType() == wkbPoint25D)
263 : {
264 : /* embedded in from non-geometry fields */
265 : }
266 :
267 : /* -------------------------------------------------------------------- */
268 : /* LineString and LinearRing */
269 : /* -------------------------------------------------------------------- */
270 143 : else if (poGeometry->getGeometryType() == wkbLineString ||
271 62 : poGeometry->getGeometryType() == wkbLineString25D)
272 : {
273 37 : AppendCoordinateList(poGeometry->toLineString(), poDS);
274 : }
275 :
276 : /* -------------------------------------------------------------------- */
277 : /* Polygon */
278 : /* -------------------------------------------------------------------- */
279 79 : else if (poGeometry->getGeometryType() == wkbPolygon ||
280 35 : poGeometry->getGeometryType() == wkbPolygon25D)
281 : {
282 18 : OGRPolygon *poPolygon = poGeometry->toPolygon();
283 36 : for (auto &&poRing : poPolygon)
284 : {
285 18 : if (!GeometryAppend(poRing))
286 0 : return FALSE;
287 : }
288 : }
289 :
290 : /* -------------------------------------------------------------------- */
291 : /* MultiPolygon */
292 : /* -------------------------------------------------------------------- */
293 26 : else if (wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon ||
294 20 : wkbFlatten(poGeometry->getGeometryType()) == wkbMultiLineString ||
295 14 : wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPoint ||
296 8 : wkbFlatten(poGeometry->getGeometryType()) ==
297 2 : wkbGeometryCollection ||
298 47 : wkbFlatten(poGeometry->getGeometryType()) == wkbMultiCurve ||
299 1 : wkbFlatten(poGeometry->getGeometryType()) == wkbMultiCurveZ)
300 : {
301 25 : OGRGeometryCollection *poGC = poGeometry->toGeometryCollection();
302 :
303 : #if 0
304 : // TODO: Why this large NOP block?
305 : if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon )
306 : {
307 : }
308 : else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiLineString )
309 : {
310 : }
311 : else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPoint )
312 : {
313 : }
314 : else
315 : {
316 : }
317 : #endif
318 :
319 62 : for (auto &&poMember : poGC)
320 : {
321 37 : if (!GeometryAppend(poMember))
322 0 : return FALSE;
323 : }
324 : }
325 1 : else if (poGeometry->getGeometryType() == wkbCompoundCurve ||
326 0 : poGeometry->getGeometryType() == wkbCompoundCurveZ)
327 : {
328 1 : AppendCompoundCurve(poGeometry->toCompoundCurve(), poDS);
329 : }
330 : else
331 : {
332 0 : CPLError(CE_Warning, CPLE_AppDefined,
333 : "Skipping unknown geometry type '%s'",
334 0 : OGRGeometryTypeToName(poGeometry->getGeometryType()));
335 0 : return FALSE;
336 : }
337 :
338 101 : return TRUE;
339 : }
340 :
341 : /************************************************************************/
342 : /* ICreateFeature() */
343 : /************************************************************************/
344 :
345 84 : OGRErr OGRILI1Layer::ICreateFeature(OGRFeature *poFeature)
346 : {
347 : // FIXME: tid not thread safe
348 : static long tid = -1; // system generated TID (must be unique within table)
349 84 : VSILFILE *fp = poDS->GetTransferFile();
350 84 : VSIFPrintfL(fp, "OBJE");
351 :
352 168 : if (poFeatureDefn->GetFieldCount() &&
353 84 : !EQUAL(poFeatureDefn->GetFieldDefn(0)->GetNameRef(), "TID"))
354 : {
355 : // Input is not generated from an Interlis 1 source
356 84 : if (poFeature->GetFID() != OGRNullFID)
357 0 : tid = (int)poFeature->GetFID();
358 : else
359 84 : ++tid;
360 84 : VSIFPrintfL(fp, " %ld", tid);
361 : // Embedded geometry
362 84 : if (poFeature->GetGeometryRef() != nullptr)
363 : {
364 46 : OGRGeometry *poGeometry = poFeature->GetGeometryRef();
365 : // 2D Point
366 46 : if (poGeometry->getGeometryType() == wkbPoint)
367 : {
368 5 : OGRPoint *poPoint = poGeometry->toPoint();
369 :
370 5 : VSIFPrintfL(fp, " %s", d2str(poPoint->getX()));
371 5 : VSIFPrintfL(fp, " %s", d2str(poPoint->getY()));
372 : }
373 : // 3D Point
374 41 : else if (poGeometry->getGeometryType() == wkbPoint25D)
375 : {
376 3 : OGRPoint *poPoint = poGeometry->toPoint();
377 :
378 3 : VSIFPrintfL(fp, " %s", d2str(poPoint->getX()));
379 3 : VSIFPrintfL(fp, " %s", d2str(poPoint->getY()));
380 3 : VSIFPrintfL(fp, " %s", d2str(poPoint->getZ()));
381 : }
382 : }
383 : }
384 :
385 : // Write all fields.
386 436 : for (int iField = 0; iField < poFeatureDefn->GetFieldCount(); iField++)
387 : {
388 352 : if (poFeature->IsFieldSetAndNotNull(iField))
389 : {
390 253 : const char *pszRaw = poFeature->GetFieldAsString(iField);
391 253 : if (poFeatureDefn->GetFieldDefn(iField)->GetType() == OFTString)
392 : {
393 : // Interlis 1 encoding is ISO 8859-1 (Latin1) -> Recode from
394 : // UTF-8
395 : char *pszString =
396 58 : CPLRecode(pszRaw, CPL_ENC_UTF8, CPL_ENC_ISO8859_1);
397 : // Replace spaces
398 234 : for (size_t i = 0; pszString[i] != '\0'; i++)
399 : {
400 176 : if (pszString[i] == ' ')
401 2 : pszString[i] = '_';
402 : }
403 58 : VSIFPrintfL(fp, " %s", pszString);
404 58 : CPLFree(pszString);
405 : }
406 : else
407 : {
408 195 : VSIFPrintfL(fp, " %s", pszRaw);
409 : }
410 : }
411 : else
412 : {
413 99 : VSIFPrintfL(fp, " @");
414 : }
415 : }
416 84 : VSIFPrintfL(fp, "\n");
417 :
418 : // Write out Geometry
419 84 : if (poFeature->GetGeometryRef() != nullptr)
420 : {
421 46 : GeometryAppend(poFeature->GetGeometryRef());
422 : }
423 :
424 84 : return OGRERR_NONE;
425 : }
426 :
427 : /************************************************************************/
428 : /* TestCapability() */
429 : /************************************************************************/
430 :
431 107 : int OGRILI1Layer::TestCapability(CPL_UNUSED const char *pszCap)
432 : {
433 107 : if (EQUAL(pszCap, OLCCurveGeometries))
434 37 : return TRUE;
435 70 : if (EQUAL(pszCap, OLCZGeometries))
436 0 : return TRUE;
437 :
438 70 : if (EQUAL(pszCap, OLCCreateField) || EQUAL(pszCap, OLCSequentialWrite))
439 32 : return poDS->GetTransferFile() != nullptr;
440 :
441 38 : return FALSE;
442 : }
443 :
444 : /************************************************************************/
445 : /* CreateField() */
446 : /************************************************************************/
447 :
448 112 : OGRErr OGRILI1Layer::CreateField(const OGRFieldDefn *poField,
449 : int /* bApproxOK */)
450 : {
451 112 : poFeatureDefn->AddFieldDefn(poField);
452 :
453 112 : return OGRERR_NONE;
454 : }
455 :
456 : /************************************************************************/
457 : /* Internal routines */
458 : /************************************************************************/
459 :
460 36 : void OGRILI1Layer::JoinGeomLayers()
461 : {
462 36 : bGeomsJoined = true;
463 36 : bool bResetConfigOption = false;
464 36 : if (EQUAL(CPLGetConfigOption("OGR_ARC_STEPSIZE", ""), ""))
465 : {
466 36 : bResetConfigOption = true;
467 36 : CPLSetThreadLocalConfigOption("OGR_ARC_STEPSIZE", "0.96");
468 : }
469 :
470 49 : for (GeomFieldInfos::const_iterator it = oGeomFieldInfos.begin();
471 62 : it != oGeomFieldInfos.end(); ++it)
472 : {
473 13 : OGRFeatureDefn *geomFeatureDefn = it->second.GetGeomTableDefnRef();
474 13 : if (geomFeatureDefn)
475 : {
476 14 : CPLDebug("OGR_ILI", "Join geometry table %s of field '%s'",
477 14 : geomFeatureDefn->GetName(), it->first.c_str());
478 : OGRILI1Layer *poGeomLayer =
479 7 : poDS->GetLayerByName(geomFeatureDefn->GetName());
480 : const int nGeomFieldIndex =
481 7 : GetLayerDefn()->GetGeomFieldIndex(it->first.c_str());
482 7 : if (it->second.iliGeomType == "Surface")
483 : {
484 5 : JoinSurfaceLayer(poGeomLayer, nGeomFieldIndex);
485 : }
486 2 : else if (it->second.iliGeomType == "Area")
487 : {
488 4 : CPLString pointField = it->first + "__Point";
489 : const int nPointFieldIndex =
490 2 : GetLayerDefn()->GetGeomFieldIndex(pointField.c_str());
491 2 : PolygonizeAreaLayer(poGeomLayer, nGeomFieldIndex,
492 : nPointFieldIndex);
493 : }
494 : }
495 : }
496 :
497 36 : if (bResetConfigOption)
498 36 : CPLSetThreadLocalConfigOption("OGR_ARC_STEPSIZE", nullptr);
499 36 : }
500 :
501 5 : void OGRILI1Layer::JoinSurfaceLayer(OGRILI1Layer *poSurfaceLineLayer,
502 : int nSurfaceFieldIndex)
503 : {
504 5 : CPLDebug("OGR_ILI", "Joining surface layer %s with geometries",
505 5 : GetLayerDefn()->GetName());
506 : OGRwkbGeometryType geomType =
507 5 : GetLayerDefn()->GetGeomFieldDefn(nSurfaceFieldIndex)->GetType();
508 :
509 10 : std::map<OGRFeature *, std::vector<OGRCurve *>> oMapFeatureToGeomSet;
510 :
511 5 : poSurfaceLineLayer->ResetReading();
512 :
513 : // First map: for each target curvepolygon, find all belonging curves
514 24 : while (OGRFeature *linefeature = poSurfaceLineLayer->GetNextFeatureRef())
515 : {
516 : // OBJE entries with same _RefTID are polygon rings of same feature
517 19 : OGRFeature *feature = nullptr;
518 19 : if (poFeatureDefn->GetFieldDefn(0)->GetType() == OFTString)
519 : {
520 19 : feature = GetFeatureRef(linefeature->GetFieldAsString(1));
521 : }
522 : else
523 : {
524 0 : GIntBig reftid = linefeature->GetFieldAsInteger64(1);
525 0 : feature = GetFeatureRef(reftid);
526 : }
527 19 : if (feature)
528 : {
529 19 : OGRGeometry *poGeom = linefeature->GetGeomFieldRef(0);
530 19 : OGRMultiCurve *curves = dynamic_cast<OGRMultiCurve *>(poGeom);
531 19 : if (curves)
532 : {
533 38 : for (auto &&curve : curves)
534 : {
535 19 : if (!curve->IsEmpty())
536 19 : oMapFeatureToGeomSet[feature].push_back(curve);
537 : }
538 : }
539 : }
540 : else
541 : {
542 0 : CPLError(CE_Warning, CPLE_AppDefined,
543 : "Couldn't join feature FID " CPL_FRMT_GIB,
544 : linefeature->GetFieldAsInteger64(1));
545 : }
546 19 : }
547 :
548 : // Now for each target polygon, assemble the curves together.
549 : std::map<OGRFeature *, std::vector<OGRCurve *>>::const_iterator oIter =
550 5 : oMapFeatureToGeomSet.begin();
551 18 : for (; oIter != oMapFeatureToGeomSet.end(); ++oIter)
552 : {
553 13 : OGRFeature *feature = oIter->first;
554 26 : std::vector<OGRCurve *> oCurves = oIter->second;
555 :
556 26 : std::vector<OGRCurve *> oSetDestCurves;
557 13 : double dfLargestArea = 0.0;
558 13 : OGRCurve *poLargestCurve = nullptr;
559 : while (true)
560 : {
561 28 : std::vector<OGRCurve *>::iterator oIterCurves = oCurves.begin();
562 28 : if (oIterCurves == oCurves.end())
563 13 : break;
564 :
565 30 : OGRPoint endPointCC;
566 15 : OGRCompoundCurve *poCC = new OGRCompoundCurve();
567 :
568 15 : bool bFirst = true;
569 : while (true)
570 : {
571 19 : bool bNewCurveAdded = false;
572 19 : const double dfEps = 1e-14;
573 20 : for (oIterCurves = oCurves.begin();
574 21 : oIterCurves != oCurves.end(); ++oIterCurves)
575 : {
576 20 : OGRCurve *curve = *oIterCurves;
577 20 : OGRPoint startPoint;
578 20 : OGRPoint endPoint;
579 20 : curve->StartPoint(&startPoint);
580 20 : curve->EndPoint(&endPoint);
581 25 : if (bFirst ||
582 5 : (fabs(startPoint.getX() - endPointCC.getX()) < dfEps &&
583 4 : fabs(startPoint.getY() - endPointCC.getY()) < dfEps))
584 : {
585 19 : bFirst = false;
586 :
587 19 : curve->EndPoint(&endPointCC);
588 :
589 : const OGRwkbGeometryType eCurveType =
590 19 : wkbFlatten(curve->getGeometryType());
591 19 : if (eCurveType == wkbCompoundCurve)
592 : {
593 : OGRCompoundCurve *poCCSub =
594 11 : curve->toCompoundCurve();
595 24 : for (auto &&subCurve : poCCSub)
596 : {
597 13 : poCC->addCurve(subCurve);
598 : }
599 : }
600 : else
601 : {
602 8 : poCC->addCurve(curve);
603 : }
604 19 : oCurves.erase(oIterCurves);
605 19 : bNewCurveAdded = true;
606 19 : break;
607 : }
608 1 : else if (fabs(endPoint.getX() - endPointCC.getX()) <
609 1 : dfEps &&
610 0 : fabs(endPoint.getY() - endPointCC.getY()) < dfEps)
611 : {
612 0 : curve->StartPoint(&endPointCC);
613 :
614 : const OGRwkbGeometryType eCurveType =
615 0 : wkbFlatten(curve->getGeometryType());
616 0 : if (eCurveType == wkbLineString ||
617 : eCurveType == wkbCircularString)
618 : {
619 : OGRSimpleCurve *poSC =
620 0 : curve->toSimpleCurve()->clone();
621 0 : poSC->reversePoints();
622 0 : poCC->addCurveDirectly(poSC);
623 : }
624 0 : else if (eCurveType == wkbCompoundCurve)
625 : {
626 : // Reverse the order of the elements of the
627 : // compound curve
628 : OGRCompoundCurve *poCCSub =
629 0 : curve->toCompoundCurve();
630 0 : for (int i = poCCSub->getNumCurves() - 1; i >= 0;
631 : --i)
632 : {
633 : OGRSimpleCurve *poSC = poCCSub->getCurve(i)
634 0 : ->toSimpleCurve()
635 0 : ->clone();
636 0 : poSC->reversePoints();
637 0 : poCC->addCurveDirectly(poSC);
638 : }
639 : }
640 :
641 0 : oCurves.erase(oIterCurves);
642 0 : bNewCurveAdded = true;
643 0 : break;
644 : }
645 : }
646 19 : if (!bNewCurveAdded || oCurves.empty() || poCC->get_IsClosed())
647 15 : break;
648 4 : }
649 :
650 15 : if (!poCC->get_IsClosed())
651 : {
652 0 : char *pszJSon = poCC->exportToJson();
653 0 : CPLError(CE_Warning, CPLE_AppDefined,
654 : "A ring %s for feature " CPL_FRMT_GIB " in layer %s "
655 : "was not closed. Dropping it",
656 0 : pszJSon, feature->GetFID(), GetName());
657 0 : delete poCC;
658 0 : CPLFree(pszJSon);
659 : }
660 : else
661 : {
662 15 : double dfArea = poCC->get_Area();
663 15 : if (dfArea >= dfLargestArea)
664 : {
665 15 : dfLargestArea = dfArea;
666 15 : poLargestCurve = poCC;
667 : }
668 15 : oSetDestCurves.push_back(poCC);
669 : }
670 15 : }
671 :
672 : // Now build the final polygon by first inserting the largest ring.
673 : OGRCurvePolygon *poPoly =
674 13 : (geomType == wkbPolygon) ? new OGRPolygon() : new OGRCurvePolygon();
675 13 : if (poLargestCurve)
676 : {
677 : std::vector<OGRCurve *>::iterator oIterCurves =
678 13 : oSetDestCurves.begin();
679 15 : for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
680 : {
681 15 : OGRCurve *poCurve = *oIterCurves;
682 15 : if (poCurve == poLargestCurve)
683 : {
684 13 : oSetDestCurves.erase(oIterCurves);
685 13 : break;
686 : }
687 : }
688 :
689 13 : if (geomType == wkbPolygon)
690 : {
691 6 : poLargestCurve = OGRCurve::CastToLinearRing(poLargestCurve);
692 : }
693 13 : OGRErr error = poPoly->addRingDirectly(poLargestCurve);
694 13 : if (error != OGRERR_NONE)
695 : {
696 0 : char *pszJSon = poLargestCurve->exportToJson();
697 0 : CPLError(CE_Warning, CPLE_AppDefined,
698 : "Cannot add ring %s to feature " CPL_FRMT_GIB
699 : " in layer %s",
700 0 : pszJSon, feature->GetFID(), GetName());
701 0 : CPLFree(pszJSon);
702 : }
703 :
704 13 : oIterCurves = oSetDestCurves.begin();
705 15 : for (; oIterCurves != oSetDestCurves.end(); ++oIterCurves)
706 : {
707 2 : OGRCurve *poCurve = *oIterCurves;
708 2 : if (geomType == wkbPolygon)
709 : {
710 1 : poCurve = OGRCurve::CastToLinearRing(poCurve);
711 : }
712 2 : error = poPoly->addRingDirectly(poCurve);
713 2 : if (error != OGRERR_NONE)
714 : {
715 0 : char *pszJSon = poCurve->exportToJson();
716 0 : CPLError(CE_Warning, CPLE_AppDefined,
717 : "Cannot add ring %s to feature " CPL_FRMT_GIB
718 : " in layer %s",
719 0 : pszJSon, feature->GetFID(), GetName());
720 0 : CPLFree(pszJSon);
721 : }
722 : }
723 : }
724 :
725 13 : feature->SetGeomFieldDirectly(nSurfaceFieldIndex, poPoly);
726 : }
727 :
728 5 : ResetReading();
729 5 : }
730 :
731 2 : OGRMultiPolygon *OGRILI1Layer::Polygonize(OGRGeometryCollection *poLines,
732 : bool
733 : #if defined(HAVE_GEOS)
734 : fix_crossing_lines
735 : #endif
736 : )
737 : {
738 2 : if (poLines->getNumGeometries() == 0)
739 : {
740 0 : return new OGRMultiPolygon();
741 : }
742 :
743 : #if defined(HAVE_GEOS)
744 2 : GEOSGeom *ahInGeoms = nullptr;
745 2 : OGRGeometryCollection *poNoncrossingLines = poLines;
746 2 : GEOSGeom hResultGeom = nullptr;
747 2 : OGRGeometry *poMP = nullptr;
748 :
749 2 : if (fix_crossing_lines && poLines->getNumGeometries() > 0)
750 : {
751 0 : CPLDebug("OGR_ILI", "Fixing crossing lines");
752 : // A union of the geometry collection with one line fixes
753 : // invalid geometries.
754 0 : OGRGeometry *poUnion = poLines->Union(poLines->getGeometryRef(0));
755 0 : if (poUnion != nullptr)
756 : {
757 0 : if (wkbFlatten(poUnion->getGeometryType()) ==
758 0 : wkbGeometryCollection ||
759 0 : wkbFlatten(poUnion->getGeometryType()) == wkbMultiLineString)
760 : {
761 0 : poNoncrossingLines =
762 0 : dynamic_cast<OGRGeometryCollection *>(poUnion);
763 0 : CPLDebug("OGR_ILI", "Fixed lines: %d",
764 0 : poNoncrossingLines->getNumGeometries() -
765 0 : poLines->getNumGeometries());
766 : }
767 : else
768 : {
769 0 : delete poUnion;
770 : }
771 : }
772 : }
773 :
774 2 : GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
775 :
776 : ahInGeoms = static_cast<GEOSGeom *>(
777 2 : CPLCalloc(sizeof(void *), poNoncrossingLines->getNumGeometries()));
778 10 : for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
779 16 : ahInGeoms[i] =
780 8 : poNoncrossingLines->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
781 :
782 4 : hResultGeom = GEOSPolygonize_r(hGEOSCtxt, ahInGeoms,
783 2 : poNoncrossingLines->getNumGeometries());
784 :
785 10 : for (int i = 0; i < poNoncrossingLines->getNumGeometries(); i++)
786 8 : GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
787 2 : CPLFree(ahInGeoms);
788 2 : if (poNoncrossingLines != poLines)
789 0 : delete poNoncrossingLines;
790 :
791 2 : if (hResultGeom == nullptr)
792 : {
793 0 : OGRGeometry::freeGEOSContext(hGEOSCtxt);
794 0 : return new OGRMultiPolygon();
795 : }
796 :
797 2 : poMP = OGRGeometryFactory::createFromGEOS(hGEOSCtxt, hResultGeom);
798 :
799 2 : GEOSGeom_destroy_r(hGEOSCtxt, hResultGeom);
800 2 : OGRGeometry::freeGEOSContext(hGEOSCtxt);
801 :
802 2 : poMP = OGRGeometryFactory::forceToMultiPolygon(poMP);
803 2 : if (poMP && wkbFlatten(poMP->getGeometryType()) == wkbMultiPolygon)
804 2 : return dynamic_cast<OGRMultiPolygon *>(poMP);
805 :
806 0 : delete poMP;
807 0 : return new OGRMultiPolygon();
808 :
809 : #else
810 : return new OGRMultiPolygon();
811 : #endif
812 : }
813 :
814 2 : void OGRILI1Layer::PolygonizeAreaLayer(OGRILI1Layer *poAreaLineLayer,
815 : int
816 : #if defined(HAVE_GEOS)
817 : nAreaFieldIndex
818 : #endif
819 : ,
820 : int
821 : #if defined(HAVE_GEOS)
822 : nPointFieldIndex
823 : #endif
824 : )
825 : {
826 : // add all lines from poAreaLineLayer to collection
827 2 : OGRGeometryCollection *gc = new OGRGeometryCollection();
828 2 : poAreaLineLayer->ResetReading();
829 10 : while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
830 8 : gc->addGeometry(feature->GetGeometryRef());
831 :
832 : // polygonize lines
833 4 : CPLDebug("OGR_ILI", "Polygonizing layer %s with %d multilines",
834 2 : poAreaLineLayer->GetLayerDefn()->GetName(),
835 : gc->getNumGeometries());
836 2 : poAreaLineLayer = nullptr;
837 2 : OGRMultiPolygon *polys = Polygonize(gc, false);
838 2 : CPLDebug("OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
839 2 : if (polys->getNumGeometries() != GetFeatureCount())
840 : {
841 0 : CPLDebug("OGR_ILI", "Feature count of layer %s: " CPL_FRMT_GIB,
842 0 : GetLayerDefn()->GetName(), GetFeatureCount());
843 0 : CPLDebug("OGR_ILI", "Polygonizing again with crossing line fix");
844 0 : delete polys;
845 0 : polys = Polygonize(gc, true); // try again with crossing line fix
846 0 : CPLDebug("OGR_ILI", "Resulting polygons: %d",
847 : polys->getNumGeometries());
848 : }
849 2 : delete gc;
850 :
851 : // associate polygon feature with data row according to centroid
852 : #if defined(HAVE_GEOS)
853 4 : OGRPolygon emptyPoly;
854 :
855 2 : CPLDebug("OGR_ILI", "Associating layer %s with area polygons",
856 2 : GetLayerDefn()->GetName());
857 : GEOSGeom *ahInGeoms = static_cast<GEOSGeom *>(
858 2 : CPLCalloc(sizeof(void *), polys->getNumGeometries()));
859 2 : GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
860 8 : for (int i = 0; i < polys->getNumGeometries(); i++)
861 : {
862 6 : ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);
863 6 : if (!GEOSisValid_r(hGEOSCtxt, ahInGeoms[i]))
864 0 : ahInGeoms[i] = nullptr;
865 : }
866 8 : for (int nFidx = 0; nFidx < nFeatures; nFidx++)
867 : {
868 6 : OGRFeature *feature = papoFeatures[nFidx];
869 6 : OGRGeometry *geomRef = feature->GetGeomFieldRef(nPointFieldIndex);
870 6 : if (!geomRef)
871 : {
872 0 : continue;
873 : }
874 : GEOSGeom point =
875 6 : reinterpret_cast<GEOSGeom>(geomRef->exportToGEOS(hGEOSCtxt));
876 :
877 6 : int i = 0;
878 12 : for (; i < polys->getNumGeometries(); i++)
879 : {
880 12 : if (ahInGeoms[i] && GEOSWithin_r(hGEOSCtxt, point, ahInGeoms[i]))
881 : {
882 6 : feature->SetGeomField(nAreaFieldIndex,
883 6 : polys->getGeometryRef(i));
884 6 : break;
885 : }
886 : }
887 6 : if (i == polys->getNumGeometries())
888 : {
889 0 : CPLDebug("OGR_ILI", "Association between area and point failed.");
890 0 : feature->SetGeometry(&emptyPoly);
891 : }
892 6 : GEOSGeom_destroy_r(hGEOSCtxt, point);
893 : }
894 8 : for (int i = 0; i < polys->getNumGeometries(); i++)
895 6 : GEOSGeom_destroy_r(hGEOSCtxt, ahInGeoms[i]);
896 2 : CPLFree(ahInGeoms);
897 2 : OGRGeometry::freeGEOSContext(hGEOSCtxt);
898 : #endif
899 2 : delete polys;
900 2 : }
901 :
902 : /************************************************************************/
903 : /* GetDataset() */
904 : /************************************************************************/
905 :
906 17 : GDALDataset *OGRILI1Layer::GetDataset()
907 : {
908 17 : return poDS;
909 : }
|
