I tried to extract a cross-section from a DEM. Somehow I get weird spikes every other point. I dont see this behavior in QGis. Somehow the extract function is not sorting the values acording to the location. The weird thing is, I get these spikes only for some cross-sections depending on the angle they are laying on the DEM. I used several different DEMs and Lines, I also tried the terra package, all with the same result. Here is the R-Code for the cross-section:

l.cut<-unlist(raster::extract(r.dem, l.trans)

Here a plot of the DEM and the two LINESTRIGS i want to use to generate the cross-section:

rel<-get_elev_raster(r.dem, z = 10)

enter image description here

Here is a graph from the profile tool from the southern cross-section in QGis:

Here is a graph from the profile tool in QGis

The same graph generated by the raster:extract function in R:

xl<-seq(from=0, to=as.integer(trans$length), by=as.numeric(trans.line$length/length(trans)))

          geom_line(aes(x=xl, y=l.cut), size=1)+
          scale_y_continuous(breaks = seq(trunc(min(l.cut)/10)*10, (trunc(max(l.cut)/10)+2)*10, by = 100))+

We see some weird spikes in the southern cross-section:

The same graph generated by the raster:extract function in R

In the northern cross-section we dont see these spikes:

enter image description here

  • 1
    Have you tried this using terra? The raster package is being depreciated and replaced by terra which has all the same functionality (same developer) just moved to C++. Even if this is a bug in raster, you would be directed to terra. To read data use terra::rast. There is a vector class in the package so to extract you would use extract(r.dem, vect(trans)) The result is a data.frame object, not a list so, you have to aggregate by ID using something like tapply. Commented Jun 19, 2023 at 19:17
  • Tried it with the suggested terra approach, with the same result.
    – sfetan
    Commented Jun 20, 2023 at 7:52
  • 1
    I was able to successfully create a cross-section profile. If you could provide a reproducible example or data we may be able to further troubleshoot. You could use dput to post your sf data. Is your line feature a LINESTRING or MULTILINESTRING? Is there only a single feature in the data. How did you create your plot (provide code please). Commented Jun 20, 2023 at 14:30
  • I updated the code for my graph. The line feature is a LINESTRING. I tried several DEM and the problem persists. One odd thing is, that this spikes only occurs when the line lays in a certain angle.
    – sfetan
    Commented Jun 20, 2023 at 16:59
  • 1
    I think I identified the issue and asked it as a separate question with as simple reproducible example here: gis.stackexchange.com/questions/465023/…
    – Adam C
    Commented Aug 10, 2023 at 18:08

1 Answer 1


I created several random profiles and cannot reproduce any type of "error" or erroneous results so, lets try to put your results into context. At a large distance, I would imagine that a cross-slope profile would exhibit these "spikes" since the elevations decreases and increases with the slope(s). I illustrate this effect in my example.

First, we add required libraries (note; spatialEco is only for the DEM), create 2 sf LINESTRING features representing down and cross slope profiles and, extract the elevation profiles which are then turned into a list object for plotting. We plot the profile lines on a shaded relief to better visualize the terrain.


elev <- rast(system.file("extdata/elev.tif", package="spatialEco"))

cross <- data.frame(id=1, cbind(c(613160.2, 645160.2),c(4918281, 4930281)), 
  names(cross)[2:3] <- c("x","y")
down <- data.frame(id=2, cbind(c(658160.2, 671160.2),c(4986281, 4959281)), 
  names(down)[2:3] <- c("x","y")

xy <- rbind(cross, down)  
  xy <- xy %>%
    st_as_sf(coords = c("x", "y"), crs = st_crs(elev)) %>%
      group_by(id) %>% 
        summarize() %>%
          filter(st_geometry_type(.) == "MULTIPOINT") %>%

z <- extract(elev, vect(xy))
  z <- lapply(unique(z$ID), \(i) as.numeric(z[z$ID == i,][,2])) 
    names(z) <- c("cross-slope", "down-slope")
slope <- terrain(elev, "slope", unit="radians")
aspect <- terrain(elev, "aspect", unit="radians")
hill <- shade(slope, aspect, 40, 270)
  plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
    plot(elev, col=rainbow(25, alpha=0.35), add=TRUE)
      plot(st_geometry(xy), lwd=2, add=TRUE)

shaded relief with profiles

Hare are comparisons of cross and down slope profiles. As you can see, the cross slope profile exhibits similar "spikes" as your result. This is due to the profile intersecting decreasing/increasing slopes.

dev.new(height=8, width=16)
lapply(1:length(z), \(i) {
  l = as.integer(st_length(xy[i,]))  
  x <- seq(from=0, to=l, by=as.integer(l/length(z[[i]])))[1:length(z[[i]])]
    plot(x, z[[i]], type="l", xlab="distance", 
         lwd=1.5, ylab="elev", main=paste0(names(z)[[i]], " profile"))

cross vs down slope

Now, as to the influence of smoothing. I have a feeling that the profile plot that you are getting from QGIS is being smoothed in someway (Python function is calling QwtPlotCurve). Often software will do things in the background that the user is unaware of, making it difficult to compare results from different software. Here is a plot that contains the original profile and the result of 3 different spline smoothing parameters. You can really see how the plot starts looking like the one produced by QGIS.

dev.new(height=8, width=16)
l = as.integer(st_length(xy[2,]))  
x <- seq(from=0, to=l, by=as.integer(l/length(z[[2]])))[1:length(z[[2]])]
  plot(x, z[[2]], type="l", xlab="distance", 
       lwd=1.5, ylab="elev", main=paste0(names(z)[[2]], " profile"))
    lines(smooth.spline(x, z[[2]], spar = 0.25), lty=2, col="red")
      lines(smooth.spline(x, z[[2]], spar = 0.4), col="darkgreen")
        lines(smooth.spline(x, z[[2]], spar = 0.6), col="blue")
    legend("bottomright", legend=c("raw", "s=2", "s=4", "s=6"),
           lty=c(1,2,1,1), col=c("black", "red", "darkgreen", "blue"))

cross-slope profile smoothing

  • I think you found the core issue, the extract approach only works on smooth landscapes. But I am still not hundred percent sure if there is not an issue with the angle of the line. I edited my question regarding your answer. I am not sure if QwtPlotCurve is not using a total different approach, the curve dosent seem very smoothed to me but more detailed. Next thing I will try is to convert the lines into points with very small distance and use them for a cross-section.
    – sfetan
    Commented Jun 23, 2023 at 16:31
  • I think I identified the issue and asked it as a separate question with as simple reproducible example here: gis.stackexchange.com/questions/465023/…
    – Adam C
    Commented Aug 10, 2023 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.