消除轮廓起点和终点之间的连接线

Question

我一直在与 ChatGPT 进行斗争，以通过 PyCharm IDE 使用 python 3.12 来生成代码，消除从声纳数据派生的轮廓线的起始和结束顶点之间的连接线（可通过

sonarlight

包访问）。

我是一个没有经验的 py 用户，已经为此奋斗了几个小时，寻找解决方案。

这是不想要的结果：

这是产生结果的当前代码：

import sys
import site

# Add your custom package path
custom_path = 'C:/Users/Rion/AppData/Local/Programs/Python/Python312/Lib/site-packages'
sys.path.append(custom_path)

import geopandas as gpd
from shapely.geometry import Polygon, LineString
from sonarlight import Sonar
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
from pyproj import Transformer
import matplotlib.contour
import pandas as pd

# Load sonar data from two files and merge them
df1 = Sonar('../data/Sonar_2024-09-03_18.08.00.sl3').df
df2 = Sonar('../data/Sonar_2024-09-03_16.12.35.sl3').df

# Merge the two DataFrames
df = pd.concat([df1, df2], ignore_index=True)

# Step 1: Compute the elevation column by subtracting water_depth from gps_altitude
df['elevation'] = df['gps_altitude'] - df['water_depth']

# Step 2: Convert lat/lon to EPSG:32736 (WGS 84 / UTM zone 36S)
transformer = Transformer.from_crs("epsg:4326", "epsg:32736", always_xy=True)
df['x'], df['y'] = transformer.transform(df['longitude'], df['latitude'])

# Step 3: Generate grid for contouring
x = df['x'].values
y = df['y'].values
z = df['elevation'].values

# Create a grid for interpolation
xi = np.linspace(x.min(), x.max(), 500)
yi = np.linspace(y.min(), y.max(), 500)
xi, yi = np.meshgrid(xi, yi)

# Interpolate the elevation data onto the grid
zi = griddata((x, y), z, (xi, yi), method='linear')

# Step 4: Create contours (both polygons and lines)
contour_levels = np.arange(z.min(), z.max(), 2)  # Adjust contour intervals as needed

# Generate contour polygons
plt.figure()
contour_poly = plt.contourf(xi, yi, zi, levels=contour_levels)

# Generate contour lines
plt.figure()
contour_line = plt.contour(xi, yi, zi, levels=contour_levels)

# Convert matplotlib contour objects to shapely geometries for polygons
contour_polygons = []
poly_elevations = []  # To store corresponding elevation for each polygon
for collection in contour_poly.collections:
    level = collection.get_label() if hasattr(collection, 'get_label') else collection.level
    for path in collection.get_paths():
        vertices = path.vertices
        poly = Polygon(vertices)
        contour_polygons.append(poly)
        poly_elevations.append(level)

# Convert matplotlib contour objects to shapely geometries for lines with discontinuity handling
contour_lines = []
line_elevations = []  # To store corresponding elevation for each line
discontinuity_threshold = 500  # Threshold for identifying discontinuities (adjust if necessary)
for collection in contour_line.collections:
    level = collection.get_label() if hasattr(collection, 'get_label') else collection.level
    for path in collection.get_paths():
        vertices = path.vertices
        # Calculate the differences between successive vertices (distance)
        diffs = np.sqrt(np.sum(np.diff(vertices, axis=0) ** 2, axis=1))
        # Insert NaN where the difference is larger than the threshold (discontinuity)
        discontinuity_indices = np.where(diffs > discontinuity_threshold)[0]
        if len(discontinuity_indices) > 0:
            new_vertices = []
            last_index = 0
            for idx in discontinuity_indices:
                new_vertices.extend(vertices[last_index:idx + 1])
                new_vertices.append([np.nan, np.nan])  # Insert NaN to break the line
                last_index = idx + 1
            new_vertices.extend(vertices[last_index:])
            vertices = np.array(new_vertices)

        # Create the LineString, ignoring NaN values (discontinuity)
        cleaned_vertices = vertices[~np.isnan(vertices).any(axis=1)]
        if len(cleaned_vertices) > 1:
            line = LineString(cleaned_vertices)
            contour_lines.append(line)
            line_elevations.append(level)

