2

I'm trying to draw a map of some regions in a shapefile in Python. My basic approach is this:

shp = fiona.open("C:/Users/nils/Documents/Maps/my_shapefile.shp")
bds = shp.bounds
ll = (bds[0], bds[1])
ur = (bds[2], bds[3])
coords = list(ll + ur)
w, h = coords[2] - coords[0], coords[3] - coords[1]

# Make figure instance, add Basemap and CCG boundaries from shapefile
fig, ax = plt.subplots(figsize=(12,10))
m = Basemap(projection="tmerc", lon_0 = -2., lat_0 = 49., ellps="WGS84",
        llcrnrlon = coords[0], llcrnrlat = coords[1],
        urcrnrlon = coords[2], urcrnrlat = coords[3],
        lat_ts = 0, resolution="i", suppress_ticks=True)
m.readshapefile("C:/Users/nils/Documents/Maps/my_shapefile.shp", "Regions")

# Extract polygon coordinates of  and names of regions to plot from shapefile
to_plot = ["region_A", "region_B", "region_C"]
poly = []; name = []
for coordinates, region in zip(m.Regions, m.Regions_info):
    if any(substr in region["name"] for substr in to_plot):
        poly.append(Polygon(coordinates))
        name.append(region["name"])

# Turn polygons into patches using descartes
patches = []
for i in poly:
    patches.append(PolygonPatch(i, facecolor='#006400', edgecolor='#787878', lw=0.25, alpha=0.5))

# Add PatchCollection to basemap
ax.add_collection(PatchCollection(patches, match_original=True))

Now my problem with this is that the shapefile covers a larger geographical area, but I only want to plot a subset of this area (think e.g. I have a UK shapefile, but want to plot a map of all regions in Wales). Now I can identify the correct regions and only add those patches as in the example above, but matplotlib will still plot the boundaries of all regions in the shapefile, and the boundaries identified by fiona's bounds method are obviously independent of the subset of patches I've chosen.

I have two questions relating to this:

  1. How can I get matplotlib to only draw the boundaries of a subset of patches defined in the shapefile?

  2. How can I get the bounds of a subset of patches, similar to what fiona's bound method does to the entire shapefile?

Nils Gudat
  • 13,222
  • 3
  • 39
  • 60

2 Answers2

2

To answer the second part as well, here's a function that achieves the desired result:

def get_bounds(patch_list):

m = Basemap()
# Read in shapefile, without drawing anything
m.readshapefile("C:/Users/ngudat/Documents/Maps/CCG/CCG_boundaries_2015", "patches", drawbounds=False)

# initialize boundaries (this is a bit of a manual step, might be better to intialize to boundaries of first patch or something)
lon_min = 0.
lon_max = -3.
lat_min = 60.
lat_max = 0.

for (shape, patch_name) in zip(m.patches, m.patches_info):
    if patches["name"] in patch_list:
        lon, lat = zip(*shape)
        if min(lon) < lon_min:
            lon_min = min(lon)
        if max(lon) > lon_max:
            lon_max = max(lon)
        if min(lat) < lat_min:
            lat_min = min(lat)
        if max(lat) > lat_max:
            lat_max = max(lat)

return lon_min, lat_min, lon_max, lat_max

Probably not the most efficient way of doing it, and the initialization step might need to be changed, but the idea should be applicable to similar situations easily.

Nils Gudat
  • 13,222
  • 3
  • 39
  • 60
0

To (partly) answer my first question, a way of achieving this would be to call readshapefile with drawbounds=False; in this case there won't be any boundaries on the map, but the boundaries of my selection of regions will be drawn with the patches in any case.

Nils Gudat
  • 13,222
  • 3
  • 39
  • 60