Sf: Interpolating locations?

[JimShady]

Hello,

Had a quick browser of the docs and couldn't find a way to fill missing geometries. Is there one? I'm thinking similar to this function except spatial.

data <- tibble(x     = c(10,NA,NA,NA,20),
               y     = c(5,NA,NA,NA,15)) %>%
        st_as_sf(coords = c('x', 'y'), na.fail = F)

data

Simple feature collection with 5 features and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 10 ymin: 5 xmax: 20 ymax: 15
epsg (SRID):    NA
proj4string:    NA
       geometry
1  POINT (10 5)
2   POINT EMPTY
3   POINT EMPTY
4   POINT EMPTY
5 POINT (20 15)

Thanks.

[tim-salabim]

You could use st_line_sample to achive this

library(sf)
library(tibble)

data <- tibble(x     = c(10,NA,NA,NA,20),
               y     = c(5,NA,NA,NA,15)) %>%
  st_as_sf(coords = c('x', 'y'), na.fail = F)

# create multipoint and cast to line
ln = st_cast(st_union(data), "LINESTRING")

# sample along line at 5 intervals
interp = st_line_sample(ln, sample = seq(0, 1, 0.2))

interp
#> Geometry set for 1 feature 
#> geometry type:  MULTIPOINT
#> dimension:      XY
#> bbox:           xmin: 10 ymin: 5 xmax: 20 ymax: 15
#> epsg (SRID):    NA
#> proj4string:    NA
#> MULTIPOINT (10 5, 12 7, 14 9, 16 11, 18 13, 20 15)

[JimShady]

Thanks @tim-salabim . I can see how that would work for the example I have given. I should actually have made my example more complicated to better illustrate my question though. My data looks more like this:

data <- tibble(x     = c(10,NA,NA,NA,20,NA,30,NA,NA,50,NA,NA,NA,NA,NA,10),
               y     = c(10,NA,NA,NA,20,NA,30,NA,NA,50,NA,NA,NA,NA,NA,10)) %>%
        st_as_sf(coords = c('x', 'y'), na.fail = F)

Simple feature collection with 16 features and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 10 ymin: 10 xmax: 50 ymax: 50
epsg (SRID):    NA
proj4string:    NA
First 10 features:
        geometry
1  POINT (10 10)
2    POINT EMPTY
3    POINT EMPTY
4    POINT EMPTY
5  POINT (20 20)
6    POINT EMPTY
7  POINT (30 30)
8    POINT EMPTY
9    POINT EMPTY
10 POINT (50 50)

Any ideas?

[tim-salabim]

You'd need to break the data into chunks. This will get you nearly there (I think), but there is still a duplicate point at the end (probably due to rounding issues):

data <- tibble(x     = c(10,NA,NA,NA,20,NA,30,NA,NA,50,NA,NA,NA,NA,NA,10),
               y     = c(10,NA,NA,NA,20,NA,30,NA,NA,50,NA,NA,NA,NA,NA,10)) %>%
  st_as_sf(coords = c('x', 'y'), na.fail = F)

idx = which(!st_is_empty(data))
st = idx[1:(length(idx)-1)]
nd = idx[2:length(idx)]

subsets = lapply(seq(st), function(i) {
  data[st[i]:nd[i], ]
})

interp = do.call(c, lapply(seq(subsets), function(i) {
  ln = st_cast(st_union(subsets[[i]]), "LINESTRING")
  int = st_line_sample(ln, sample = seq(0, 1, length.out = nrow(subsets[[i]])))
  return(int)
}))

interp = st_cast(interp, "POINT")

interp = interp[!duplicated(interp)]

interp

st_coordinates(interp)
#           X        Y
# 1  10.00000 10.00000
# 2  12.50000 12.50000
# 3  15.00000 15.00000
# 4  17.50000 17.50000
# 5  20.00000 20.00000
# 6  25.00000 25.00000
# 7  30.00000 30.00000
# 8  36.66667 36.66667
# 9  43.33333 43.33333
# 10 50.00000 50.00000
# 11 16.66667 16.66667
# 12 23.33333 23.33333
# 13 36.66667 36.66667

[JimShady]

Thanks @tim-salabim . It's a bit of a pain isn't it ! Great work with the code.

I've just been thinking I might be best converting the X and Y into normal numeric columns, then using this function, then converting back to a sf object.

[tim-salabim]

That's also an option :-)

[tim-salabim]

I'm closing this here, feel free to reopen if necessary

[JimShady]

Incase anyone finds it useful when future googling, I ended up doing this:

library(zoo)

data <- tibble(x  = c(10,NA,NA,NA,20,NA,30,NA,NA,50,NA,NA,NA,NA,NA,10),
               y  = c(10,NA,NA,NA,20,NA,30,NA,NA,50,NA,NA,NA,NA,NA,10)) %>%
                      st_as_sf(coords = c('x', 'y'), na.fail = F)

st_geometry(data) <- as_tibble(st_coordinates(data)) %>% 
                      mutate(x = na.approx(x), y = na.approx(y)) %>% 
                      st_as_sf(coords = c('x', 'y')) %>%
                      st_geometry()

Cheers.