Index: ogr/ogrsf_frmts/shape/makefile.vc =================================================================== --- ogr/ogrsf_frmts/shape/makefile.vc (révision 12190) +++ ogr/ogrsf_frmts/shape/makefile.vc (copie de travail) @@ -1,7 +1,7 @@ OBJ = shape2ogr.obj shpopen.obj dbfopen.obj ogrshapedriver.obj \ ogrshapedatasource.obj ogrshapelayer.obj shptree.obj -EXTRAFLAGS = -I.. -I..\.. /DSHAPELIB_DLLEXPORT -DUSE_CPL +EXTRAFLAGS = -I.. -I..\.. /DSHAPELIB_DLLEXPORT -DUSE_CPL $(GEOS_CFLAGS) GDAL_ROOT = ..\..\.. Index: ogr/ogrsf_frmts/shape/GNUmakefile =================================================================== --- ogr/ogrsf_frmts/shape/GNUmakefile (révision 12190) +++ ogr/ogrsf_frmts/shape/GNUmakefile (copie de travail) @@ -5,6 +5,10 @@ OBJ = shape2ogr.o shpopen.o dbfopen.o shptree.o \ ogrshapedriver.o ogrshapedatasource.o ogrshapelayer.o +ifeq ($(HAVE_GEOS),yes) +CPPFLAGS := -DHAVE_GEOS=1 $(GEOS_CFLAGS) $(CPPFLAGS) +endif + CPPFLAGS := -DUSE_CPL -I.. -I../.. $(GDAL_INCLUDE) $(CPPFLAGS) default: $(O_OBJ) Index: ogr/ogrsf_frmts/shape/shape2ogr.cpp =================================================================== --- ogr/ogrsf_frmts/shape/shape2ogr.cpp (révision 12190) +++ ogr/ogrsf_frmts/shape/shape2ogr.cpp (copie de travail) @@ -306,8 +306,247 @@ return ( poRing ); } - + + /************************************************************************/ +/* MakeMultiPolygon() */ +/************************************************************************/ + + +enum +{ + CONTAINS, + IS_CONTAINED_BY, + NOT_RELATED +}; + +static OGRGeometry* MakeMultiPolygon (OGRPolygon const * const * polygons, + int nbPolygons, + int* pIsValidGeometry) +{ + int i, j; + + if (nbPolygons == 1) + { + return polygons[0]->clone(); + } + else + { +#ifndef HAVE_GEOS + static int firstTime = 1; + if (firstTime) + { + CPLError(CE_Warning, CPLE_AppDefined, + "GDAL should be built with GEOS support enabled in order the " + "SHAPE driver to provide reliable results on complex polygons."); + firstTime = 0; + } +#endif + + OGREnvelope* envelopes = new OGREnvelope[nbPolygons]; + int** relations = new int*[nbPolygons]; + double* areas = new double[nbPolygons]; + int go_on = 1; + for(i=0;igetEnvelope(&envelopes[i]); + areas[i] = polygons[i]->get_Area(); + } + + /* This a several step algorithm : + 1) Compute in a matrix how polygons relate to each other + (this is the moment for detecting pathological intersections and exiting) + 2) For each polygon, find the smallest enclosing polygon + 3) For each polygon, compute its inclusion depth (0 means toplevel) + 4) For each polygon of odd depth (= inner ring), add it to its outer ring + + Complexity : O(nbPolygons^2) + */ + + /* Compute how each polygon relate to the other ones + To save a bit of computation we always begin the computation by a test on the + enveloppe. We also take into account the areas to avoid some useless tests. + (A contains B implies envelop(A) contains envelop(B) and area(A) > area(B)) + In practise, we can hope that few full geometry intersection of inclusion test is done: + * if the polygons are well separated geographically + (a set of islands for example), no full geometry intersection or inclusion + test is done. (the envelopes don't intersect each other) + * if the polygons are 'lake inside an island inside a lake inside an area' and + that each polygon is much smaller than its enclosing one, their bounding boxes + are stricly contained into each oter, and thus, no full geometry intersection or inclusion + test is done + */ + for(i=0;go_on && i areas[j] + && envelopes[i].Contains(envelopes[j]) +#ifdef HAVE_GEOS + && polygons[i]->Contains(polygons[j]) +#endif + ) + { + relations[i][j] = CONTAINS; + relations[j][i] = IS_CONTAINED_BY; + } + else if (areas[j] > areas[i] + && envelopes[j].Contains(envelopes[i]) +#ifdef HAVE_GEOS + && polygons[j]->Contains(polygons[i]) +#endif + ) + { + relations[j][i] = CONTAINS; + relations[i][j] = IS_CONTAINED_BY; + } + else if (envelopes[i].Intersects(envelopes[j]) == FALSE +#ifdef HAVE_GEOS + /* We use Overlaps instead of Intersects to be more tolerant about touching polygons */ + || polygons[i]->Overlaps(polygons[j]) == FALSE +#endif + ) + { + relations[i][j] = NOT_RELATED; + relations[j][i] = NOT_RELATED; + } + else + { + /* Bad... The polygons are intersecting but no one is + contained inside the other one. This is a really broken + case. We just make a multipolygon with the whole set of + polygons */ + go_on = 0; +#if 0 + fprintf(stderr, "bad intersection for polygons %d and %d\n", i, j); + char* wkt1; + char* wkt2; + polygons[i]->exportToWkt(&wkt1); + polygons[j]->exportToWkt(&wkt2); + fprintf(stderr, "geom %d : %s\n", i, wkt1); + fprintf(stderr, "geom %d : %s\n", j, wkt2); + CPLFree(wkt1); + CPLFree(wkt2); +#endif + } + } + } + + if (pIsValidGeometry) + *pIsValidGeometry = go_on; + + OGRGeometry* geom = NULL; + + if (go_on == 0) + { + geom = new OGRMultiPolygon(); + for(i=0;iaddGeometry(polygons[i]); + } + } + else + { + int* directContainerIndex = new int[nbPolygons]; + + /* Find the smallest enclosing polygon of each polygon */ + for(i=0;i= 0) + { + j = directContainerIndex[j]; + depth++; + } + containedDepth[i] = depth; +// fprintf(stderr, "%d is of depth %d\n", i, depth); + } + + int nbTopLevelPolygons = 0; + OGRPolygon** tempPolygons = new OGRPolygon*[nbPolygons]; + + /* Create a copy of toplevel polygons */ + for(i=0;iclone(); + if (nbTopLevelPolygons == 1) + geom = tempPolygons[i]; + } + } + + /* Add interior rings to toplevel polygons */ + for(i=0;iaddRing((OGRLinearRing*)polygons[i]->getExteriorRing()); + } + } + + if (nbTopLevelPolygons > 1) + { + geom = new OGRMultiPolygon(); + + /* Add toplevel polygons to the multipolygon */ + for(i=0;iaddGeometryDirectly(tempPolygons[i]); + } + } + } + + delete[] tempPolygons; + delete[] directContainerIndex; + delete[] containedDepth; + } + + for(i=0;inParts]; + for( iRing = 0; iRing < psShape->nParts; iRing++ ) + { + tabPolygons[iRing] = new OGRPolygon(); + tabPolygons[iRing]->addRingDirectly(CreateLinearRing ( psShape, iRing )); + } + + int isValidGeometry; + poOGR = MakeMultiPolygon (tabPolygons, + psShape->nParts, + &isValidGeometry); + + if (!isValidGeometry) + { + CPLError(CE_Warning, CPLE_AppDefined, + "Geometry of polygon of fid %d cannot be translated to Simple Geometry. " + "All polygons will be contained in a multipolygon.\n", + iShape); + } + + for( iRing = 0; iRing < psShape->nParts; iRing++ ) + { + delete tabPolygons[iRing]; + } + delete[] tabPolygons; +#else /* Multipart polygon. */ int nOuter = 0; /* Number of outer rings */ @@ -637,7 +903,7 @@ CPLFree( direction ); CPLFree( outer ); CPLFree( outside ); - +#endif } /* End of multipart polygon processing. */ if( psShape->nSHPType == SHPT_POLYGON )