使用matplotlib路径计算等高线面积答案

【问题标题】：Contour area calculation using matplotlib path使用matplotlib路径计算等高线面积
【发布时间】：2018-02-06 03:17:14
【问题描述】：

我正在尝试使用 matplotlib Path 和 shapely Polygon 函数计算特定轮廓线覆盖的区域，这样做的问题是我的多边形是不规则的并且其中有孔，并且 shapely Polygpon 函数不知道哪些顶点是外部的，哪些是内部的，因此哪些应该被视为多边形的一部分，哪些应该被视为间隙的一部分。有没有办法解决这个问题？

import numpy as np
import numpy.ma as ma
from shapely.geometry import Polygon
from mpl_toolkits.basemap import Basemap

m = Basemap(projection='cyl',llcrnrlat=-30,urcrnrlat=30,llcrnrlon=30%360,\
            urcrnrlon=-130%360,lat_ts=20,resolution='c');
lons, lats = np.meshgrid(lon_sst, lat_sst);
x, y = m(lons, lats);
levels=[28.5,np.amax(data)]
cs = m.contourf(x,y,data,levels)
cn = cs.collections[0].get_paths()
vs=[]
for i in range(len(cn)):
    vs.append(cn[i].vertices)
area_of_individual_polygons = []
for i in range(len(vs)):
    area_of_individual_polygons.append(Polygon(vs[i]).area)
total_area = np.sum(area_of_individual_polygons)

【问题讨论】：

能否请您也发布您的shapely 代码？
嗨，Thomas，我的匀称代码包含在 Polygon.area 函数中

标签： python-3.x matplotlib-basemap shapely vertices contourf

【解决方案1】：

你的代码有两个问题

您处理的顶点不是您认为的那样
您的等高线图显示了“洞”，您需要从总面积中移除而不是添加，以使绿色阴影区域正确

但是让我们从头开始。 Axes.contourf() 返回一个QuadContourSet（您将其存储在cs 中。这包含matplotlib 绘制填充轮廓所需的所有数据。由于您只有两个级别，cs.contours 仅包含一项，即PathCollection。使用cs.collections[0].get_paths()，您将获得描述多边形的路径列表。但是，如果您仔细查看第一张图像，则会看到较大的绿色形状（我们称其为陆地），并带有一些较小的白色“洞” '在其中（我们称它们为湖泊，也可能是山谷）和主要陆地上的一些较小的绿色斑块（我们称它们为岛屿）。现在的问题是陆地和湖泊都是一条路，而岛屿是不同的路径。

首先定义一个方向的轮廓（比如顺时针），然后定义相反方向的湖泊轮廓（然后是逆时针方向）来定义陆地的轮廓。为了让 matplotlib 知道如何处理这个问题，使用 codes 将路径分成几部分。这里使用了三个代码：Path.MOVETO (=1)、Path.LINETO (=2) 和 Path.CLOSEPOLY (=79)。您可以使用path.codes 获得这些代码的列表（见下文）。每个段都以Path.MOVETO 开头，以Path.CLOSEPOLY 结尾，中间有很多PATH.LINETOs（至少在这种情况下，没有弯曲路径）。每个代码都有一个相应的顶点（您使用path.vertices 获得，只有与Path.CLOSEPOLY 对应的顶点被忽略（通常设置为（0,0）），因此您必须将它们删除以进行最终计算。在帮助下在这些代码中，每条路径都被分割成段。

告诉你这一切之后，你需要做的不仅是从每条路径中检索顶点，还要检索相应的代码。然后你需要将代码和顶点分割成段（使用代码这很容易），然后分别计算每个段的面积。最重要的是，要正确计算阴影区域的面积，您必须将各个 first 路径段的面积相加并减去其他段的面积。这是一个如何做到这一点的示例：

from matplotlib import pyplot as plt
import numpy as np
import numpy.ma as ma
from shapely.geometry import Polygon
from mpl_toolkits.basemap import Basemap

from matplotlib.path import Path
from matplotlib.patches import PathPatch

##setup
lonmin, lonmax = 30%360,-130%360
latmin, latmax = -30, 30
lon_sst = np.linspace(lonmin, lonmax, 50)
lat_sst = np.linspace(latmin, latmax, 50)
lons, lats = np.meshgrid(lon_sst, lat_sst);

##generating some example data
lonmid = 0.5*(lonmin+lonmax)
latmid = 0.5*(latmin+latmax)
data = (
    10*np.cos(np.deg2rad(lons-lonmid))**2*np.cos(3*np.deg2rad(lats-latmid))**2
    -10*np.exp(-0.01*((lons-(lonmin+2*lonmid)/3)**2+2*(lats-latmid)**2))
    -10*np.exp(-0.01*((lons-(lonmax+2*lonmid)/3)**2+(lats-latmid)**2))
    +10*np.exp(-0.1*((lons-(lonmid+2*lonmax)/3)**2+2*(lats-(latmid+2*latmax)/3)**2))
    )

##opening figure and axes
fig,ax = plt.subplots()

##do the Basemap stuff
m = Basemap(
    ax=ax,projection='cyl',
    llcrnrlat=-30,    urcrnrlat=30,
    llcrnrlon=30%360, urcrnrlon=-130%360,
    lat_ts=20,resolution='c'
)
x, y = m(lons, lats);
levels=[2,np.amax(data)]
cs = ax.contourf(x,y,data,levels)

##organizing paths and computing individual areas
paths = cs.collections[0].get_paths()
help(paths[0])
area_of_individual_polygons = []
for p in paths:
    sign = 1  ##<-- assures that area of first(outer) polygon will be summed
    verts = p.vertices
    codes = p.codes
    idx = np.where(codes==Path.MOVETO)[0]
    vert_segs = np.split(verts,idx)[1:]
    code_segs = np.split(codes,idx)[1:]
    for code, vert in zip(code_segs,vert_segs):

        ##visualising (not necessary for the calculation)
        new_path = Path(vert,code)
        patch = PathPatch(
            new_path,
            edgecolor = 'black' if sign == 1 else 'red',
            facecolor = 'none',
            lw =1
        )
        ax.add_patch(patch)

        ##computing the area of the polygon
        area_of_individual_polygons.append(sign*Polygon(vert[:-1]).area)
        sign = -1 ##<-- assures that the other (inner) polygons will be subtracted

##computing total area        
total_area = np.sum(area_of_individual_polygons)
formstring = ''.join(['{:+.2}' for a in area_of_individual_polygons])+'={:.2}'
print(formstring.format(*area_of_individual_polygons,total_area))

plt.show()

这是与上述代码对应的图像。请注意，我在 black 中列出了添加的区域以及在 red 中减去的区域。另请注意，仍有一些极端情况会出错（例如湖中的岛屿）。

这里仍然是print命令的输出，只是代码中发生的计算的可视化：

+4.5e+03-3.2e+02-1.9e+02+4.3e+01=4e+03

希望这可以解决问题。

【讨论】：