4

I have a table with xy coordinates in R and the following packages loaded: spatstat, raster, rgdal, sp and zoo. I do know how to use this to make a point shapefile, a line shapefile that goes through these points and a smoothed line shapefile using rollmean on my original xy coordinates. My workflow it is something similar to this:

  1. Table with xy coordinates
  2. Build a spatial object from table
  3. save the shapefile shape points giving the name of the spatial object, the directory where the file will be saved and the name of the file without extension.
  4. Use same table with xy coordinates to do a polyline shapefile from a SpatialLinesDataFrame
  5. Using rollmean do a smoothing of the xy coordinates, and save it as a smoothed polyline shapefile.

Ok, now I have all the objects I want, including SpatialLines and SpatialLinesDataFrame for the smoothed coordinates.

Now I want to do the following 2 main transformations using R (I know how to do that if I use ArcGIS but I don't want to go back and forth from R to ArcGIS):

a. Transform the smoothed shapefile (or any of the above "smoothed" objects) into points at 1 m apart, and add to the shapefile an attribute table with xy coordinates for each point along the smoothed line

b. Transform the smoothed shapefile in p points equally distributed along the smoothed line (that means for me equal distance – although I don't know what that distance is) and add to the attribute table xy coordinates for each point.

# make shapefiles - point, line, smoothed line
#tabl is a table with the first 2 columns x and y coordinates which is what is used 
#to actually do the shapefile
#rast - the raster from which we take the projection and the resolution
# smooth.win = how big you want the smoothing window, the real smooth param is actually the 
# no. of cells you want to smooth over
#dir.name = path of the directory where i want to save the shape files.
#Need to use double back slash when writing the directory path name.
# file name – a name that will be added to the name of the raster that is analysed

make.shape <- function(tabl, dir.name, file.name, smooth.win, rast){

require(spatstat)
require(raster)
require(rgdal)
require(sp)
require(zoo)

crs <- projection(rast)
res.x<-round(xres(rast),1)
res.y <- round(yres(rast),1)
rast.name <- names(rast)
file.name <- paste(rast.name, "_", file.name, sep = "")

if(is.data.frame(tabl)) tabl <- tabl else tabl <- as.data.frame(tabl)

# build spatial objects that can be saved as ESRI shapefiles
shape.sp <- SpatialPointsDataFrame(tabl[,1:2], tabl, proj4string=CRS(crs))

# save the shapefile shape points giving the name of the spatial object,
#the directory where the file will be saved and the name of the file
#without extension. 

writeOGR(shape.sp, dir.name, paste(file.name, "_points", sep = ""), 
   driver = "ESRI Shapefile")

# PolyLine shapefile

xy <- cbind(tabl[,1], tabl[,2])
cl <- Line(xy)
cl1 <- Lines(list(cl), ID = 1)
row.names(tabl) = seq(1, dim(tabl)[1])

shape.l <- SpatialLines(list(cl1), proj4string = CRS(crs))
shape.ln <- SpatialLinesDataFrame(shape.l, tabl)

writeOGR(shape.ln, dir.name, paste(file.name, "_line", sep = ""), 
     driver = "ESRI Shapefile")

# do the smoothing of the shape line
rmy <- rollmean(tabl[,2], smooth.win)
rmx <- rollmean(tabl[,1], smooth.win)

x1y1 <- matrix(c(rmx,rmy), ncol=2, nrow = length(rmy), byrow = FALSE)
cl1a <- Line(x1y1)
cl1b <- Lines(list(cl1a), ID = 1)
row.names(x1y1) <- seq(1, dim(x1y1)[1])
shape.l1 <- SpatialLines(list(cl1b), proj4string = CRS(crs))
x1y1 <- as.data.frame(x1y1)
shape.ln1 <- SpatialLinesDataFrame(shape.l1, x1y1)

writeOGR(shape.ln1, dir.name, paste(file.name, "_smooth_line_",
     round(smooth.win*res.x,0), "m", sep = ""), driver = "ESRI Shapefile")
}
  • What have you tried so far? If you have the code for steps 1 to 5 then it would be good to have it - we can use it as the base to attempt the next two steps a and b. Remember that a minimal, "fake" example is all we need, so you could create a deliberately simplified version of steps 1 to 5 if you feel that helps. – SlowLearner Dec 17 '13 at 21:00
  • Hi,I am posting my R code to do shapefiles when using a table or a matrix that has xy coordinates in the first 2 columns. – Monica Palaseanu-Lovejoy Dec 17 '13 at 21:31
  • Can you supply some (fake or real) data for tabl as well otherwise the make.shape function can't be called and we are no closer to being able to copy/paste the code. – SlowLearner Dec 18 '13 at 8:33
  • OK, here it is a little sample of my data. I will provide rezolution and projection string instead of the raster rast - since i have no idea how to provide that. Just run the code under the function declaration that the code as a function itsel. It will work - hopefully. Code next cooment. OK the code looks very bad as a comment - i hope you can make something out of it ..... it is not clear to me how to inset code in comments section – Monica Palaseanu-Lovejoy Dec 18 '13 at 19:56
  • Xi <- c(590798.8, 590796.8, 590796.8, 590796.9, 590796.9, 590794.9, 590792.9, 590793.0, 590793.0, 590793.0, 590791.0 5,90791.0, 590791.1, 590791.1, 590789.1, 590789.1, 590789.2, 590789.2, 590789.2, 590789.2) Yi <- c(4147491, 4147492, 4147493, 4147494, 4147495, 4147496, 4147497, 4147498, 4147499, 4147500, 4147501, 4147502, 4147504, 4147505, 4147506, 4147507, 4147508, 4147509, 4147510, 4147511) tabl <- cbind(Xi, Yi) crs <- "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0" res.x <- 1 res.y <- 1 rast.name <- "my_raster" file.name <- - "lines_to_points" – Monica Palaseanu-Lovejoy Dec 18 '13 at 19:57
6

I am a bit unclear on your question but if you have an sp line object then you can use "spsample" in a loop to create a systematic distance based sample for each line. Here is a function that implements a distance based sample along a line.

The arguments are x = sp SpatialLinesDataFrame and sdist = Sampling distance.

# Function for systematic distance line sample
sample.line <- function(x, sdist=100)
   {
    if (!require(sp)) stop("sp PACKAGE MISSING")
     if (!inherits(x, "SpatialLinesDataFrame")) stop("MUST BE SP SpatialLinesDataFrame OBJECT")
      lgth <- SpatialLinesLengths(x) 
      lsub <- x[1,]
        ns <- round( (lgth[1] / sdist), digits=0)
          lsamp <- spsample(lsub, n=ns, type="regular", offset=c(0.5,0.5))
            results <- SpatialPointsDataFrame(lsamp, data=data.frame(ID=rep(rownames(x@data[1,]),ns))) 
    for (i in 2:dim(x)[1] ) 
      {    
       lsub <- x[i,]
         ns <- round( (lgth[i] / sdist), digits=0)
           lsamp <- spsample(lsub, n=ns, type="regular")
       lsamp <- SpatialPointsDataFrame(lsamp, data=data.frame(ID=rep(rownames(x@data[i,]),ns)))
       results <- rbind(results, lsamp)     
     }
  ( results )
}

require(sp)
# Create some data
l1 = cbind(c(1,2,3),c(3,2,2))
  l1a = cbind(l1[,1]+.05,l1[,2]+.05)
    l2 = cbind(c(1,2,3),c(1,1.5,1))
Sl1 = Line(l1)
  Sl1a = Line(l1a)
    Sl2 = Line(l2)
      S1 = Lines(list(Sl1, Sl1a), ID="a")
        S2 = Lines(list(Sl2), ID="b")
sp.lines = SpatialLines(list(S1,S2))
  ldf <- data.frame(ID=c(1,2), Y=c(0.5,0.25), row.names=c("a","b"))
    sp.lines <- SpatialLinesDataFrame(sp.lines, ldf)   

# Run sample function and display results   
lsamp <- sample.line(sp.lines, sdist=0.5)   
  plot(sp.lines, col = c("red", "blue"))
    plot(lsamp, pch=20, add=TRUE)
  • 1
    Thank you so much Jeffrey. This is a very nice answer. I discovered that if i have only 1 polyline shapefile, only the first part of the function is needed. Your solution is certainly much more elegant than mine. I have problems with fully understanding objects in sp, although i do use them. One interesting observation is that my code uses the first vertex of the polyline as the first point in the new generated points, while your results are exactly inbetween my points. Both codes give same number of points, yours being the most elegant solution. Thanks again. Monica – Monica Palaseanu-Lovejoy Dec 20 '13 at 14:49
  • If i modify the offset parameter to c(0,0) then the points obtained from this function are identical with the points i obtained previously through my code or by XToolsPro. Again thanks so much. Monica. – Monica Palaseanu-Lovejoy Dec 20 '13 at 19:24
1

I think i solved the problem in the end. I realized that using my function "make.shape" i actually have access to the xy coord of the smoothed line vertices.

That means i can calculate each distance between each pair of vertices, and the slope of each segment in-between pairs of vertices. Knowing that the equation of a straight line is y = m*x + b i can calculate b for each segment as well.

To move from a quadratic equation to a linear one i decided that actually i can project everything on a horizontal line, get my new x coordinates on this line and for each segment knowing m and b i can calculate the corresponding y. For now i have a huge for loop to go through pairs of vertices, but it runs pretty fast, although it is somewhat slower than the ArcGIS tool that does this (xToolsPro). I hope this will help somebody else who is faced with a similar problem.

Thanks, Monica

-1

@Jeffrey Evans has a nice idea but his function does not keep attribute of line to point shp. Please see my solution which is made some adjustments based his script.

### Function for systematic distance line sample nguồn: Stakeflow
sample.line <- function(x, sdist=100)
{
  if (!require(sp)) stop("sp PACKAGE MISSING")
  if (!inherits(x, "SpatialLinesDataFrame")) stop("MUST BE SP SpatialLinesDataFrame OBJECT")
  lgth <- SpatialLinesLengths(x) 
  lsub <- x[1,]
  ns <- round( (lgth[1] / sdist), digits=0)
  lsamp <- spsample(lsub, n=ns, type="regular", offset=c(0.5,0.5))
  data = x@data[1,]
  df = as.data.frame(lapply(data, rep, ns))
  results <- SpatialPointsDataFrame(lsamp, df, match.ID = TRUE)

  for (i in 2:dim(x)[1] ) 
  {    
    lsub <- x[i,]
    ns <- round( (lgth[i] / sdist), digits=0)
    lsamp <- spsample(lsub, n=ns, type="regular")
    data = x@data[i,]
    df = as.data.frame(lapply(data, rep, ns))
    lsamp <- SpatialPointsDataFrame(lsamp, df, match.ID = TRUE)
    results <- rbind(results, lsamp)     
  }
  ( results )
}

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.