Устранить соединительные линии между начальными и конечными вершинами контуров. ⇐ Python

[Anonymous]

Я боролся с ChatGPT за создание кода, который устранял бы соединительные линии между начальными и конечными вершинами контурных линий, полученных из данных сонара (доступных через пакет sonarlight), используя Python 3.12 через PyCharm IDE.< /p>
Я неопытный пользователь 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")


Подробнее здесь: https://stackoverflow.com/questions/790 ... f-contours