Search code examples
rgeospatial

How to create a best fit linestring through a set of spatial points using R


I have a set of unstructured coordinates along a coastline that I would like to approximate with a spatial linestring.

Coords:

x<-structure(list(anid = c(1, 1, 1, 1, 1, 1, 1,1,1,1,1,1,1, 1,1,1,1,1,1,1)
                  , longitude = c(-96.67,   -96.7,  -96.72, -96.5,  -96.97, -96.72, -96.55, -97.07, -96.48, -96.39, -96.45, -96.44, -96.98, -97.02, -96.78, -96.57, -96.36, -97.07, -96.53, -96.45)
                  , latitude = c(28.13, 28.11,  28.09,  28.27,  27.89,  28.09,  28.23,  27.81,  28.28,  28.36,  28.31,  28.32,  27.88,  27.85,  28.05,  28.21,  28.39,  27.8,   28.24,  28.31))
             , class = "data.frame", row.names = c(NA,   -20L))

The coords plotted

library(sf)
library(leaflet)
TX_coord <- st_as_sf(x, coords = c(2:3))
st_crs(TX_coord) <- 4326
leaflet(TX_coord) %>% addTiles() %>% addCircleMarkers()

enter image description here

I'd like to approximate these points with a simple linestring. However if I cast to a linestring it's a mess as the locations are unordered.

library(sfheaders)
lsTX<- sf_linestring( obj = x
                      , x = "longitude"
                      , y = "latitude"
                      , linestring_id = "anid"
                      , keep = T
)
leaflet(lsTX) %>% addTiles() %>% addPolylines(weight=3)

enter image description here

What I really want is just a spatial version of a best fit line through the points. Like sfit below but as some kinda spatial line that I can plot with leaflet or similar. I can't work out that step. Thanks.

sfit <- smooth.spline(x[,2] ~ x[,3], spar=0.80)
plot(x[,2], x[,3], pch=20, main="best fit through coordinates", 
     xlab="longitude", ylab="latitude", cex=0.50)
lines(sfit$y, sfit$x, lwd=1.25, col="blue") 

enter image description here


Solution

  • You could run smooth.spline with the coordinates from the spatial data and turn the results into a linestring:

    sfit <- smooth.spline(st_coordinates(TX_coord)[,1] ~ st_coordinates(TX_coord)[,2], spar=0.80)
    cbind(sfit$y, sfit$x) %>% st_linestring() %>% leaflet() %>% addTiles() %>% addPolylines()
    

    enter image description here