1 : /******************************************************************************
2 : * $Id: ogrili1layer.cpp 16925 2009-05-03 14:16:49Z rouault $
3 : *
4 : * Project: Interlis 1 Translator
5 : * Purpose: Implements OGRILI1Layer class.
6 : * Author: Pirmin Kalberer, Sourcepole AG
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Pirmin Kalberer, Sourcepole AG
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 "ogr_ili1.h"
31 : #include "cpl_conv.h"
32 : #include "cpl_string.h"
33 : #include "ogr_geos.h"
34 :
35 : CPL_CVSID("$Id: ogrili1layer.cpp 16925 2009-05-03 14:16:49Z rouault $");
36 :
37 : /************************************************************************/
38 : /* OGRILI1Layer() */
39 : /************************************************************************/
40 :
41 : OGRILI1Layer::OGRILI1Layer( const char * pszName,
42 : OGRSpatialReference *poSRSIn, int bWriterIn,
43 : OGRwkbGeometryType eReqType,
44 0 : OGRILI1DataSource *poDSIn )
45 :
46 : {
47 0 : if( poSRSIn == NULL )
48 0 : poSRS = NULL;
49 : else
50 0 : poSRS = poSRSIn->Clone();
51 :
52 0 : poDS = poDSIn;
53 :
54 0 : poFeatureDefn = new OGRFeatureDefn( pszName );
55 0 : poFeatureDefn->Reference();
56 0 : poFeatureDefn->SetGeomType( eReqType );
57 :
58 0 : nFeatures = 0;
59 0 : papoFeatures = NULL;
60 0 : nFeatureIdx = 0;
61 :
62 0 : poSurfacePolyLayer = 0;
63 0 : poAreaLineLayer = 0;
64 :
65 0 : bWriter = bWriterIn;
66 0 : }
67 :
68 : /************************************************************************/
69 : /* ~OGRILI1Layer() */
70 : /************************************************************************/
71 :
72 0 : OGRILI1Layer::~OGRILI1Layer()
73 : {
74 : int i;
75 :
76 0 : for(i=0;i<nFeatures;i++)
77 : {
78 0 : delete papoFeatures[i];
79 : }
80 0 : CPLFree(papoFeatures);
81 :
82 0 : if( poFeatureDefn )
83 0 : poFeatureDefn->Release();
84 :
85 0 : if( poSRS != NULL )
86 0 : poSRS->Release();
87 0 : }
88 :
89 :
90 0 : OGRErr OGRILI1Layer::AddFeature (OGRFeature *poFeature)
91 : {
92 0 : nFeatures++;
93 :
94 : papoFeatures = (OGRFeature **)
95 0 : CPLRealloc( papoFeatures, sizeof(void*) * nFeatures );
96 :
97 0 : papoFeatures[nFeatures-1] = poFeature;
98 :
99 0 : return OGRERR_NONE;
100 : }
101 :
102 : /************************************************************************/
103 : /* ResetReading() */
104 : /************************************************************************/
105 :
106 0 : void OGRILI1Layer::ResetReading(){
107 0 : nFeatureIdx = 0;
108 0 : }
109 :
110 : /************************************************************************/
111 : /* GetNextFeature() */
112 : /************************************************************************/
113 :
114 0 : OGRFeature *OGRILI1Layer::GetNextFeature()
115 : {
116 : OGRFeature *poFeature;
117 :
118 0 : if (poSurfacePolyLayer != 0) JoinSurfaceLayer();
119 0 : if (poAreaLineLayer != 0) PolygonizeAreaLayer();
120 :
121 0 : while(nFeatureIdx < nFeatures)
122 : {
123 0 : poFeature = GetNextFeatureRef();
124 0 : if (poFeature)
125 0 : return poFeature->Clone();
126 : }
127 0 : return NULL;
128 : }
129 :
130 0 : OGRFeature *OGRILI1Layer::GetNextFeatureRef() {
131 0 : OGRFeature *poFeature = NULL;
132 0 : if (nFeatureIdx < nFeatures)
133 : {
134 0 : poFeature = papoFeatures[nFeatureIdx++];
135 : //apply filters
136 0 : if( (m_poFilterGeom == NULL
137 : || FilterGeometry( poFeature->GetGeometryRef() ) )
138 : && (m_poAttrQuery == NULL
139 : || m_poAttrQuery->Evaluate( poFeature )) )
140 0 : return poFeature;
141 : }
142 0 : return NULL;
143 : }
144 :
145 : /************************************************************************/
146 : /* GetFeatureRef() */
147 : /************************************************************************/
148 :
149 0 : OGRFeature *OGRILI1Layer::GetFeatureRef( long nFID )
150 :
151 : {
152 : OGRFeature *poFeature;
153 :
154 0 : ResetReading();
155 0 : while( (poFeature = GetNextFeatureRef()) != NULL )
156 : {
157 0 : if( poFeature->GetFID() == nFID )
158 0 : return poFeature;
159 : }
160 :
161 0 : return NULL;
162 : }
163 :
164 : /************************************************************************/
165 : /* GetFeatureCount() */
166 : /************************************************************************/
167 :
168 0 : int OGRILI1Layer::GetFeatureCount( int bForce )
169 : {
170 0 : if (m_poFilterGeom == NULL && m_poAttrQuery == NULL)
171 : {
172 0 : return nFeatures;
173 : }
174 : else
175 : {
176 0 : return OGRLayer::GetFeatureCount(bForce);
177 : }
178 : }
179 :
180 0 : static char* d2str(double val)
181 : {
182 : static char strbuf[255];
183 0 : if( val == (int) val )
184 0 : sprintf( strbuf, "%d", (int) val );
185 0 : else if( fabs(val) < 370 )
186 0 : sprintf( strbuf, "%.16g", val );
187 0 : else if( fabs(val) > 100000000.0 )
188 0 : sprintf( strbuf, "%.16g", val );
189 : else
190 0 : sprintf( strbuf, "%.3f", val );
191 0 : return strbuf;
192 : }
193 :
194 0 : static void AppendCoordinateList( OGRLineString *poLine, OGRILI1DataSource *poDS)
195 : {
196 0 : int b3D = (poLine->getGeometryType() & wkb25DBit);
197 :
198 0 : for( int iPoint = 0; iPoint < poLine->getNumPoints(); iPoint++ )
199 : {
200 0 : if (iPoint == 0) VSIFPrintf( poDS->GetTransferFile(), "STPT" );
201 0 : else VSIFPrintf( poDS->GetTransferFile(), "LIPT" );
202 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poLine->getX(iPoint)) );
203 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poLine->getY(iPoint)) );
204 0 : if (b3D) VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poLine->getZ(iPoint)) );
205 0 : VSIFPrintf( poDS->GetTransferFile(), "\n" );
206 : }
207 0 : VSIFPrintf( poDS->GetTransferFile(), "ELIN\n" );
208 0 : }
209 :
210 0 : int OGRILI1Layer::GeometryAppend( OGRGeometry *poGeometry )
211 : {
212 : /* -------------------------------------------------------------------- */
213 : /* 2D Point */
214 : /* -------------------------------------------------------------------- */
215 0 : if( poGeometry->getGeometryType() == wkbPoint )
216 : {
217 : /* embedded in from non-geometry fields */
218 : }
219 : /* -------------------------------------------------------------------- */
220 : /* 3D Point */
221 : /* -------------------------------------------------------------------- */
222 0 : else if( poGeometry->getGeometryType() == wkbPoint25D )
223 : {
224 : /* embedded in from non-geometry fields */
225 : }
226 :
227 : /* -------------------------------------------------------------------- */
228 : /* LineString and LinearRing */
229 : /* -------------------------------------------------------------------- */
230 0 : else if( poGeometry->getGeometryType() == wkbLineString
231 : || poGeometry->getGeometryType() == wkbLineString25D )
232 : {
233 0 : AppendCoordinateList( (OGRLineString *) poGeometry, poDS );
234 : }
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Polygon */
238 : /* -------------------------------------------------------------------- */
239 0 : else if( poGeometry->getGeometryType() == wkbPolygon
240 : || poGeometry->getGeometryType() == wkbPolygon25D )
241 : {
242 0 : OGRPolygon *poPolygon = (OGRPolygon *) poGeometry;
243 :
244 0 : if( poPolygon->getExteriorRing() != NULL )
245 : {
246 0 : if( !GeometryAppend( poPolygon->getExteriorRing() ) )
247 0 : return FALSE;
248 : }
249 :
250 0 : for( int iRing = 0; iRing < poPolygon->getNumInteriorRings(); iRing++ )
251 : {
252 0 : OGRLinearRing *poRing = poPolygon->getInteriorRing(iRing);
253 :
254 0 : if( !GeometryAppend( poRing ) )
255 0 : return FALSE;
256 : }
257 : }
258 :
259 : /* -------------------------------------------------------------------- */
260 : /* MultiPolygon */
261 : /* -------------------------------------------------------------------- */
262 0 : else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon
263 : || wkbFlatten(poGeometry->getGeometryType()) == wkbMultiLineString
264 : || wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPoint
265 : || wkbFlatten(poGeometry->getGeometryType()) == wkbGeometryCollection )
266 : {
267 0 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeometry;
268 : int iMember;
269 :
270 0 : if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon )
271 : {
272 : }
273 0 : else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiLineString )
274 : {
275 : }
276 0 : else if( wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPoint )
277 : {
278 : }
279 : else
280 : {
281 : }
282 :
283 0 : for( iMember = 0; iMember < poGC->getNumGeometries(); iMember++)
284 : {
285 0 : OGRGeometry *poMember = poGC->getGeometryRef( iMember );
286 :
287 0 : if( !GeometryAppend( poMember ) )
288 0 : return FALSE;
289 : }
290 :
291 : }
292 :
293 : else
294 0 : return FALSE;
295 :
296 0 : return TRUE;
297 : }
298 :
299 : /************************************************************************/
300 : /* CreateFeature() */
301 : /************************************************************************/
302 :
303 0 : OGRErr OGRILI1Layer::CreateFeature( OGRFeature *poFeature ) {
304 : static long tid = -1; //system generated TID (must be unique within table)
305 0 : VSIFPrintf( poDS->GetTransferFile(), "OBJE" );
306 :
307 0 : if ( !EQUAL(poFeatureDefn->GetFieldDefn(0)->GetNameRef(), "TID") )
308 : {
309 : //Input is not generated from an Interlis 1 source
310 0 : if (poFeature->GetFID() != OGRNullFID)
311 0 : tid = poFeature->GetFID();
312 : else
313 0 : ++tid;
314 0 : VSIFPrintf( poDS->GetTransferFile(), " %ld", tid );
315 : //Embedded geometry
316 0 : if( poFeature->GetGeometryRef() != NULL )
317 : {
318 0 : OGRGeometry *poGeometry = poFeature->GetGeometryRef();
319 : // 2D Point
320 0 : if( poGeometry->getGeometryType() == wkbPoint )
321 : {
322 0 : OGRPoint *poPoint = (OGRPoint *) poGeometry;
323 :
324 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poPoint->getX()) );
325 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poPoint->getY()) );
326 : }
327 : // 3D Point
328 0 : else if( poGeometry->getGeometryType() == wkbPoint25D )
329 : {
330 0 : OGRPoint *poPoint = (OGRPoint *) poGeometry;
331 :
332 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poPoint->getX()) );
333 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poPoint->getY()) );
334 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", d2str(poPoint->getZ()) );
335 : }
336 : }
337 : }
338 :
339 : // Write all fields.
340 0 : for(int iField = 0; iField < poFeatureDefn->GetFieldCount(); iField++ )
341 : {
342 0 : if ( !EQUAL(poFeatureDefn->GetFieldDefn(iField)->GetNameRef(), "ILI_Geometry") )
343 : {
344 0 : if ( poFeature->IsFieldSet( iField ) )
345 : {
346 0 : const char *pszRaw = poFeature->GetFieldAsString( iField );
347 0 : VSIFPrintf( poDS->GetTransferFile(), " %s", pszRaw );
348 : }
349 : else
350 : {
351 0 : VSIFPrintf( poDS->GetTransferFile(), " @" );
352 : }
353 : }
354 : }
355 0 : VSIFPrintf( poDS->GetTransferFile(), "\n" );
356 :
357 : // Write out Geometry
358 0 : if( poFeature->GetGeometryRef() != NULL )
359 : {
360 0 : if (EQUAL(poFeatureDefn->GetFieldDefn(poFeatureDefn->GetFieldCount()-1)->GetNameRef(), "ILI_Geometry"))
361 : {
362 : //Write original ILI geometry
363 0 : VSIFPrintf( poDS->GetTransferFile(), "%s", poFeature->GetFieldAsString( poFeatureDefn->GetFieldCount()-1 ) );
364 : }
365 : else
366 : {
367 : //Convert to ILI geometry
368 0 : GeometryAppend(poFeature->GetGeometryRef());
369 : }
370 : }
371 :
372 0 : return OGRERR_NONE;
373 : }
374 :
375 : /************************************************************************/
376 : /* TestCapability() */
377 : /************************************************************************/
378 :
379 0 : int OGRILI1Layer::TestCapability( const char * pszCap ) {
380 0 : return FALSE;
381 : }
382 :
383 : /************************************************************************/
384 : /* CreateField() */
385 : /************************************************************************/
386 :
387 0 : OGRErr OGRILI1Layer::CreateField( OGRFieldDefn *poField, int bApproxOK ) {
388 0 : poFeatureDefn->AddFieldDefn( poField );
389 :
390 0 : return OGRERR_NONE;
391 : }
392 :
393 : /************************************************************************/
394 : /* GetSpatialRef() */
395 : /************************************************************************/
396 :
397 0 : OGRSpatialReference *OGRILI1Layer::GetSpatialRef() {
398 0 : return poSRS;
399 : }
400 :
401 :
402 : /************************************************************************/
403 : /* Internal routines */
404 : /************************************************************************/
405 :
406 0 : void OGRILI1Layer::SetSurfacePolyLayer(OGRILI1Layer *poSurfacePolyLayerIn) {
407 0 : poSurfacePolyLayer = poSurfacePolyLayerIn;
408 0 : }
409 :
410 0 : void OGRILI1Layer::JoinSurfaceLayer()
411 : {
412 0 : if (poSurfacePolyLayer == 0) return;
413 :
414 0 : CPLDebug( "OGR_ILI", "Joining surface layer %s with geometries", GetLayerDefn()->GetName());
415 0 : GetLayerDefn()->SetGeomType(poSurfacePolyLayer->GetLayerDefn()->GetGeomType());
416 0 : ResetReading();
417 0 : while (OGRFeature *feature = GetNextFeatureRef())
418 : {
419 0 : OGRFeature *polyfeature = poSurfacePolyLayer->GetFeatureRef(feature->GetFID());
420 0 : if (polyfeature) {
421 0 : feature->SetGeometry(polyfeature->GetGeometryRef());
422 : }
423 : }
424 :
425 0 : ResetReading();
426 0 : poSurfacePolyLayer = 0;
427 : }
428 :
429 0 : void OGRILI1Layer::SetAreaLayers(OGRILI1Layer *poReferenceLayer, OGRILI1Layer *poAreaLineLayerIn) {
430 0 : poAreaReferenceLayer = poReferenceLayer;
431 0 : poAreaLineLayer = poAreaLineLayerIn;
432 0 : }
433 :
434 0 : OGRMultiPolygon* OGRILI1Layer::Polygonize( OGRGeometryCollection* poLines, bool fix_crossing_lines )
435 : {
436 0 : OGRMultiPolygon *poPolygon = new OGRMultiPolygon();
437 :
438 0 : if (poLines->getNumGeometries() == 0) return poPolygon;
439 :
440 : #if defined(HAVE_GEOS)
441 0 : GEOSGeom *ahInGeoms = NULL;
442 0 : int i = 0;
443 0 : OGRGeometryCollection *poNoncrossingLines = poLines;
444 0 : GEOSGeom hResultGeom = NULL;
445 0 : OGRGeometry *poMP = NULL;
446 :
447 0 : if (fix_crossing_lines && poLines->getNumGeometries() > 0)
448 : {
449 : #if (GEOS_VERSION_MAJOR >= 3)
450 : CPLDebug( "OGR_ILI", "Fixing crossing lines");
451 : //A union of the geometry collection with one line fixes invalid geometries
452 : poNoncrossingLines = (OGRGeometryCollection*)poLines->Union(poLines->getGeometryRef(0));
453 : CPLDebug( "OGR_ILI", "Fixed lines: %d", poNoncrossingLines->getNumGeometries()-poLines->getNumGeometries());
454 : #else
455 : #warning Interlis 1 AREA cleanup disabled. Needs GEOS >= 3.0
456 : #endif
457 : }
458 :
459 0 : ahInGeoms = (GEOSGeom *) CPLCalloc(sizeof(void*),poNoncrossingLines->getNumGeometries());
460 0 : for( i = 0; i < poNoncrossingLines->getNumGeometries(); i++ )
461 0 : ahInGeoms[i] = poNoncrossingLines->getGeometryRef(i)->exportToGEOS();
462 :
463 : hResultGeom = GEOSPolygonize( ahInGeoms,
464 0 : poNoncrossingLines->getNumGeometries() );
465 :
466 0 : for( i = 0; i < poNoncrossingLines->getNumGeometries(); i++ )
467 0 : GEOSGeom_destroy( ahInGeoms[i] );
468 0 : CPLFree( ahInGeoms );
469 0 : if (poNoncrossingLines != poLines) delete poNoncrossingLines;
470 :
471 0 : if( hResultGeom == NULL )
472 0 : return NULL;
473 :
474 0 : poMP = OGRGeometryFactory::createFromGEOS( hResultGeom );
475 :
476 0 : GEOSGeom_destroy( hResultGeom );
477 :
478 0 : return (OGRMultiPolygon *) poMP;
479 :
480 : #endif
481 :
482 : return poPolygon;
483 : }
484 :
485 :
486 0 : void OGRILI1Layer::PolygonizeAreaLayer()
487 : {
488 0 : if (poAreaLineLayer == 0) return;
489 :
490 : //add all lines from poAreaLineLayer to collection
491 0 : OGRGeometryCollection *gc = new OGRGeometryCollection();
492 0 : poAreaLineLayer->ResetReading();
493 0 : while (OGRFeature *feature = poAreaLineLayer->GetNextFeatureRef())
494 0 : gc->addGeometry(feature->GetGeometryRef());
495 :
496 : //polygonize lines
497 0 : CPLDebug( "OGR_ILI", "Polygonizing layer %s with %d multilines", poAreaLineLayer->GetLayerDefn()->GetName(), gc->getNumGeometries());
498 0 : OGRMultiPolygon* polys = Polygonize( gc , false);
499 0 : CPLDebug( "OGR_ILI", "Resulting polygons: %d", polys->getNumGeometries());
500 0 : if (polys->getNumGeometries() != poAreaReferenceLayer->GetFeatureCount())
501 : {
502 0 : CPLDebug( "OGR_ILI", "Feature count of layer %s: %d", poAreaReferenceLayer->GetLayerDefn()->GetName(), GetFeatureCount());
503 0 : CPLDebug( "OGR_ILI", "Polygonizing again with crossing line fix");
504 0 : delete polys;
505 0 : polys = Polygonize( gc, true ); //try again with crossing line fix
506 : }
507 0 : delete gc;
508 :
509 : //associate polygon feature with data row according to centroid
510 : #if defined(HAVE_GEOS)
511 : int i;
512 0 : OGRPolygon emptyPoly;
513 0 : GEOSGeom *ahInGeoms = NULL;
514 :
515 0 : CPLDebug( "OGR_ILI", "Associating layer %s with area polygons", GetLayerDefn()->GetName());
516 0 : ahInGeoms = (GEOSGeom *) CPLCalloc(sizeof(void*),polys->getNumGeometries());
517 0 : for( i = 0; i < polys->getNumGeometries(); i++ )
518 : {
519 0 : ahInGeoms[i] = polys->getGeometryRef(i)->exportToGEOS();
520 0 : if (!GEOSisValid(ahInGeoms[i])) ahInGeoms[i] = NULL;
521 : }
522 0 : poAreaReferenceLayer->ResetReading();
523 0 : while (OGRFeature *feature = poAreaReferenceLayer->GetNextFeatureRef())
524 : {
525 0 : GEOSGeom point = (GEOSGeom)feature->GetGeometryRef()->exportToGEOS();
526 0 : for (i = 0; i < polys->getNumGeometries(); i++ )
527 : {
528 0 : if (ahInGeoms[i] && GEOSWithin(point, ahInGeoms[i]))
529 : {
530 0 : OGRFeature* areaFeature = feature->Clone();
531 0 : areaFeature->SetGeometry( polys->getGeometryRef(i) );
532 0 : AddFeature(areaFeature);
533 0 : break;
534 : }
535 : }
536 0 : if (i == polys->getNumGeometries())
537 : {
538 0 : CPLDebug( "OGR_ILI", "Association between area and point failed.");
539 0 : feature->SetGeometry( &emptyPoly );
540 : }
541 0 : GEOSGeom_destroy( point );
542 : }
543 0 : for( i = 0; i < polys->getNumGeometries(); i++ )
544 0 : GEOSGeom_destroy( ahInGeoms[i] );
545 0 : CPLFree( ahInGeoms );
546 : #endif
547 0 : poAreaReferenceLayer = 0;
548 0 : poAreaLineLayer = 0;
549 : }
|
