6

I have a SpatialLineDataFrame with 20 observations of 4 variables (in WGS84). I want to create waypoints every 10 meters. Every waypoint will has the x y coordinates. I tried to do it manually like:

point_first=SpatialPoints(cbind(x,y)) 

for the first point, but there are many points I want to do, so it is not possible to do it that way. I think I need a loop but I don't know how to do it in order to have the coordinates for every waypoint.

Is there any command that I should use?

7

There is a nice tutorial on segmentation of spatial lines which you can find here, and you could make use of CreateSegment to achieve your goal. But first of all, let's create some sample data.

## german borders
library(rgdal)
library(rworldmap)

data(countriesCoarse)
spy_germany <- subset(countriesCoarse, GEOUNIT == "Germany")

## polygons to lines
sln_germany <- as(spy_germany, "SpatialLinesDataFrame")

Since you'd like to work with meters rather than degrees, reproject the data (e.g. to UTM) and extract the coordinates which we will require later on as input for CreateSegment.

## reprojection (utm-32n) and extraction of coordinates
sln_germany_utm <- spTransform(sln_germany, CRS("+init=epsg:32632"))
mat_germany_crd <- coordinates(sln_germany_utm)[[1]][[1]]

Next, let's create a sequence representing the desired segmentation. rgeos::gLength comes in quite handy for that. I chose segments of 100 km in length.

## border length and sequence of segments (100,000 m = 100 km)
library(rgeos)
num_len <- gLength(sln_germany_utm)

int_len_seq <- seq(0, num_len, 100000)

You can now loop over this sequence and extract the start and end coordinate of each single segment.

## loop over segments of 100 km
source("CreateSegment.R")

ls_crd_seg <- lapply(2:length(int_len_seq), function(i) {

  # extract coordinates of current line segment
  mat_segment <- CreateSegment(mat_germany_crd, 
                               from = int_len_seq[i-1], to = int_len_seq[i])

  # end coordinate
  crd_out <- matrix(mat_segment[nrow(mat_segment), ], ncol = 2)

  # during the first iteration, also return the start coordinate
  if (i == 2) {
    crd_start <- matrix(mat_segment[1, ], ncol = 2)
    crd_out <- rbind(crd_start, crd_out)
  } 

  return(crd_out)
})

Finally, rbind the coordinates and assign the projection of our initial shapefile. That's it!

## coordinates to SpatialPoints
crd_seg <- do.call("rbind", ls_crd_seg)
spt_seg <- SpatialPoints(crd_seg)
proj4string(spt_seg) <- proj4string(sln_germany_utm)

## visualize
plot(sln_germany_utm, lwd = 1.5)
points(spt_seg)

segmented_borders

  • This is a clever approach! – Jeffrey Evans Jul 2 '15 at 15:31
  • Very nice approach! – geo_dd Jul 3 '15 at 8:14
6

This function will create a systematic point sample of a SpatialLinesDataFrame. The LID column in the results corresponds to the row.names of the source line object and can be used to relate the points to the lines.

   # Function for creating a systematic point sample of a SpatialLinesDataFrame
   #   x        SpatialLinesDataFrame
   #   sdist    Distance of sample spacing 
   #   offset   offset for start of line
   #   longlat  TRUE/FALSE, is data in geographic units. If TRUE distance is 
   #              great circle with units (sdist argument) in kilometers.    
    systematic.lineSample <- function(x, sdist=100, offset = 0.5, longlat = FALSE)
       {
        if (!require(sp)) stop("sp PACKAGE MISSING")
         if (!inherits(x, "SpatialLinesDataFrame")) stop("MUST BE SP SpatialLinesDataFrame OBJECT")
        lids <- rownames(x@data)    
        lsub <- x[rownames(x@data) == lids[1] ,]
         LLength <- SpatialLinesLengths(lsub, longlat = longlat)
          ns <- round( (LLength[1] / sdist), digits=0)
          lsamp <- spsample(lsub, n=ns, type="regular", offset=c(offset, offset))
          results <- SpatialPointsDataFrame(lsamp, data=data.frame(LID=rep(as.numeric(lids[1]), ns))) 
        if( length(lids) > 1) { 
          for (i in 2:length(lids) ) 
            {    
             lsub <- x[rownames(x@data) == lids[i] ,]
              LLength <- SpatialLinesLengths(lsub, longlat = longlat)
               ns <- round( (LLength[1] / sdist), digits=0)
               lsamp <- spsample(lsub, n=ns, type="regular")
             lsamp <- SpatialPointsDataFrame(lsamp, data=data.frame(LID=rep(as.numeric(lids[i]), ns)))
             results <- rbind(results, lsamp)     
           }
        }
      ( results )
    }

Here is an example of the functions usage

# Create example data
require(sp)
sp.lines <- SpatialLines(list(Lines(list(Line(cbind(c(1,2,3),c(3,2,2)))), ID="1"),
                         Lines(list(Line(cbind(c(1,2,3),c(1,1.5,1)))), ID="2")))
sp.lines <- SpatialLinesDataFrame( sp.lines, data.frame(ID=1:2, row.names=c(1,2)) )

# Systematic line sample
lsample <- systematic.lineSample(sp.lines, sdist = 20, offset = 0, longlat = TRUE) 

# Plot results
plot(sp.lines)
  plot(lsample, pch = 20, add = TRUE)
  • 1
    I haven't heard of sp::SpatialLinesLengths before, thanks for that! Always better to rely on less packages... anyway, when I tried systematic.lineSample(sln_germany_utm, sdist = 100000), an error was thrown during the for loop since only one row name was found (i.e. the initial polygon record). Any way to fix this? – fdetsch Jul 2 '15 at 13:47
  • @fdetsch Done, all the code needed was an if statement proceeding to the loop if there is more than one line. Thanks for pointing this out. – Jeffrey Evans Jul 2 '15 at 15:25
  • Note; I also added a longlat argument to indicate geographic units. If the data is in geographic units (longlat=TRUE) then the distance units are in kilometers. – Jeffrey Evans Jul 2 '15 at 15:51
2

The SpatStat package has a function for creating points along lines at specified intervals.

# maptools to read shapefile
require(maptools)
require(spatstat)

# Load a line shapefile
myLine <- readShapeLines("~/workspace/TEMP/myLine.shp")
plot(myLine)

enter image description here

# create a PSP object to use with spatstat
myPSP <- as.psp(myLine)

# Create points at 10 km intervals
myPoints <- pointsOnLines(myPSP, 10000)
plot(myPoints, add = T)

enter image description here

  • 1
    Note that the point placement using pointsOnLines is realized for each segment separately (see documentation) instead of considering the line as a whole. Applying your code to my example (with eps = 100000) results in a far too narrow point placing, e.g. seen from the small jags along the eastern and southern border. – fdetsch Jul 2 '15 at 13:26
  • Aha, had never used this on a multi-segment line (only transects). – Simbamangu Jul 3 '15 at 8:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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