Creating bordering polygons from spatial point data for plotting

I have had a good hunt around and can't find a solution to what seems to be a tricky problem.

I have basic XY (coordinate) data:

What I want to do is create neighbouring polygons based on the coordinate data that do not overlap and have a certain size limit (so they don't extend forever into the ocean).

Forgive my poor MS Paint skills, but the desired outcome would be something like this:

I have a polygon that marks the land/sea interface so polygons cannot overlap that either.

I am using Leaflet to make these maps interactive, its not for any statistical analysis but to provide an overview.

The ultimate aim is to have each polygon coloured by a variable (e.g. Temperature) with ecological data overlaid.

Some example data:

    > data[1:10,]
   Station  Lat_dec  Long_dec Surface_T
1      247 50.33445 -2.240283     15.19
2      245 50.58483 -2.535217     14.11
3      239 50.16883 -2.509250     15.41
4      225 50.32848 -2.765967     15.34
5      229 50.63900 -2.964800     14.09
6      227 50.33757 -3.303217     15.12
7      217 50.16657 -3.563817     15.13
8      207 49.66683 -3.556550     15.04
9      213 50.16512 -3.824667     14.97
10     219 49.83707 -3.815483     14.78

the Code to produce figure 1 is a basic leaflet script:

leaflet() %>% 
  addProviderTiles('Esri.OceanBasemap'
  ) %>% 
  addCircleMarkers(data = data,
                   lng = ~Long_dec,
                   lat = ~Lat_dec,
                   radius = 2
  ) %>%
  addPolygons(data = Land,
              weight = 1,
              color = 'black')

Most of the examples use downloaded polygons (e.g. what seems to be the classic US States rather than making them).

Here is something to get started with:

library(sf)
library(dplyr)

#create sf object with points
stations <- st_as_sf( df, coords = c( "Long_dec", "Lat_dec" ) ) 

#create voronoi/thiessen polygons
v <- stations %>% 
  st_union() %>%
  st_voronoi() %>%
  st_collection_extract()

library(leaflet)
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers( data = stations ) %>%
  addPolygons( data = v )