使用Matplotlib生成等高线后，如何导出含Z值的Shapefile？

修改后的完整代码 #

#Convert UTM list into numpy array         
list_utm_array = np.array(list_utm_rot, dtype = float)

#Create equal distance grid spacing     
# Create x,y,z arrays for gridding
x = list_utm_array[:,0]
y = list_utm_array[:,1]
z = list_utm_array[:,2]

# Determine min and max x,y,z values
min_x = int(min(x))
max_x = int(max(x))
min_y = int(min(y))
max_y = int(max(y))
min_z = int(min(z))
max_z = int(max(z)) * 1.2 # multiply by 1.2 for greater interpolated z values

# Create empty x and y axes for meshgrid
xi = np.linspace(min_x, max_x, 100)
yi = np.linspace(min_y, max_y, 100)

# Create a empty meshgrid 
X, Y = np.meshgrid(xi, yi)

# Interpolate rotated utm values onto the grid
from scipy.interpolate import griddata
Z = griddata((x, y), z, (X, Y), method='cubic')

#Create contour plot using Matlibplot
contour_set = plt.contourf(X, Y, Z, levels = 5)

# ------------------- 关键修改部分 -------------------
# 获取当前坐标轴上的所有等高线填充集合，以及对应的层级值
contour_collections = plt.gca().collections
levels = contour_set.levels

# 存储多边形和对应Z属性的列表
contour_data = []

# 遍历每个填充区域（对应一个Z区间）
for i, collection in enumerate(contour_collections):
    # 获取当前填充区域对应的Z区间上下限
    z_lower = levels[i]
    z_upper = levels[i+1]
    # 可根据需求选择存储的Z值（比如区间下限、上限或中间值）
    z_representative = z_lower
    
    # 提取当前填充区域的所有路径
    paths = collection.get_paths()
    for path in paths:
        # 将路径转为多边形（处理带洞的复杂多边形）
        for polygon_coords in path.to_polygons():
            polygon = Polygon(polygon_coords)
            # 将多边形和Z属性存入列表
            contour_data.append({
                'geometry': polygon,
                'z_lower': z_lower,
                'z_upper': z_upper,
                'z_value': z_representative
            })

# 创建带属性的GeoDataFrame
gdf = gpd.GeoDataFrame(contour_data, crs = "EPSG:32614")

# Save to shapefile
gdf.to_file("out_with_z.shp")

关键说明 #

1. 获取填充集合与层级值：

   • 新版Matplotlib中，QuadContourSet不再直接提供collections属性，需从当前坐标轴（plt.gca()）获取所有等高线填充集合。
   • contour_set.levels返回等高线的所有分界值，若设置levels=5，则会生成6个分界值、5个填充区域（对应5个collection）。

2. 关联Z值与多边形：

   • 每个collection对应一个Z值区间（levels[i]到levels[i+1]），可根据需求存储区间下限、上限或中间值作为属性。
   • 每个路径（path）可能包含多个多边形（比如带内洞的区域），因此需遍历path.to_polygons()的结果。

3. Shapefile属性支持：

   • GeoDataFrame中的z_lower、z_upper、z_value列会被导出为Shapefile的属性字段，可直接在GIS软件中用于分类、渲染等操作。

内容的提问来源于stack exchange，提问作者Rathpr