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:
- Table with xy coordinates
- Build a spatial object from table
- 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.
- Use same table with xy coordinates to do a polyline shapefile from a SpatialLinesDataFrame
- 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")
}
a
andb
. 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.tabl
as well otherwise themake.shape
function can't be called and we are no closer to being able to copy/paste the code.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"