How do I apply functions to each element of a SpatialLinesDataFrame, in particular using other data also stored as an attribute (e.g. coordinates or elevation per line vertex)?
I've a number of SpatialLine objects which are grouped together in a single SpatialLinesDataFrame (SLDF) object with an "Id" name. Using the SLDF, I can extract elevations to the points along each line, and calculate things like the length of each line as in the example below.
library(raster)
library(sp)
library(rgeos)
# made up elevations somewhere in England
dtm <- raster(xmn=412115, xmx=439324, ymn=352217, ymx=391748)
dtm <- setValues(dtm,runif(ncell(dtm),min=1,max=150))
projection(dtm)<-CRS("+init=epsg:27700")
# create some lines as SLDFs
x1 <- spLines(rbind(c(420000,360000), c(420300,365000), c(420600,370000)),crs=("+init=epsg:27700"),attr=data.frame(id=1))
x2 <- spLines(rbind(c(425000,360000), c(425300,365000), c(430600,370000)),crs=("+init=epsg:27700"),attr=data.frame(id=2))
list_sldf=c(x1,x2) # in the main code I'm using, there are 100s of thesefirst added to a list
merged.lines <- do.call(rbind, list_sldf) # combine to single SLDF
# extract elevation to each line
merged.lines$elev = raster::extract(dtm, merged.lines)
# calculate path length of each line
merged.lines$path_length=rgeos::gLength(merged.lines, byid=TRUE)
I now want to calculate additional metrics for each line, initally cumulative length per vertex and then eventually elevation specific derivatives e.g. along line curvature extract.
As an example, I've tried the following for calcualting cumulative distance (giving me an x axis to my "y" axis elevations):
If I try looping (just to start - ideally I want to vectorise this) the spatial dataframe to access each row, then I can't add new variables. For example, the function I have below to calculate the cumulative length between the vertices of each line is based on the coordinates of a line (e.g. cordinates[merged.lines[1,]
) but I can't work out how to apply this sequentially across the SLDF.
#' Calcualte euclidean distance between 2 points
#' @example euc.dist(c(0,0),c(2,2))
euc.dist <- function(x1, x2){
sqrt(sum((x1 - x2) ^ 2))
}
#' @param xys - a matrix of 2 columns (x/y)
#' @example
#' xys=merged.lines@lines[[1]]@Lines[[1]]@coords # a matrix
#' cumDist(xys)
cumulative.Dist<-function(xys,x='x',y='y', longlat=FALSE){
colnames(xys)<-c(x,y)
Dist <- 0
for(i in 1:length(xys[,x])){
if (i==1){
print("first point")
Dist[i]=0
}
else{
print("other points")
Dist[i] = euc.dist(c(xys[i,x],xys[i,y]), c(xys[i-1,x],xys[i-1,y]))
}
}
DistCumulative <- 0
for(i in 2:length(xys[,x])) {
DistCumulative[i] = Dist[i] + DistCumulative[i-1]
}
#xys_cumulativeDist <- cbind(xys, Dist, DistCumulative)
return(DistCumulative)
}
This works to get the values on an individual basis:
xy_df_=as.data.frame(coordinates(merged.lines[1,]))
cumulative.Dist(xy_df_)
How do I get it to run over each line in the SLDF?
Can I vectorise this?
for (i in 1:length(merged.lines)){
merged.lines[i,]$cumulative_length=cumulative.Dist(as.data.frame(coordinates(merged.lines[i,])), x='X1', y='X2')
}
Perhaps I am setting the variable incorrectly.
1 Answer 1
I was hoping to be able to vectorise the process but for now, I can achieve what is needed by looping the slots of the SpatialLinesDataFrame (one for each line).
Also, extracting elevation using merged.lines$elev = raster::extract(dtm, merged.lines)
causes issues if you then want to write out the SpatialLinesDatFrame to a shapefile
(I didn't check any other extensions). This appears to be because the elevation values become an attribute of the @data
slot.
My solution is below and was helped following this Q/A.
At each iteration, I get hold of the coordinates of the line and use them to calculate what I need which are then add summary values as attributes to the @data slot.
The position of the value in the relevant attribute list within the @data slot corresponds to the line e.g. 1st Id value in data at merged.lines@data$id[1]
relates to the first line object i.e. merged.lines@lines[1]
# Create data
library(raster)
library(sp)
library(rgeos)
# made up elevations somewhere in England
dtm <- raster(xmn=412115, xmx=439324, ymn=352217, ymx=391748)
dtm <- setValues(dtm,runif(ncell(dtm),min=1,max=150))
projection(dtm)<-CRS("+init=epsg:27700")
# create some lines as SLDFs
x1 <- spLines(rbind(c(420000,360000), c(420300,365000), c(420600,370000)),crs=("+init=epsg:27700"),attr=data.frame(id=1))
x2 <- spLines(rbind(c(425000,360000), c(425300,365000), c(430600,370000)),crs=("+init=epsg:27700"),attr=data.frame(id=2))
list_sldf=c(x1,x2)
# combine to single SLDF
merged.lines <- do.call(rbind, list_sldf)
# Check structure to see how things are laid out
str(merged.lines)
# Apply function to each line in SpatialLineDataFrame and add as new data slot
for (i in 1:length(merged.lines@lines)){
# get coords
xys=merged.lines@lines[[i]]@Lines[[1]]@coords
# get elevations from line vertex coordinates
elev=raster::extract(dtm, xys)
# Add as new attribute (a summary attribute for this specific line instance)
#merged.lines@data$length[[i]]=max(cumulative.Dist(xys)) # << calculate distance along line (see function in original question)
merged.lines@data$mean_elev[[i]]=mean(elev)
merged.lines@data$elev_rng[[i]]=max(elev)-min(elev)
}
# Check structure of SLDF again to see how it's changed
str(merged.lines)
#Write to file
opath="s4_testing.shp"
writeOGR(obj=merged.lines, dsn=opath, layer="merged.lines", driver="ESRI Shapefile")
When you then access the shapefile, each line within the vector file (in this case a shapefile) will have it's own specific attributes.