0

试图摆脱 Basemap 以在地图上绘制轮廓，特别是从 Jupyter 笔记本上这样做，我遇到了与 folium 相关的讨论：https ://github.com/python-visualization/folium/issues/958 。嘿！我想 - 这很棒，让我用 ipyleaflet 试试。

下面的代码让我接近了一个可行的解决方案，但我看到了一些模糊的轮廓，如下图所示。

注意：读取数据的代码被跳过

import matplotlib.pyplot as plt
from ipyleaflet import Map, basemaps, Polygon

# generate matplotlib contours
cs = plt.contourf(lat, lon, val.T)  # note transposition to obtain correct coordinate order
# set-up map
zoom = 4
center = [50., 0.]
map = Map(basemap=basemaps.CartoDB.Positron, center=center, zoom=zoom)
# add contours as polygons
# get the pathlines of the contours with cs.allsegs
for clev in cs.allsegs:
    polygons = Polygon(
                locations=[p.tolist() for p in clev],
                color="green",
                weight=1,
                opacity=0.4,
                fill_color="green",
                fill_opacity=0.3
    )
    map.add_layer(polygons);
map

对于下图，我从 ECMWF 的 ERA-5 再分析中加载了数据，并绘制了 2020 年 2 月的水平风速图。仅选择了一个等高线来演示该问题。下图中的左侧面板显示了plt.contourf(lon, lat, val)结果，这是正确的。右侧的面板显示了用传单绘制的轮廓。虽然大部分等高线显示正确（请仅关注浅绿色区域，即第二高等高线级别），但线段的排序似乎存在一些问题。有谁知道如何解决这个问题？

等高线问题图解

4

1 回答 1

1

在下面的函数split_contours的帮助下，代码现在可以按预期工作。Path的文档表明轮廓属性allkinds可能为 None。这由下面的代码处理。不包括的是路径段代码 3 和 4，它们对应于贝塞尔曲线。

cs = plt.contourf(lat, lon, val.T)
plt.close()

# get the pathlines of the contours:
# print(cs.allsegs[0][0][0:12])  # returns list of polygons per contour level; each polygon is a list of [x,y]tuples
# i.e. [level][polygon][x,y]
# print(len(cs.allsegs), len(cs.allsegs[0]))

# other useful information from the contour object:
# print(cs.get_array())  # get contour levels

# plot map
from ipyleaflet import Map, basemaps, Polygon

zoom = 4
center = [50., 0.]
# map = Map(basemap=basemaps.NASAGIBS.ViirsTrueColorCR, center=center, zoom=zoom) # loads current satellite image
# map = Map(basemap=basemaps.NASAGIBS.ViirsEarthAtNight2012, center=center, zoom=zoom)
map = Map(basemap=basemaps.CartoDB.Positron, center=center, zoom=zoom)

def split_contours(segs, kinds=None):
    """takes a list of polygons and vertex kinds and separates disconnected vertices into separate lists.
    The input arrays can be derived from the allsegs and allkinds atributes of the result of a matplotlib
    contour or contourf call. They correspond to the contours of one contour level.
    
    Example:
    cs = plt.contourf(x, y, z)
    allsegs = cs.allsegs
    allkinds = cs.allkinds
    for i, segs in enumerate(allsegs):
        kinds = None if allkinds is None else allkinds[i]
        new_segs = split_contours(segs, kinds)
        # do something with new_segs
        
    More information:
    https://matplotlib.org/3.3.3/_modules/matplotlib/contour.html#ClabelText
    https://matplotlib.org/3.1.0/api/path_api.html#matplotlib.path.Path
    """
    if kinds is None:
        return segs    # nothing to be done
    # search for kind=79 as this marks the end of one polygon segment
    # Notes: 
    # 1. we ignore the different polygon styles of matplotlib Path here and only
    # look for polygon segments.
    # 2. the Path documentation recommends to use iter_segments instead of direct
    # access to vertices and node types. However, since the ipyleaflet Polygon expects
    # a complete polygon and not individual segments, this cannot be used here
    # (it may be helpful to clean polygons before passing them into ipyleaflet's Polygon,
    # but so far I don't see a necessity to do so)
    new_segs = []
    for i, seg in enumerate(segs):
        segkinds = kinds[i]
        boundaries = [0] + list(np.nonzero(segkinds == 79)[0])
        for b in range(len(boundaries)-1):
            new_segs.append(seg[boundaries[b]+(1 if b>0 else 0):boundaries[b+1]])
    return new_segs
        
# add contours as polygons
# hardwired colors for now: these correspons to the 8-level default of matplotlib with an added orange color
colors = ["#48186a", "#424086", "#33638d", "#26828e", "#1fa088", "#3fbc73", "#84d44b", "#d8e219", "#fcae1e"]
allsegs = cs.allsegs
allkinds = cs.allkinds

for clev in range(len(cs.allsegs)):
    kinds = None if allkinds is None else allkinds[clev]
    segs = split_contours(allsegs[clev], kinds)
    polygons = Polygon(
                    locations=[p.tolist() for p in segs],
                    # locations=segs[14].tolist(),
                    color=colors[min(clev, 4)],
                    weight=1,
                    opacity=0.8,
                    fill_color=colors[clev],
                    fill_opacity=0.5
    )
    map.add_layer(polygons);
map

结果如下所示：在此处输入图像描述 现在当然可以使用传单选项并使用不同的背景地图等。

于 2021-01-10T16:33:09.867 回答