I ran HDBSCAN on these coordinates and got some clusters but some are too large. HDBSCAN has a minimum cluster size parameter, but no maximum size. All I want is to intuitively divide larger clusters so that none are larger than a maximum cluster size. Essentially, something that creates the second plot below. I could rerun HDBSCAN on the cluster, but it isn't guaranteed to reduce the size. I can do k-means, but I'm using clustering on a shortest-path distance matrix, something that I don't think k-means can handle.
cl = gpd.read_file("cluster.gpkg")
df = pd.concat([cl.geometry.x, cl.geometry.y], axis = 1)
mins = 40
min_size = 400
clusterer = hdbscan.HDBSCAN(min_cluster_size=min_size, min_samples=mins, allow_single_cluster=True)
clusterer.fit(df.iloc[:, 0:2])
pd.Series(clusterer.labels_).value_counts()
df['cluster'] = clusterer.labels_
df['subcluster'] = np.nan
# Plot
plt.figure(figsize=(10,6))
unique_labels = np.unique(clusterer.labels_)
noise_data = df[df['cluster'] == -1]
fig, ax = plt.subplots(1, 2, figsize=(20,6)) # Create 1 row, 2 columns of plots
for label in unique_labels:
subset = df[df['cluster'] == label]
if label == -1:
ax[0].scatter(subset[0], subset[1], label="Noise", color='gray', marker=".", alpha=0.3)
else:
ax[0].scatter(subset[0], subset[1], label=f'Cluster {label}', marker = ".", alpha=1)
ax[0].set_xlabel('Longitude')
ax[0].set_ylabel('Latitude')
ax[0].set_title('Clusters')
ax[0].legend()

Have a look at this PR which adds to the source code this parameter. They seem to be what you are after though they are hard to find, as noted in this issue.
Hope this helps.