Skip to content

Commit

Permalink
Handle case where too few points
Browse files Browse the repository at this point in the history
  • Loading branch information
jonaraphael committed Nov 10, 2023
1 parent b1dd4fc commit 88bdbc4
Showing 1 changed file with 34 additions and 27 deletions.
61 changes: 34 additions & 27 deletions cerulean_cloud/cloud_function_ais_analysis/utils/associate.py
Original file line number Diff line number Diff line change
Expand Up @@ -202,34 +202,41 @@ def slick_to_curves(
y_fit_sort_x = [c[1] for c in coords_unique_x]
y_fit_sort_y = [c[1] for c in coords_unique_y]

# fit a B-spline to the centerline
tck_sort_x, fp_sort_x, _, _ = scipy.interpolate.splrep(
x_fit_sort_x, y_fit_sort_x, k=3, s=smoothing_factor, full_output=True
)
tck_sort_y, fp_sort_y, _, _ = scipy.interpolate.splrep(
y_fit_sort_y, x_fit_sort_y, k=3, s=smoothing_factor, full_output=True
)

# choose the spline that has the lowest fit error
if fp_sort_x <= fp_sort_y:
tck = tck_sort_x
x_fit = x_fit_sort_x
y_fit = y_fit_sort_x

num_points = max(round((x_fit[-1] - x_fit[0]) / 100), 5)
x_new = np.linspace(x_fit[0], x_fit[-1], 10)
y_new = scipy.interpolate.BSpline(*tck)(x_new)
else:
tck = tck_sort_y
x_fit = x_fit_sort_y
y_fit = y_fit_sort_y

num_points = max(round((y_fit[-1] - y_fit[0]) / 100), 5)
y_new = np.linspace(y_fit[0], y_fit[-1], num_points)
x_new = scipy.interpolate.BSpline(*tck)(y_new)
# Check if there are enough points for spline fitting
min_points_required = 4 # for cubic spline, k=3, need at least 4 points
if len(coords_unique_x) >= min_points_required:
# fit a B-spline to the centerline
tck_sort_x, fp_sort_x, _, _ = scipy.interpolate.splrep(
x_fit_sort_x, y_fit_sort_x, k=3, s=smoothing_factor, full_output=True
)
tck_sort_y, fp_sort_y, _, _ = scipy.interpolate.splrep(
y_fit_sort_y, x_fit_sort_y, k=3, s=smoothing_factor, full_output=True
)

# store as LineString
curve = shapely.geometry.LineString(zip(x_new, y_new))
# choose the spline that has the lowest fit error
if fp_sort_x <= fp_sort_y:
tck = tck_sort_x
x_fit = x_fit_sort_x
y_fit = y_fit_sort_x

num_points = max(round((x_fit[-1] - x_fit[0]) / 100), 5)
x_new = np.linspace(x_fit[0], x_fit[-1], 10)
y_new = scipy.interpolate.BSpline(*tck)(x_new)
else:
tck = tck_sort_y
x_fit = x_fit_sort_y
y_fit = y_fit_sort_y

num_points = max(round((y_fit[-1] - y_fit[0]) / 100), 5)
y_new = np.linspace(y_fit[0], y_fit[-1], num_points)
x_new = scipy.interpolate.BSpline(*tck)(y_new)

# store as LineString
curve = shapely.geometry.LineString(zip(x_new, y_new))
else:
curve = shapely.geometry.LineString(
[coords_unique_x[0], coords_unique_x[-1]]
)
slick_curves.append(curve)

slick_curves_gdf = gpd.GeoDataFrame(geometry=slick_curves, crs=slick_gdf.crs)
Expand Down

0 comments on commit 88bdbc4

Please sign in to comment.