# Step 5: Create GeoDataFrames for polygons and lines with 'elevation' as an attribute
gdf_poly = gpd.GeoDataFrame({'elevation': poly_elevations, 'geometry': contour_polygons}, crs="EPSG:32736")
gdf_lines = gpd.GeoDataFrame({'elevation': line_elevations, 'geometry': contour_lines}, crs="EPSG:32736")

# Step 6: Set the index (row ID) to the 'elevation' value
gdf_poly.set_index('elevation', inplace=True)
gdf_lines.set_index('elevation', inplace=True)

# Step 7: Export to shapefiles without clipping
gdf_poly.to_file("../results/blydepoort_bathymetry_20240903_poly.shp")
gdf_lines.to_file("../results/blydepoort_bathymetry_20240903_contours.shp")

Answer 1

好吧，虽然我没有在OP上提供可重现的示例，但我设法通过修改后的代码产生了所需的结果，如下所示，感谢@Cincinnatus：

import sys
import site

# Add your custom package path
custom_path = 'C:/Users/Rion/AppData/Local/Programs/Python/Python312/Lib/site-packages'
sys.path.append(custom_path)

import geopandas as gpd
from shapely.geometry import LineString
from sonarlight import Sonar
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
from pyproj import Transformer
import pandas as pd

# Load sonar data from two files and merge them
df1 = Sonar('../data/Sonar_2024-09-03_18.08.00.sl3').df
df2 = Sonar('../data/Sonar_2024-09-03_16.12.35.sl3').df

# Merge the two DataFrames
df = pd.concat([df1, df2], ignore_index=True)

# Step 1: Compute the elevation column by subtracting water_depth from gps_altitude
df['elevation'] = df['gps_altitude'] - df['water_depth']

# Step 2: Convert lat/lon to EPSG:32736 (WGS 84 / UTM zone 36S)
transformer = Transformer.from_crs("epsg:4326", "epsg:32736", always_xy=True)
df['x'], df['y'] = transformer.transform(df['longitude'], df['latitude'])

# Step 3: Generate grid for contouring
x = df['x'].values
y = df['y'].values
z = df['elevation'].values

# Create a grid for interpolation
xi = np.linspace(x.min(), x.max(), 500)
yi = np.linspace(y.min(), y.max(), 500)
xi, yi = np.meshgrid(xi, yi)

# Interpolate the elevation data onto the grid
zi = griddata((x, y), z, (xi, yi), method='linear')

# Step 4: Create contours (lines only)
contour_levels = np.arange(z.min(), z.max(), 2)  # Adjust contour intervals as needed

# Generate contour lines
plt.figure()
contour_line = plt.contour(xi, yi, zi, levels=contour_levels)

# Convert matplotlib contour objects to shapely geometries for lines with discontinuity handling
contour_lines = []
line_elevations = []  # To store corresponding elevation for each line

# Use updated method to access contour collections
for i in range(len(contour_line.levels)):
    level = contour_line.levels[i]
    for path in contour_line.allsegs[i]:  # Use allsegs to get the correct line segments
        line = LineString(path)
        contour_lines.append(line)
        line_elevations.append(level)

# Step 5: Create GeoDataFrame for lines only with 'elevation' as an attribute
gdf_lines = gpd.GeoDataFrame({'elevation': line_elevations, 'geometry': contour_lines})

# Step 6: Set the index (row ID) to the 'elevation' value
gdf_lines.set_index('elevation', inplace=True)

# Step 7: Export to shapefile without polygons
gdf_lines.to_file("../results/blydepoort_bathymetry_20240903_contours.shp")

虽然不是代码或主要问题的一部分，但还是进行了后验。

@Cincinnatus 这对于一些人来说可以接受学习吗？