2

我正在尝试在 Python 的 shapefile 中绘制某些区域的地图。我的基本方法是这样的：

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))

现在我的问题是 shapefile 覆盖了更大的地理区域，但我只想绘制该区域的一个子集（例如，我有一个英国 shapefile，但想绘制威尔士所有地区的地图）。现在我可以识别正确的区域并只添加上面示例中的那些补丁，但是 matplotlib 仍然会绘制 shapefile 中所有区域的边界，并且由 fiona 的bounds方法识别的边界显然与我的补丁子集无关选择。

我有两个与此相关的问题：

  1. 如何让 matplotlib 仅绘制 shapefile 中定义的补丁子集的边界？

  2. 如何获得补丁子集的边界，类似于 fiona 的bound方法对整个 shapefile 所做的那样？

4

2 回答 2

2

为了回答第二部分，这里有一个可以达到预期结果的函数：

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

可能不是最有效的方法，并且可能需要更改初始化步骤，但该想法应该很容易适用于类似情况。

于 2016-02-17T10:13:36.870 回答
0

为了（部分）回答我的第一个问题，实现此目的的一种方法是readshapefile使用drawbounds=False; 在这种情况下，地图上不会有任何边界，但无论如何，我选择的区域的边界都会用补丁绘制。

于 2016-01-07T14:25:08.970 回答