关于skimage.measure.find_contours返回值及等高线绘制导出的技术问题

1. skimage.measure.find_contours的返回值 #

skimage.measure.find_contours返回一组numpy数组组成的列表，每个数组对应一条独立的等高线。数组内每一行的格式为(行索引, 列索引)，代表该点在输入网格（如你的grid_z.T）中的位置，并非直接的地理经纬度坐标。

2. 等高线与插值网格同图绘制 #

你的代码当前将网格图和等高线分别绘制在两个独立的画布中，要合并到同一张图，只需移除第二个plt.figure()调用，并将等高线的网格索引坐标转换为实际地理坐标：

修改后的完整绘制代码片段：

import rasterio.transform
import snowgis_hn_utils as hn
import glob
import os
import matplotlib.pyplot as plt
import pandas
import geopandas
import numpy
from shapely import Point
from scipy.interpolate import griddata
from skimage import measure

points = [[7.481479858, 46.33465909],
[7.481419827, 46.33467655],
[7.481358293, 46.33469339],
[7.481294507, 46.33471043],
[7.481231878, 46.33472677],
[7.481173172, 46.33474505],
[7.481648, 46.33465],
[7.4816235, 46.33466],
[7.4815875, 46.33467],
[7.481547, 46.33469],
[7.4815065, 46.33470],
[7.481467, 46.33471],
[7.4814275, 46.334715],
[7.48139, 46.33473],
[7.4813545, 46.33474],
[7.481321, 46.33475],
[7.481288, 46.33476],
[7.4812555, 46.33477],
[7.481221, 46.33478],
[7.481184, 46.33479],
[7.481649023, 46.33467052],
[7.481596643, 46.33468609],
[7.481545527, 46.33470014],
[7.481492227, 46.33471225]]
values = [1.3, 1.0, 0.95, 0.90, 0.80, 0.70, 1.40, 1.35, 1.20, 1.30, 1.40, 1.25, 1.35, 1.45, 1.45, 1.35, 1.30, 1.20, 1.25, 1.30, 1.05, 0.80, 0.90, 0.95]
    
grid_x, grid_y = numpy.mgrid[
    7.481173172:7.481649023:100j,
    46.33465000:46.33479000:100j
]

grid_z = griddata(points, values, (grid_x, grid_y), method='linear')

plt.figure()
# 绘制插值网格
plt.imshow(grid_z.T, extent=(7.481173172, 7.481649023, 46.33465000, 46.33479000), origin='lower', cmap="RdYlGn")

# 地理坐标范围
x_min, x_max = 7.481173172, 7.481649023
y_min, y_max = 46.33465000, 46.33479000
# 获取网格的行列数
n_cols, n_rows = grid_z.T.shape

contours = measure.find_contours(grid_z.T, 1)
if len(contours) > 0:
    for contour in contours:
        # 将网格索引转换为地理经纬度
        lon = x_min + (contour[:, 1] / (n_cols - 1)) * (x_max - x_min)
        lat = y_min + (contour[:, 0] / (n_rows - 1)) * (y_max - y_min)
        plt.plot(lon, lat, linewidth=2, color="k")
else:
    print("无对应等高线") 

plt.colorbar(label="高程值")
plt.xlabel("经度")
plt.ylabel("纬度")
plt.title("插值地形与等高线")
plt.show()

3. 导出等高线为GeoJSON #

借助geopandas将等高线转换为地理空间要素后导出：

# 在上述代码获取contours之后执行
import geopandas as gpd
from shapely.geometry import LineString

# 存储等高线的LineString对象
contour_geoms = []
x_min, x_max = 7.481173172, 7.481649023
y_min, y_max = 46.33465000, 46.33479000
n_cols, n_rows = grid_z.T.shape

for contour in contours:
    # 转换为地理坐标
    lon = x_min + (contour[:, 1] / (n_cols - 1)) * (x_max - x_min)
    lat = y_min + (contour[:, 0] / (n_rows - 1)) * (y_max - y_min)
    # 创建LineString
    line = LineString(zip(lon, lat))
    contour_geoms.append(line)

# 创建GeoDataFrame，指定坐标系为WGS84（EPSG:4326）
gdf = gpd.GeoDataFrame(
    data={"elevation": [1.0]*len(contour_geoms)},  # 记录等高线对应的高程
    geometry=contour_geoms,
    crs="EPSG:4326"
)

# 导出为GeoJSON文件
gdf.to_file("terrain_contours.geojson", driver="GeoJSON")

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