Pyshp draws polygons out of order after intersection of two shapefiles

Question

I have two shapefiles. One contains electoral borders but not coastline and the other contains provinces and coastlines. I am plotting the electoral districts using Python's pyshp and filling in the colors according to some data. That all works on the original shapefile:

And here is the important part of my plotting script:

sf = shp.Reader(shapefile)
records = []
shapes = []
for i,j in zip(sf.records(),sf.shapeRecords()):
    records.append(i)
    shapes.append(j)

plt.figure()
for rec,shape in zip(records,shapes):
    xy = [[i[0],i[1]] for i in shape.shape.points]
    plt.fill([i[0] for i in xy], [i[1] for i in xy], lw=0.2, edgecolor='white', facecolor=cmap(<some numbers from a list>))

I wanted the electoral districts to have nicer looking borders on the coasts so I used the intersect feature in QGIS to create a new shapefile of the electoral district file as an input layer and the coastline file as the intersecting layer. The result of that displays just fine in QGIS:

Here are the original shapefiles that I intersected shown overlapping in QGIS:

As you can hopefully see from the screen grab, some electoral districts get broken into multiple polygons because they contain islands, for example, but maintain their attributes (which is what I want).

However, when I use the same Python script above it seems that the plotting order of the points is all messed up resulting in oddly shaped polygons. It looks like the drawing is jumping across points previously part of the same electoral district and the same polygon but now apart of different polygons but still the same district.

Can anyone help? I can provide more details upon request.

EDIT:

I've narrowed the problem down to the way I'm handling multipart polygons. The intersection didn't break up polygons into new polygons, it just created multiple parts for each one and I'm not sure how to use pyshp to draw the parts separately.

Answer 1

Okay, I figured it out after my last edit. It was indeed a problem with polygons containing multiple parts. Here is my fixed figure.

So my first attempt at a solution was to apply a Multipart-to-Single-Polygon function in QGIS. That improved things but did not eliminate the fractured appearance entirely. A quick look at sf.shapeRecords[i].shape.parts showed that there were still multipart shapes in my shapefile. Fortunately the shape.parts attribute is a list containing indices of where new parts begin in the shape.points list. So you simply need to break up your shape drawing at the appropriate index and restart drawing the next part until the next index. If the shape.parts attribute contains the list [0] then there is only one part. Here's a minimal example of my implementation:

for shape in sf.shapeRecords:
    if len(shape.shape.parts) == 1:
        x = [i[0] for i in shape.shape.points]             
        y = [i[1] for i in shape.shape.points]
        plt.fill(x, y)
    else:
        for index in range(len(shape.shape.parts)):
            try:
                x = [i[0] for i in shape.shape.points[shape.shape.parts[index]:shape.shape.parts[index+1]]]
                y = [i[1] for i in shape.shape.points[shape.shape.parts[index]:shape.shape.parts[index+1]]]
            except:
                x = [i[0] for i in shape.shape.points[shape.shape.parts[index]:]]
                y = [i[1] for i in shape.shape.points[shape.shape.parts[index]:]]
            plt.fill(x, y)

In words, you iterate through your shapes and check each one to see if its parts list has only one entry. If it has one entry it means there is only one part so you can plot as normal. If it has more than one entry then you break the points list into sections according to the indices given in the parts list and plot each one separately.

Hopefully this saves someone some trouble later on.