Ticket #1217: gdal_svn_trunk_ogr_shape_bug1217.patch

ogr/ogrsf_frmts/shape/makefile.vc

@@ -1 +1 @@
-
-OBJ     =       shape2ogr.obj shpopen.obj dbfopen.obj ogrshapedriver.obj \
-                ogrshapedatasource.obj ogrshapelayer.obj shptree.obj
-
+ $(GEOS_CFLAGS)
-
-GDAL_ROOT       =       ..\..\..
-

ogr/ogrsf_frmts/shape/GNUmakefile

@@ -5 +5 @@
-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)

ogr/ogrsf_frmts/shape/shape2ogr.cpp

@@ -306 +306 @@
-
-    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;i<nbPolygons;i++)
+        {
+            relations[i] = new int[nbPolygons];
+            polygons[i]->getEnvelope(&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<nbPolygons;i++)
+        {
+            for(j=i+1;go_on && j<nbPolygons;j++)
+            {
+                if (areas[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;i<nbPolygons;i++)
+            {
+                ((OGRMultiPolygon*)geom)->addGeometry(polygons[i]);
+            }
+        }
+        else
+        {
+            int* directContainerIndex = new int[nbPolygons];
+
+            /* Find the smallest enclosing polygon of each polygon */
+            for(i=0;i<nbPolygons;i++)
+            {
+                int jSmallestContainer = -1;
+                double areaSmallestContainer = 0;
+                for(j=0;j<nbPolygons;j++)
+                {
+                    if (i != j)
+                    {
+                        if (relations[i][j] == IS_CONTAINED_BY)
+                        {
+                            if (jSmallestContainer < 0 || areas[j] < areaSmallestContainer)
+                            {
+                                jSmallestContainer = j;
+                                areaSmallestContainer = areas[j];
+                            }
+                        }
+                    }
+                }
+                directContainerIndex[i] = jSmallestContainer;
+            }
+
+            /* Compute the inclusion depth of each polygon */
+            int* containedDepth = new int [nbPolygons];
+            for(i=0;i<nbPolygons;i++)
+            {
+                int depth = 0;
+                int j = directContainerIndex[i];
+                while (j >= 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;i<nbPolygons;i++)
+            {
+                if ((containedDepth[i] % 2) == 0)
+                {
+                    nbTopLevelPolygons ++;
+                    tempPolygons[i] = (OGRPolygon*)polygons[i]->clone();
+                    if (nbTopLevelPolygons == 1)
+                        geom = tempPolygons[i];
+                }
+            }
+
+            /* Add interior rings to toplevel polygons */
+            for(i=0;i<nbPolygons;i++)
+            {
+                if ((containedDepth[i] % 2) == 1)
+                {
+                    tempPolygons[directContainerIndex[i]]->addRing((OGRLinearRing*)polygons[i]->getExteriorRing());
+                }
+            }
+
+            if (nbTopLevelPolygons > 1)
+            {
+                geom = new OGRMultiPolygon();
+
+                /* Add toplevel polygons to the multipolygon */
+                for(i=0;i<nbPolygons;i++)
+                {
+                    if ((containedDepth[i] % 2) == 0)
+                    {
+                        ((OGRMultiPolygon*)geom)->addGeometryDirectly(tempPolygons[i]);
+                    }
+                }
+            }
+
+            delete[] tempPolygons;
+            delete[] directContainerIndex;
+            delete[] containedDepth;
+        }
+
+        for(i=0;i<nbPolygons;i++)
+        {
+            delete[] relations[i];
+        }
+        delete[] relations;
+        delete[] areas;
+        delete[] envelopes;
+
+        return geom;
+    }
+}
+
+
+
+
+
+/************************************************************************/
-/*                          SHPReadOGRObject()                          */
-/*                                                                      */
-/*      Read an item in a shapefile, and translate to OGR geometry      */
@@ -460 +699 @@
-        }
-        else
-        {
+#ifdef HAVE_GEOS
+            OGRPolygon** tabPolygons = new OGRPolygon*[psShape->nParts];
+            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 +903 @@
-            CPLFree( direction );
-            CPLFree( outer );
-            CPLFree( outside );
-
+
-        } /* End of multipart polygon processing. */
-
-        if( psShape->nSHPType == SHPT_POLYGON )