Write Your Own Geospatial Utilities - Open Source Geospatial Tools: Applications in Earth Observation - page 219

#using the GetExtent() method, we extract the layer's
→ spatial extent (bounding box)
layer_extent = layer.GetExtent()
xoff = (layer_ext[1]-layer_ext[0])/50
yoff = (layer_ext[3]-layer_ext[2])/50

# Prepare figure
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlim(ext[0]-xoff,ext[1]+xoff)
ax.set_ylim(ext[2]-yoff,ext[3]+yoff)

We initialize the paths as a list variables and reset the layer using ResetReading().

paths = []
lyr.ResetReading()

The next code block forms the crux of the program. It begins by looping through all
of the features stored within the layer and contains a nested for loop that iterates over
the geometries of the individual features. In this instance, we will write the program
to handle polygon features and so we check that the correct OGR GeometryType has
been passed, namely ogr.wkbPolygon .

# Read all features in layer and store as paths
for feat in lyr:
for geom in feat.GetGeometryRef():
# check if geom is polygon
if geom.GetGeometryType() == ogr.wkbPolygon:
codes = []
all_x = []
all_y = []
for i in range(geom.GetGeometryCount()):
# Read ring geometry and create path
r = geom.GetGeometryRef(i)
x = [r.GetX(j) for j in range(r.GetPointCount())

→ ]

y = [r.GetY(j) for j in range(r.GetPointCount())

→ ]

# skip boundary between individual rings
codes += [mpath.Path.MOVETO] + \
(len(x)-1)*[mpath.Path.LINETO]

all_x += x
all_y += y
path = mpath.Path(np.column_stack((all_x,all_y)),

→ codes)

paths.append(path)

# Add paths as patches to axes
for path in paths:
patch = mpatches.PathPatch(path, \
facecolor=color, edgecolor='black')
ax.add_patch(patch)