在底图中部分填充 shapefile 多边形

Question

我想使用 basemap/shapefile 数据填充多边形，但只能填充一定的百分比。例如，在下面的示例中，我们根据值进行填充，但假设我想根据这些值填充 % 的多边形（来自的代码）：

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np

fig= plt.figure()
ax= fig.add_subplot(111)
m=Basemap(projection='cyl',llcrnrlat=34.5,llcrnrlon=19,
                           urcrnrlat=42,urcrnrlon=28.5,resolution='h')
m.drawmapboundary(fill_color='aqua')
m.fillcontinents(color='w',lake_color='aqua')
m.drawcoastlines()
m.readshapefile('data/nomoi/nomoi','nomoi')

dict1={14464: 1.16, 14465: 1.35, 14466: 1.28, 14467: 1.69, 14468: 1.81, 14418: 1.38}
colvals = dict1.values()

cmap=plt.cm.RdYlBu
norm=plt.Normalize(min(colvals),max(colvals))

patches   = []

for info, shape in zip(m.nomoi_info, m.nomoi):
    if info['ID_2'] in list(dict1.keys()):
        color=cmap(norm(dict1[info['ID_2']]))
        patches.append( Polygon(np.array(shape), True, color=color) )

pc = PatchCollection(patches, match_original=True, edgecolor='k', linewidths=1., zorder=2)
ax.add_collection(pc)

#colorbar
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array(colvals)
fig.colorbar(sm, ax=ax)

plt.show()

谢谢。

Answer 1

import math
from shapely.geometry import Polygon as shpoly

#shapefile of main massachusetts shape
iowpoly = state_shapes['Massachusetts'][32]

def return_xy(coords):
    
    return [np.asarray([i[0] for i in coords]), np.asarray([i[1] for i in coords])]

def return_area(coords):
    
    x, y = return_xy(coords)
    
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

def return_bounding_box(coords):
    x, y = return_xy(coords)
    
    return [[min(x), min(y)], [max(x), max(y)]]

def split_x_wise(bbox, weights, split = 2):
    
    lleft = bbox[0]
    uright = bbox[1]
    
    dx = abs(uright[0] - lleft[0])
    
    weights = np.cumsum(sorted(weights, reverse=True))
    
    xcoords = [lleft[0]+weights[x-1]*dx for x in range(1, split)]

    return xcoords

def generate_splits_by_area(coords, bbox, weights, tolerance = 0.03, div = 100):
    
    xareasplits = {}
    
    weights = np.cumsum(sorted(weights, reverse=True))[:-1]
    
    lleft = bbox[0]
    uright = bbox[1]
    dx = abs(uright[0] - lleft[0])

    xsplit = [lleft[0]+(dx/div)*x for x in range(1, div)]
    
    for w in weights:
        xareasplits[str(w)] = None
    
    mainarea = shpoly(coords).area
    
    for i, s in enumerate(xsplit):
    
        poly = []
    
        if i == 0:
            continue
        
        for ip, p in enumerate(coords):
            if p[0] < s:
                poly.append(p)
                
        shpl = shpoly(poly).area

        frac = shpl/mainarea

        for w in weights:
            if abs(w-frac) <= tolerance:
                if xareasplits[str(w)] == None:
                    xareasplits[str(w)] = s
                    
    return list(xareasplits.values())
        
def return_split(coords, weights, split = 2, by_area = False, tolerance = 0.03, div = 100):

    polys = {}
    
    for x in range(0, split):
        polys[str(x+1)] = {'points':[], 'maxit' : None}
        
    bbox = return_bounding_box(coords)
        
    if not by_area:
        xsplit = split_x_wise(bbox, weights, split)
        #test = generate_splits_by_area(coords, bbox, weights, tolerance=tolerance, div=div)
    else:
        xsplit = generate_splits_by_area(coords, bbox, weights, tolerance=tolerance, div=div)
    
    xsplit.append(bbox[0][0])
    xsplit.append(bbox[1][0])
    
    xsplit = sorted(xsplit)
    
    #print(xsplit)
    #print(test)
    
    for ip, p in enumerate(coords):
        for i, splt in enumerate(xsplit):
            if i > 0:
                if (p[0] > xsplit[i-1]) & (p[0] < splt):

                    if len(polys[str(i)]['points']) == 0:
                        polys[str(i)]['points'].append(coords[ip-1])

                    polys[str(i)]['points'].append(p)
                    polys[str(i)]['maxit'] = ip
                    
    for poly, data in polys.items():
        
        tmaxit = data['maxit']+1
        
        if tmaxit >= len(coords):
            data['points'].append(coords[0])
        else:
            data['points'].append(coords[tmaxit])
    
    return polys
    #return [p for p in coords if p[0] > xsplit[0]]

#bboxiowa = return_bounding_box(iowpoly)
splitpoly = return_split(iowpoly, weights = [0.2780539772727273, 0.1953716856060606, 0.19513494318181818, 0.18329782196969696, 0.14814157196969696],by_area = True,split = 5)

for k, v in splitpoly.items():
    print (k, len(v['points']))
    print (v['maxit'])
    
test = shpoly(splitpoly["1"]['points'])
test

我设法编写了自己的代码来从 shapefile 中拆分和填充 shape 多边形。上面的代码示例将马萨诸塞州的 shapefile 分成 5 段，根据权重和面积加权。

拆分的前两部分如下所示：

在底图中部分填充 shapefile 多边形

Filling shapefile polygons partially in basemap

python

matplotlib

shapefile

fill

matplotlib-basemap