raster::rasterToContour; contour lines are not continuous

Question

I am trying to extract contour lines from a raster object using the raster package in R.

rasterToContour appears to work well and plots nicely, but when investigated it appears the contour lines are broken up into irregular segments. Example data from ?rasterToContour

library(raster)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
x <- rasterToContour(r)
class(x)
plot(r)
plot(x, add=TRUE)

I am trying to extract the contour line of a sample site in the raster. So, we choose a random site, extract its elevation, and run rasterToContour() again, specifying the elevation for the contour line level.

# our sample site - a random cell chosen on the raster
xyFromCell(r, 5000) %>%
  SpatialPoints(proj4string = crs(r)) %>%
  {. ->> site_sp} %>%
  st_as_sf %>%
  {. ->> site_sf}

# find elevation of sample site, and extract contour lines
extract(r, site_sf) %>%
  {. ->> site_elevation}

# extract contour lines
r %>%
  rasterToContour(levels = site_elevation) %>%
  {. ->> contours_sp} %>%
  st_as_sf %>%
  {. ->> contours_sf}

# plot the site and new contour lines (approx elevation 326)
plot(r)
plot(contours_sf, add = TRUE)
plot(site_sf, add = TRUE)

# plot the contour lines and sample site - using sf and ggplot
ggplot()+
  geom_sf(data = contours_sf)+
  geom_sf(data = site_sf, color = 'red')

Then we use st_intersects to find the lines that intersect the site (with a buffer width of 100 to ensure it touches the line). But, this returns all of the contour lines.

contours_sf %>%
  filter(
    st_intersects(., site_sf %>% st_buffer(100), sparse = FALSE)[1,]
  ) %>%
  ggplot()+
  geom_sf()

I assume all lines are returned because they appear to be a single MULTILINESTRING sf object.

contours_sf

# Simple feature collection with 1 feature and 1 field
# geometry type:  MULTILINESTRING
# dimension:      XY
# bbox:           xmin: 178923.1 ymin: 329720 xmax: 181460 ymax: 333412.3
# CRS:            +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs
# level                       geometry
# C_1 326.849822998047 MULTILINESTRING ((179619.3 ...

So, I have split the contours_sf MULTILINESTRING into individual lines using ngeo::st_segments() (I couldn't find any sf way to do this, but am open to using alternative methods, especially if this is the problem).

Unexpectedly this returns 394 features; from looking at the figure I would expect approximately 15 separate lines.

contours_sf %>%
  nngeo::st_segments()

# Simple feature collection with 394 features and 1 field
# geometry type:  LINESTRING
# dimension:      XY
# bbox:           xmin: 178923.1 ymin: 329720 xmax: 181460 ymax: 333412.3
# CRS:            +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs
# First 10 features:
#   level                         result
# 1  326.849822998047 LINESTRING (179619.3 329739...
# 2  326.849822998047 LINESTRING (179580 329720.4...
# 3  326.849822998047 LINESTRING (179540 329720, ...
# 4  326.849822998047 LINESTRING (179500 329735.8...
# 5  326.849822998047 LINESTRING (179495.3 329740...
# 6  326.849822998047 LINESTRING (179460 329764, ...
# 7  326.849822998047 LINESTRING (179442.6 329780...
# 8  326.849822998047 LINESTRING (179420 329810, ...
# 9  326.849822998047 LINESTRING (179410.2 329820...
# 10 326.849822998047 LINESTRING (179380 329847.3...

Then, when we filter to keep only the lines which intersect the site (with a buffer width of 100), only a small section of the expected contour line is returned (red section of line, I assume reflective of the 100 buffer width).

contours_sf %>%
  nngeo::st_segments() %>%
  filter(
    # this syntax used as recommended by this answer https://stackoverflow.com/a/57025700/13478749
    st_intersects(., site_sf %>% st_buffer(100), sparse = FALSE)
  ) %>%

  ggplot()+
  geom_sf(colour = 'red', size = 3)+
  geom_sf(data = contours_sf)+
  geom_sf(data = site_sf, colour = 'cyan')+
  geom_sf(data = site_sf %>% st_buffer(100), colour = 'cyan', fill = NA)

Anyone got ideas for the following points:

Explain why the contour lines are 'broken'
Provide an efficient method for 'joining' the broken pieces together
An alternative to nngeo::st_segments(), if this is in fact the source of the 394 lines not ~15

score 3 · Accepted Answer · answered May 18 '21 at 09:45

Converting the MULTILINESTRING to a LINESTRING seems to do what you need:

contours_sf %>% st_cast("LINESTRING") %>% 
  filter(st_intersects(., st_buffer(site_sf, 100), sparse=FALSE)[,1]) %>%

ggplot()+
  geom_sf(data = contours_sf)+
  geom_sf(data = site_sf, color = 'red') +
  geom_sf(color = 'pink')

Robert Hijmans · Answer 2 · 2021-05-19T04:47:51.233

Perhaps it works better if you start by disaggregating the lines

library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
x <- rasterToContour(r)
x <- disaggregate(x)

Or with terra

library(terra)
r <- rast(f)
x <- as.contour(r)

x
# class       : SpatVector 
# geometry    : lines 
# dimensions  : 8, 1  (geometries, attributes)

x <- disaggregate(x)
x
# class       : SpatVector 
# geometry    : lines 
# dimensions  : 56, 1  (geometries, attributes)

And you can continue like this

y <- st_as_sf(x)

Or like this

r <- rast(system.file("ex/meuse.tif", package="terra"))
site <- vect(xyFromCell(r, 5000), crs=crs(r))
elevation <- extract(r, site)
v <- disaggregate(as.contour(r, levels=elevation))
i <- which.min(distance(site, v))
vv <- v[i]

plot(r)
lines(v)
lines(vv, col="red", lwd=2)
points(site, col="blue", cex=2)

Thanks @Robert, both of these options work well – hugh-allan May 24 '21 at 05:05 — hugh-allan, May 24 '21 at 05:05

raster::rasterToContour; contour lines are not continuous

2 Answers2

Linked