skimage.measure.find_contours 方法返回什么？

Question

我有一个包含点类型几何数据的 CSV 文件。我想获得这个地形模型的等高线。该算法很好地给了我结果，但我根本不理解结果。我无法绘制任何东西。

这是我的代码：

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

#data = pandas.read_csv(file, sep=",")
#data_geo = geopandas.GeoDataFrame(data, geometry=geopandas.points_from_xy(data.longitude, data.latitude))

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")
    
contours = measure.find_contours(grid_z.T, 1)
if len(contours) > 0:
    plt.figure()
    for contour in contours:
        plt.plot(contour[:, 1], contour[:, 0], linewidth=2, color="k")
        
else:
    print("Not values") 

我希望将等高线和我的点 (grid_z) 放在同一个图上，并以 GeoJSON 格式导出等高线。

Answer 1

扩展我上面的评论：

import numpy as np
import matplotlib.pyplot as plt

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]

xmin, xmax, ymin, ymax = 7.481173172, 7.481649023, 46.33465000, 46.33479000
x_resolution, y_resolution = 100, 100
dx = (xmax - xmin) / x_resolution # pixel width
dy = (ymax - ymin) / y_resolution # pixel height

grid_x, grid_y = np.meshgrid(
    np.linspace(xmin, xmax, x_resolution),
    np.linspace(ymin, ymax, y_resolution),
)

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

fig, ax = plt.subplots()
ax.imshow(grid_z.T, extent=(xmin, xmax, ymin, ymax), origin='lower', cmap="RdYlGn")

contours = measure.find_contours(grid_z.T, 1)

def map_to_geographic_coordinates(image_contour):
    rows, columns = image_contour.T
    return np.c_[columns * dx + xmin, rows * dy + ymin]

if len(contours) > 0:
    for image_contour in contours:
        geographic_contour = map_to_geographic_coordinates(image_contour)
        ax.plot(geographic_contour[:, 0], geographic_contour[:, 1], linewidth=1, color="k")

plt.show()