R/06makeLado.R

In darrennorris/parcelareadev: To test and develop parcel area calculations

#' makelado
#' 
#' @description
#' Should not be called directly by users.
#' \code{makelado} creates buffers to the left and right of central line.
#'
#' @param list.lines List with central lines (class SpatialLines) generated by \code{doline}.
#' @param faixa_lado Vector of plot survey widths.
#' @param lado_epsg EPSG code to use.
#'
#' @return List with SpatialPolygons (buffers) to left and right of central line.
#' @export
#'
#' @examples
#' \dontrun{
#' }
makelado <- function(list.lines, faixa_lado = c(0.5, 1, 1.5, 2, 10, 20),
                     lado_epsg = 3395) {
  # http://stackoverflow.com/questions/5726297/cut-polygons-using-contour-line-beneath-the-polygon-layers
require(sp)
require(rgdal)
require(rgeos)
  
  mycrs <- paste("+init=epsg:", lado_epsg,sep="")
  
l <- list.lines$SpatialLines_proj
mynames <- paste("ladobuff", names(l),sep="_")
names(l) <- mynames
lsplit <- list.lines$SpatialLines_proj$SpatialLinesAll
lsec <- subset(lsplit, seg_id ==1)

myline <- function(x) {
  require(sp)
  require(rgdal)
  require(rgeos)
  lsec2 <- gLineMerge(x, byid=FALSE, id = NULL)
lfull <- gLineMerge(lsplit, byid=FALSE, id = NULL)
dfb <- data.frame(bufnames = paste("ladobuf_", faixa_lado,"m", sep=""),
                  faixa_dist = faixa_lado)
mb <- function(x, myl = lsec, myl2 = lfull, myl3 = lsec2){
  require(sp)
  require(rgdal)
  require(rgeos)
  gb <- gBuffer(myl, width=x$faixa_dist, capStyle="FLAT")
  lpi2 <- gIntersection(gb, lsec)
  blpi2 <- gBuffer(lpi2, width = 0.000001)
  dpi2 <- gDifference(gb, blpi2)
  p1 <- Polygons(list(dpi2@polygons[[1]]@Polygons[[1]]), "1")
  p2 <- Polygons(list(dpi2@polygons[[1]]@Polygons[[2]]), "2")
  SpP <- SpatialPolygons(list(p1,p2), 1:2)                      
  mydata2 <- data.frame(id_lado = getSpPPolygonsIDSlots(SpP))
  rownames(mydata2) <- mydata2$id_lado
  spdf2 <- SpatialPolygonsDataFrame(SpP, data= mydata2)
  proj4string(spdf2) <- CRS(mycrs)
  
  spdf2wgs84 <- spTransform(spdf2, CRS("+init=epsg:4326"))
  lsplitwgs84 <- spTransform(myl, CRS("+init=epsg:4326"))
  # extract 20 points, uniform
  sp1 <- spsample(spdf2wgs84, n = 20, type = "stratified")
  sp1 <- SpatialPointsDataFrame(sp1, data = data.frame(id_lado = over(sp1,spdf2wgs84)))
  proj4string(sp1) <- CRS("+init=epsg:4326")
  
  dfs <- data.frame(myl)[1,]
  dfp <-  data.frame(plot_id = dfs$plot_id, long_x = dfs$long_x, lat_y = dfs$lat_y,
                     lead_x = dfs$lead_x, lead_y = dfs$lead_y ,id_lado = sp1$id_lado,
                     point_x = coordinates(sp1)[, 1], point_y = coordinates(sp1)[, 2]
  )
  dfp$aid <- row.names(dfp)
  
  pdist <- function(x){
    #source("05dist2gcS.R")
    myp <- parcelareadev::dist2gcS(c(x$long_x, x$lat_y), c(x$lead_x, x$lead_y), 
                    c(x$point_x, x$point_y), r=6378137)
    dfout <- data.frame(p_dist = myp)
    dfout
  }
  
  dout <- plyr::ddply(dfp, c("plot_id", "aid", "id_lado",
                            "point_x", "point_y"),
                      .fun = pdist)
  dout$lado <- ifelse(dout$p_dist < 0 ,"left", "right")
  
  dout <- plyr::ddply(dout, c("plot_id", "lado"), plyr::summarise, 
                point_x = mean(na.omit(point_x)), 
                point_y = mean(na.omit(point_y)) 
  )
  
  coordinates(dout) <- ~point_x + point_y
  proj4string(dout) <- CRS("+init=epsg:4326")
  dout <- spTransform(dout, CRS(mycrs))
  dout <- gBuffer(dout, width= 0.2, byid=TRUE)
    
  # symetric for left and right
  gb <- gBuffer(myl2, width= x$faixa_dist, capStyle="FLAT", byid=FALSE)
  lpi <- gIntersection(gb, myl2)
  blpi <- gBuffer(lpi, width = 0.000001)
  dpi <- gDifference(gb, blpi)
  
  if (length(dpi@polygons[[1]]@plotOrder) == 2) {
    require(sp)
    require(rgdal)
    require(rgeos)
    p1 <- Polygons(list(dpi@polygons[[1]]@Polygons[[1]]), "1")
    p2 <- Polygons(list(dpi@polygons[[1]]@Polygons[[2]]), "2")
    SpP <- SpatialPolygons(list(p1,p2), 1:2)                      
    mydata <- data.frame(id_lado = getSpPPolygonsIDSlots(SpP))
    rownames(mydata) <- mydata$id_lado
    spdfull <- SpatialPolygonsDataFrame(SpP, data= mydata)  
    proj4string(spdfull) <- CRS(mycrs)
    spdfull$plot_id <- over(spdfull, dout)[, 1]
    spdfull$lado <- over(spdfull, dout)[, 2]
    
    # asimetric for correct size and shape
    gba <- gBuffer(myl3, width= x$faixa_dist, capStyle="FLAT", byid=FALSE)
    lpia <- gIntersection(gba, myl3)
    blpia <- gBuffer(lpia, width = 0.000001)
    dpia <- gDifference(gba, blpia)
       
    require(raster)
    spdf_clip <- raster::intersect(spdfull, dpia)
    
    require(maptools)
    spdf_union <- unionSpatialPolygons(spdf_clip, IDs = spdf_clip$lado)
    mydata <- data.frame(lado = getSpPPolygonsIDSlots(spdf_union))
    rownames(mydata) <- mydata$lado
    spdf_union <- SpatialPolygonsDataFrame(spdf_union, data= mydata) 
    spdf_union$mycol <- ifelse(spdf_union$lado=="left","grey80", "grey60" )
    proj4string(spdf_union) <- CRS(mycrs)
    
    } else {spdf_union <- NA}
  
  lado_buff <- spdf_union
}

lado_all <- plyr::dlply(dfb, c("bufnames"), .fun = mb)
}

outbuflado <- lapply(l, FUN = myline)

}