[gdal-dev] How can I support polygons with islands in my APPLICATION USING gdal/ogr library?

[Even Rouault]

Le jeudi 12 septembre 2013 19:56:11, sepideh a écrit :
> MyTest.rar <http://osgeo-org.1560.x6.nabble.com/file/n5077579/MyTest.rar> 
> I have uploaded my MyTest.shp created in ArcGIS. If you open it, you will
> see that it has nine polygons and in the attribute table of the layer, I
> have explained about the features of each polygon.  My problem is:  I'm
> developing a Visual C++ MFC application with GDAL/OGR library that will
> show shapefiles in an OpenGL window. And I have created this shapefile to
> test some parts of my application. *My application will only support ESRI
> shapefiles.*   As you know, the OGRPolygon class has only a member
> function *getNumInteriorRings* and not such function for exterior rings. I
> mean it supports multiple interior rings but not multiple exterior rings.
> So I decided to fetch one of the rings of each polygon with
> *getExteriorRing()* and the others with *getInteriorRing(int)* and decide
> with *isClockwise()* method whether it is an exterior or interior ring.
> This approch is best for the polygons with *FIDs : 0,1,2,3,4,5,6 * in my
> MyTest.shp. When I use this code for them:  /int NumberOfInnerRings =
> poPolygon
> ->getNumInteriorRings();/  I get these results:  FID = 0 -> 5FID = 1 ->
> 1FID = 2 -> 0FID = 3 -> 1FID = 4 -> 2FID = 5 -> 0FID = 6 -> 2  showing
> that one of the rings has been accepted as exterior and the others as
> interior polygons. You see that the polygons with the above FIDs are
> simple polygons, multipart polygons and polygons with holls inside.  But
> for the two other polygons this is not the case. These are polygons that
> have islands in as you see in this link:
> http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Creating_new_
> donut_holes_and_island_polygons/001t0000003p000000/
> <http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Creating_ne
> w_donut_holes_and_island_polygons/001t0000003p000000/> I get these wrong
> results for them as the number of Interior rings:  FID = 7 -> 1FID = 8 ->
> 2  meaning that islands are not known niether as interior nor exterior
> polygons!!!  Can you tell me how can I support islands with GDAL/OGR
> library or introduce me another open source that does this?

The shapefile displays correctly in QGIS which uses OGR to read Shapefile which 
makes me think that the problem is not in OGR.
Your above description is a bit confusing but my understanding is that your 
problem comes from the fact that the shapefile has a mix of 2 distinct geometry 
type :
- POLYGON: they have always one exterior ring, and 0 or more interior ring
- MULTIPOLYGON: that are maded of 0 or more polygons. On multipolygons, 
getNumInteriorRings() cannot be used. You must iterate first over each polygon 
part of the multipolygon

You don't need to check the clockwise orientation to determine if a ring is 
exterior or interior. This is already done by the OGR Shapefile driver.

ogrinfo with summary mode will give you the structure of the geometries :

$ ogrinfo MyTest -al -geom=summary
INFO: Open of `MyTest'
      using driver `ESRI Shapefile' successful.

Layer name: MyTest
Geometry: Polygon
Feature Count: 9
Extent: (-14519050.256578, -31065929.994238) - (26123782.221403, 
8007498.521952)
Layer SRS WKT:
PROJCS["WGS_1984_UTM_Zone_39N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",500000.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",51.0],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0]]
Id: Integer (6.0)
type_poly: String (50.0)
OGRFeature(MyTest):0
  Id (Integer) = 0
  type_poly (String) = polygonWithFiveHolesInside
  POLYGON : 22 points, 5 inner rings (11 points, 11 points, 17 points, 9 
points, 6 points)

OGRFeature(MyTest):1
  Id (Integer) = 0
  type_poly (String) = polygonWithTwoSeperateParts
  MULTIPOLYGON : 2 geometries:
POLYGON : 11 points
POLYGON : 11 points

OGRFeature(MyTest):2
  Id (Integer) = 0
  type_poly (String) = singlePolygon
  POLYGON : 10 points

OGRFeature(MyTest):3
  Id (Integer) = 0
  type_poly (String) = polygonWithOneHoleInside
  POLYGON : 11 points, 1 inner rings (17 points)

OGRFeature(MyTest):4
  Id (Integer) = 0
  type_poly (String) = polygonWithTwoHolesInside
  POLYGON : 9 points, 2 inner rings (8 points, 6 points)

OGRFeature(MyTest):5
  Id (Integer) = 0
  type_poly (String) = singlePolygon
  POLYGON : 6 points

OGRFeature(MyTest):6
  Id (Integer) = 0
  type_poly (String) = polygonWithThreeSeperateParts
  MULTIPOLYGON : 3 geometries:
POLYGON : 7 points
POLYGON : 9 points
POLYGON : 11 points

OGRFeature(MyTest):7
  Id (Integer) = 0
  type_poly (String) = PolygonWithAnIslandInside
  MULTIPOLYGON : 2 geometries:
POLYGON : 9 points, 1 inner rings (9 points)
POLYGON : 6 points

OGRFeature(MyTest):8
  Id (Integer) = 0
  type_poly (String) = PolygonWithTwoIslandsInsideOneOfTheIslandsHasAHole
  MULTIPOLYGON : 3 geometries:
POLYGON : 7 points, 1 inner rings (6 points)
POLYGON : 6 points, 1 inner rings (4 points)
POLYGON : 6 points

On the later,
- the "main" polygon is "POLYGON : 7 points, 1 inner rings (6 points)"
- the island with a hole is the "POLYGON : 6 points, 1 inner rings (4 points)" 
- the other island is "POLYGON : 6 points"


-- 
Geospatial professional services
http://even.rouault.free.fr/services.html