[gdal-dev] Shapefile - Multipolygon with interior rings

[Jay L.]

Good afternoon,
I am working on a program to convert a shapefile from ocentric to ographic
latitude.  The program reads an ESRI shapefile, creates a new shapefile,
copying the attribute data, and performs a latitude "shift" for each point
in the geometry.  OGR is working wonderfully on all geometry types, except
for multipolygons with donuts.

Below is the code snippet.  Above this the geometry type and feature
attributes are read from the input shapefile.

When printing the point_count (in blue) variable I am able to see that my
multipolygon consists of two polygons, each with 2 rings (one interior and
one exterior).  I iterate through these rings and attempt to get a point
count, which I use to iterate over each linearring.  The point count for the
exterior ring runs without a problem and returns the proper number of
points.  The point count on the interior ring returns an error.

Stepping up one "level" I can print the geometry (green text) for the
multipolygon and the description of 2 polygons (in proper WKT format with an
exterior and interior ring) are written.  It seems that once I am "into" the
linearrings of a multipolygon OGR is only "seeing" the exterior ring.

Thoughts?

Thanks!
Jay


Code Snippet:

elif geom_type == "MULTIPOLYGON":

    while feature_in is not None:

        #Create the output feature from the input
        geometry = feature_in.GetGeometryRef()

        #Dig into the input geometry
        geom_count = geometry.GetGeometryCount()

        #Pack the multipolygon lunch for later
        multipoly = ogr.Geometry(ogr.wkbMultiPolygon)

        #Iterate through the polygons
        for g in range(geom_count):
            geom = geometry.GetGeometryRef(g)
            geom_count = geom.GetGeometryCount()
            geom_name = geom.GetGeometryName()

            #For individual polygons in the multipolygon layer
            if geom_name == "LINEARRING":

                #Create a linear ring to populate with the points.
                ring = ogr.Geometry(ogr.wkbLinearRing)

                #Get a point count
                point_count = geom.GetPointCount()

                #Iterate through the points
                for h in range (point_count):
                    pntx = geom.GetX(h)
                    pnty = geom.GetY(h)
                    pntz = geom.GetZ(h)

                    #perform the conversion
                    if pntz == 0.0:
                        new_y = convert(pnty)
                        ring.AddPoint_2D(pntx, new_y)


                    elif pntz != 0.0:
                        new_y = convert(pnty)
                        ring.AddPoint(pntx, new_y, pntz)

                #Add the linearring
                poly.AddGeometry(ring)

                #Assign the spatial reference to the feature
                poly.AssignSpatialReference(feature_in_spatref)

            #For multipolygons in the multipolygon layer
            elif geom_name == "POLYGON":

                #Create a polygon to populate with linear ring(s)
                poly = ogr.Geometry(ogr.wkbPolygon)

                #Iterate through the polygons
                for h in range (geom_count):
                    geom = geom.GetGeometryRef(h)

                    #Check for multipolyons with donuts
                    try:
                        #Get a point count
                        point_count = geom.GetPointCount()

                        #Create a linear ring to populate with the points.
                        ring = ogr.Geometry(ogr.wkbLinearRing)

                        #Iterate through the points
                        for j in range (point_count):
                            pntx = geom.GetX(j)
                            pnty = geom.GetY(j)
                            pntz = geom.GetZ(j)

                            #perform the conversion
                            if pntz == 0.0:
                                new_y = convert(pnty)
                                ring.AddPoint_2D(pntx, new_y)


                            elif pntz != 0.0:
                                new_y = convert(pnty)
                                ring.AddPoint(pntx, new_y, pntz)

                    #Exit gracefully if multipolygons with donuts are
encountered
                    except:
                        print "This shapefile is not supported."

                #Add the linearring
                poly.AddGeometry(ring)

            #Add the polygon to the multipolygon
            multipoly.AddGeometry(poly)

            #Set the spatial reference of the feature
            multipoly.AssignSpatialReference(feature_in_spatref)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100809/827a46e9/attachment-0001.html