Merge branch 'master' of https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee into devel

Add master updates to devel branch

.Rbuildignore

-^.*\.Rproj$
-^\.Rproj\.user$

.gitignore

 .Rproj.user
 .Rhistory
 .RData
+.Rbuildignore
 .Ruserdata
+lidaRtRee.Rproj
+lidaRtRee.Rcheck/
+lidaRtRee_*
+*~
+..Rcheck/
\ No newline at end of file

R/common.R

@@ -233,7 +233,7 @@ speciesColor <- function()
              c("Pinus mugo","salmon2", "C"),
              c("Pinus nigra","salmon1", "C"),
              c("Pinus sp.","salmon1", "C"),
-             c("Pinus uncinata","salmon1", "C"),
+             c("Pinus uncinata","salmon4", "C"),
              c("Pinus sylvestris","salmon", "C"),
              c("Populus alba","slateblue", "B"),
              c("Populus nigra","slateblue", "B"),

R/coregistration.R

@@ -375,7 +375,7 @@ coregistration <- function(chm, trees, mask= NULL, buffer=19, step=0.5, dm=2, pl
     # display initial tree positions
     graphics::points(trees[,1], trees[,2], cex=trees[,3]/40)
     # display coregistered tree positions
-    graphics::points(trees[,1]+resul$dx1, trees[,2]+resul$dy1, cex=trees[,3]/20, col="red")
+    graphics::points(trees[,1]+resul$dx1, trees[,2]+resul$dy1, cex=trees[,3]/40, col="red")
     graphics::legend("topleft", c("Initial", "Coregistered"), pch=1, col=c("black", "red"))
   }
   list(correlation.raster=r.correlation, local.max=resul)

R/tree.metrics.R

+#' # package lidaRtRee
+#' # Copyright Irstea
+#' # Author(s): Jean-Matthieu Monnet
+#' # Licence: LGPL-3
+#' ###################### FUNCTIONS FOR TREE-LEVEL METRICS COMPUTATION #################
+#' #
+#' ##############################
+#' #' segment-wise computation of point metrics
+#' #'
+#' #' @param p A dataframe with point attributes, including the segment ID field
+#' #' @param columnId A string: name of the field containing the segment ID
+#' #' @param FUN Function to compute for each segment ID. Default: number of points
+#' #' @return A vector of values obtainted by applying the function to the point cloud in each segment
+#' #' @export
+#' computeMetric <- function(p,columnId,FUN=nrow)
+#' {
+#'   by(p,as.factor(p[,columnId]),FUN)
+#' }
+#' ##############################
+#' #' replace values in a raster based on segment ID and tree attributes
+#' #'
+#' #' @param r.dem.w a raster with segment ID, integers starting from 0 to max(r.dem.w)
+#' #' @param segms A dataframe with tree attributes
+#' #' @param attr the attribute to use for rast
+#' #' @return a raster
+#' #' @export
+#' rasterValueFromSegment <- function(r.dem.w, segms, attr)
+#' {
+#'   dummy <- r.dem.w
+#'   # creer le tableau de correspondance id / attribut
+#'   # tableau avec autant de lignes que de segments 'max + 1' en comptant le 0
+#'   fh <- matrix(nrow=max(raster::values(r.dem.w))+1,ncol=1)
+#'   # mettre la hauteur correspondant aux id
+#'   fh[segms$id+1] <- segms[,attr]
+#'   # convertir id en h dans segmentation
+#'   raster::values(dummy) <- as.vector(fh[raster::values(dummy)+1])
+#'   dummy
+#' }
+#' ###############################
+#' #' add trunk altitude information
+#' #'
+#' #' @param segms a data.frame of segment id
+#' #' @param dtm a raster with terrain altitude
+#' #' @param p a point cloud data.frame
+#' #' @return a list with two elements: the segments with additional attribute (terrain altitude in the cell below the maximum) and the point cloud with additional attribute (point height relatively to terrain altitude below maximum)
+#' #' @export
+#' pointHeightAboveTrunk <- function(segms, dtm, p)
+#' {
+#'   # extract terrain altitude
+#'   segms$alt.dtm <- pointsInSegments(segms[,c("x","y")],dtm)
+#'   dummy <- merge(p,segms[,c("id","alt.dtm")],by.x="seg.id",by.y="id",all.x=T)
+#'   dummy$h.trunk <- dummy$z-dummy$alt.dtm
+#'   list(segms,dummy)
+#' }
\ No newline at end of file

lidaRtRee.Rcheck/00_pkg_src/lidaRtRee/DESCRIPTION

-Package: lidaRtRee
-Version: 1.0
-Title: Forest Analysis with Airborne Laser Scanning (Lidar) Data
-Date: 2018-04-27
-Author: Jean-Matthieu Monnet [aut, cre]
-Maintainer: Jean-Matthieu Monnet <jean-matthieu.monnet@irstea.fr>
-Description: Provides functions for forest analysis using airborne laser scanning data. It includes complementary steps for foret mapping: extraction of both physical and statistical features from lidar data, model calibration with ground reference, and maps export. 
-URL: https://gitlab.irstea.fr/jean-matthieu.monnet/lidaRtRee
-Imports: graphics, stats, grDevices, sp, raster, imager, akima, leaps,
-        gvlma, car, foreach, doParallel, reldist, lidR (>= 1.4-2)
-License: LGPL-3
-LazyData: TRUE
-RoxygenNote: 6.0.1
-NeedsCompilation: no
-Packaged: 2018-04-27 11:56:04 UTC; jean-matthieu
-Depends: R (>= 2.10)

lidaRtRee.Rcheck/00_pkg_src/lidaRtRee/NAMESPACE

-# Generated by roxygen2: do not edit by hand
-
-export(circle2Raster)
-export(coregistration)
-export(createDisk)
-export(demFiltering)
-export(heightRegression)
-export(histDetection)
-export(histStack)
-export(maximaDetection)
-export(maximaSelection)
-export(plot2Dmatched)
-export(plotTreeInventory)
-export(points2DSM)
-export(points2DTM)
-export(polar2Projected)
-export(rasterChullMask)
-export(rasterLocalmax)
-export(rasterXYMask)
-export(rasters2Cor)
-export(rastersMovingCor)
-export(segAdjust)
-export(segmentation)
-export(speciesColor)
-export(treeExtraction)
-export(treeMatching)
-export(treeSegmentation)

lidaRtRee.Rcheck/00_pkg_src/lidaRtRee/R/chmchablais3-data.R

-#' @name chmchablais3
-#' 
-#' @title Canopy height model (Chablais 3 plot)
-#' 
-#' @description Canopy height model computed from airborne laser scanning data acquired in July 2010.
-#' 
-#' @docType data
-#'
-#' @usage data(chmchablais3)
-#'
-#' @format A raster object
-#'
-#' @keywords datasets
-#'
-#' @references Monnet, J.-M. 2011. Using airborne laser scanning for mountain forests mapping: Support vector regression for stand parameters estimation and unsupervised training for treetop detection. Ph.D. thesis. University of Grenoble, France. pp. 21-22 & 34 \url{https://tel.archives-ouvertes.fr/tel-00652698/document}
-#'
-#' @examples
-#' data(chmchablais3)
-#' chmchablais3
-#' \dontrun{
-#' raster::plot(chmchablais3)}
-NULL
-"chmchablais3"

lidaRtRee.Rcheck/00_pkg_src/lidaRtRee/R/common.R

-# package lidaRtRee
-# Copyright Irstea
-# Author(s): Jean-Matthieu Monnet
-# Licence: LGPL-3
-################################
-#' Digital Surface Model
-#' 
-#' Creates a Digital Surface Model from LAS object. Raster extent is specified by the coordinates of lower left and upper right corners. Default extent covers the full range of points, and aligns on multiple values of the resolution. Cell value is the maximum height of points contained in the cell.
-#'
-#' @param .las LAS object (e.g. from \code{\link[lidR]{LAS}} function)
-#' @param res numeric. raster resolution
-#' @param xmin numeric. lower left corner easting coordinate for output raster.
-#' @param xmax numeric. upper right corner easting coordinate for output raster.
-#' @param ymin numeric. lower left corner northing coordinate for output raster.
-#' @param ymax numeric. upper right corner northing coordinate for output raster.
-#' @return A raster object.
-#' @seealso \code{\link{points2DTM}} for Digital Terrain Model computation.
-#' @examples
-#' data(laschablais3)
-#' 
-#' # create digital surface model with first-return points, resolution 0.5 m
-#' dsm <- points2DSM(lidR::lasfilterfirst(laschablais3), res=0.5)
-#' 
-#' \dontrun{
-#' # display raster
-#' raster::plot(dsm,asp=1)}
-#' @export
-points2DSM <- function(.las, res=1, xmin, xmax, ymin, ymax)
-{
-  #
-  if (missing(xmin) | missing(xmax) | missing(ymin) | missing(ymax)) # if no extent info
-  {
-    xmin <- floor(min(.las@data$X)/res)*res
-    xmax <- ceiling(max(.las@data$X)/res)*res
-    ymin <- floor(min(.las@data$Y)/res)*res
-    ymax <- ceiling(max(.las@data$Y)/res)*res
-  }
-  # create empty raster
-  r <- raster::raster()
-  raster::extent(r) <- c(xmin,xmax,ymin,ymax)
-  raster::res(r) <- res
-  raster::crs(r) <- NA
-  # convert LAS coordinates to spatial data
-  points <- as.data.frame(.las@data[,1:2])
-  val <- .las@data[,3]
-  sp::coordinates(points)=c(1,2)
-  # rasterize with max function
-  raster::rasterize(points,r,val,fun=max)
-}
-###################################
-#' Digital Terrain Model
-#' 
-#' Creates a Digital Terrain Model from LAS object. Raster extent is specified by the coordinates of lower left and upper right corners. Default extent covers the full range of points, and aligns on multiple values of the resolution. Cell value is compute as the Delaunay interpolation at the cell center. Relies on the \code{\link[akima]{interp}} function, which is slow, while waiting for new release of \code{\link[lidR]{grid_terrain}} which will integrate an additional argument to specify the output extent.
-#'
-#' @param .las LAS object (e.g. from \code{\link[lidR]{LAS}} function) containing only ground points
-#' @param res numeric. raster resolution
-#' @param xmin numeric. lower left corner easting coordinate for output raster.
-#' @param xmax numeric. upper right corner easting coordinate for output raster.
-#' @param ymin numeric. lower left corner northing coordinate for output raster.
-#' @param ymax numeric. upper right corner northing coordinate for output raster.
-#' @seealso \code{\link{points2DSM}} for Digital Surface Model computation.
-#' @return A raster object
-#' @examples
-#' data(laschablais3)
-#' 
-#' # create digital terrain model with points classified as ground
-#' dtm <- points2DTM(lidR::lasfilter(laschablais3, Classification==2))
-#' 
-#' \dontrun{
-#' # display raster
-#' raster::plot(dtm,asp=1)}
-#' @export
-points2DTM = function(.las, res = 1, xmin, xmax, ymin, ymax)
-{
-  if (missing(xmin) | missing(xmax) | missing(ymin) | missing(ymax)) # if no extent info
-  {
-    xmin <- floor(min(.las@data$X)/res)*res
-    xmax <- ceiling(max(.las@data$X)/res)*res
-    ymin <- floor(min(.las@data$Y)/res)*res
-    ymax <- ceiling(max(.las@data$Y)/res)*res
-  }
-    points <- as.matrix(.las@data[,c("X","Y","Z")])
-    # interpolation: value estimated at pixel center by bilinear interpolation
-    mnt <- akima::interp(points[,1],points[,2],points[,3], xo=seq(xmin+res/2, xmax-res/2,by=res), yo=seq(ymax-res/2,ymin+res/2,by=-res),linear=TRUE,extrap=FALSE, duplicate="user", dupfun=function(x){min(x)})
-    #
-  dtm = raster::raster(t(mnt[[3]]), xmn = xmin, xmx = xmax, ymn = ymin, ymx = ymax)
-  return(dtm)
-}
-# ################################
-# #' creates Digital Elevation Model from LAS object and extent info (parallelised version)
-# #' tiles extent and calls points2DTM or points2DSM
-# #'
-# #' @param .las LAS object (from \code{lidR} package)
-# #' @param type model to return ("dsm" or "dtm")
-# #' @param res the raster resolution
-# #' @param xmin the raster lower left corner easting coordinate
-# #' @param xmax the raster upper right corner easting coordinate
-# #' @param ymin the raster lower left corner northing coordinate
-# #' @param ymax the raster upper right corner northing coordinate
-# #' @param method a string specifying the interpolation method
-# #' @param buffer the size of buffer area for tiling
-# #' @param tile.size the size of tiles for parallelisation
-# #' @param n.cores the number of cores to use
-# #' @return a raster object
-# #' @export
-# points2DEM = function(.las, type="dtm", res = 1, xmin, xmax, ymin, ymax, method="delaunay", buffer=5, tile.size=250, n.cores=4)
-# {
-#   doParallel::registerDoParallel(cores=n.cores)
-#   if (missing(xmin) | missing(xmax) | missing(ymin) | missing(ymax)) # if no extent info
-#   {
-#     xmin <- floor(min(.las@data$X)/res)*res
-#     xmax <- ceiling(max(.las@data$X)/res)*res
-#     ymin <- floor(min(.las@data$Y)/res)*res
-#     ymax <- ceiling(max(.las@data$Y)/res)*res
-#   }
-#   # low left corner
-#   xlow <- floor(xmin/tile.size)*tile.size
-#   ylow <- floor(ymin/tile.size)*tile.size
-#   # tiles indices which intersect the extent
-#   i <- 0
-#   j <- 0
-#   n <- 1
-#   l <- list()
-#   while(xlow+i*tile.size<xmax)
-#   {
-#     while(ylow + j * tile.size < ymax)
-#     {
-#       l[[n]] <- c(i,j)
-#       n <- n+1
-#       j <- j+1
-#     }
-#     i <- i+1
-#     j <- 0
-#   }
-#   l <- foreach::foreach(i=1:length(l), .errorhandling="remove") %dopar%
-#   {
-#     coord <- c(max(xmin, xlow+l[[i]][1]*tile.size), min(xmax, xlow+(l[[i]][1]+1)*tile.size), max(ymin, ylow+l[[i]][2]*tile.size), min(ymax, ylow+(l[[i]][2]+1)*tile.size))
-#     dummy <- .las %>% lidR::lasfilter(X >= coord[1]-buffer & X <= coord[2]+buffer & Y >= coord[3]-buffer & Y <= coord[4]+buffer)
-#     if (type=="dtm") {return(points2DTM(dummy, res = res, coord[1], coord[2], coord[3], coord[4], method=method))}
-#     else {return(points2DSM(dummy, res = res, coord[1], coord[2], coord[3], coord[4]))}
-#   }
-#   l$fun <- mean
-#   l$na.rm <- TRUE
-#   do.call(raster::mosaic, l)
-# }
-################################
-#' Polar to cartesian coordinates conversion
-#' 
-#' Computes projected coordinates (Easting, Northing, Altitude) from polar coordinates (Azimuth, Slope, Distance) and center position (Easting, Northing, Altitude). Magnetic declination and meridian convergence are optional parameters. In case distance is measured to the border of objects (e.g. trees), the diameter can be added to compute the coordinates of object center.
-#'
-#' @param x vector. easting coordinates of centers in meter
-#' @param y vector. northing coordinates of centers in meter
-#' @param z vector. altitudes of centers in meters
-#' @param declination vector. magnetic declination values in radian
-#' @param convergence vector. meridian convergence values in radian
-#' @param azimuth vector. azimuth values from centers in radian
-#' @param slope vector. slope values from centers in radian
-#' @param dist vector. distances between centers and objects in meter
-#' @param diameter vector. diameters in meter (e.g. in case a radius should be added to the distance)
-#' @seealso \code{\link{plotTreeInventory}} for tree inventory display
-#' @return A data.frame with easting, northing and altitude coordinates, and horizontal distance from centers to objects centers
-#' @examples
-#' # create data.frame of trees with polar coordinates and diameters
-#' trees <- data.frame (x=rep(c(0,10), each=2),
-#'  y = rep(c(0,10),each=2),
-#'  z=rep(c(0,2),each=2),
-#'  azimuth=rep(c(0,pi/3)),
-#'  dist=rep(c(2,4)),
-#'  slope=rep(c(0,pi/6)),
-#'  diameter.cm=c(15,20,25,30))
-#' trees
-#'  
-#'  # compute projected coordinates
-#' polar2Projected(trees$x, trees$y, trees$z, trees$azimuth, trees$dist,
-#'  trees$slope, declination=0.03, convergence=0.02, trees$diameter.cm/200)
-#' @export
-polar2Projected <- function(x, y, z=0, azimuth, dist, slope=0, declination=0, convergence=0, diameter=0)
-{
-  d <- dist*cos(slope)+diameter/2
-  data.frame(x = x + d*sin(azimuth+convergence+declination),
-             y = y + d*cos(azimuth+convergence+declination),
-             z = z + (dist*sin(slope)+diameter/2),
-             d)
-}
-###############################
-#' Table of species names, abreviations and display colors
-#' 
-#' table for species names, abreviations and type (coniferous/broadleaf), and display color
-#' @return A data frame with species name, color, coniferous (C) / broadleaf (B) type, and name abreviation GESP of GEnus and SPecies 
-#' @seealso \code{\link{plotTreeInventory}} for tree inventory display
-#' @examples
-#' # load table
-#' tab.species <- speciesColor()
-#' head(tab.species)
-#' summary(tab.species)
-#' @export
-speciesColor <- function()
-{
-  d <- rbind(c("Abies alba","purple", "C"),
-             c("Acer","orange", "B"),
-             c("Acer campestre","orange1", "B"),
-             c("Acer opalus","orange2", "B"),
-             c("Acer platanoides","orange3", "B"),
-             c("Acer pseudoplatanus","orange4", "B"),
-             c("Acer sp.","orange", "B"),
-             c("Alnus incana","violet", "B"),
-             c("Alnus viridis","violet", "B"),
-             c("Betula pendula","darkgreen", "B"),
-             c("Betula pubescens","darkgreen", "B"),
-             c("Betula sp.","darkgreen", "B"),
-             c("Buxus sempervirens","black", "B"),
-             c("Carpinus betulus","seagreen", "B"),
-             c("Castanea sativa","pink", "B"),
-             c("Corylus avellana","cadetblue2", "B"),
-             c("Cornus mas","black", "B"),
-             c("Cotinus sp.","black", "B"),
-             c("Crataegus sp.","black", "B"),
-             c("Crataegus monogyna","black", "B"),
-             c("Euonymus latifolius","black", "B"),
-             c("Fagus sylvatica","green", "B"),
-             c("feuillus","black", "B"),
-             c("Fraxinus excelsior","yellow", "B"),
-             c("Ilex aquifolium","cyan", "B"),
-             c("inconnu","black", NA),
-             c("Juniperus communis","black", "C"),
-             c("Juniperus sp.","black", "C"),
-             c("Laburnum anagyroides","black", "B"),
-             c("Larix decidua","pink", "C"),
-             c("Larix kaempferi","pink", "C"),
-             c("Malus sylvestris","plum", "B"),
-             c("Picea abies","blue", "C"),
-             c("Pinus cembro","salmon3", "C"),
-             c("Pinus mugo","salmon2", "C"),
-             c("Pinus nigra","salmon1", "C"),
-             c("Pinus sp.","salmon1", "C"),
-             c("Pinus sylvestris","salmon", "C"),
-             c("Populus alba","slateblue", "B"),
-             c("Populus nigra","slateblue", "B"),
-             c("Populus tremula","slateblue", "B"),
-             c("Pseudotsuga menziesii","darkblue", "C"),
-             c("Prunus avium","grey", "B"),
-             c("Prunus sp.","grey", "B"),
-             c("Pyrus communis","grey", "B"),
-             c("Quercus petraea","turquoise", "B"),
-             c("Quercus pubescens","turquoise", "B"),
-             c("Quercus robur","turquoise", "B"),
-             c("Salix caprea","darkgoldenrod2", "B"),
-             c("Salix sp.","darkgoldenrod2", "B"),
-             c("Salix nigra","darkgoldenrod2", "B"),
-             c("Sorbus aria","red3", "B"),
-             c("Sorbus aucuparia","red", "B"),
-             c("Taxus baccata","burlywood3", "C"),
-             c("Tilia cordata","chocolate4", "B"),
-             c("Tilia platyphyllos","chocolate4", "B"),
-             c("Tilia sp.","chocolate4", "B"),
-             c("Ulmus glabra","brown", "B"),
-             c("Ulmus sp.","brown", "B"),
-             c("Viburnum lantana","black", "B"))
-  d <- as.data.frame(d,stringsAsFactors = F)
-  names(d) <- c("name","col","broad.conif")
-  d$broad.conif <- factor(d$broad.conif)
-  # abreviation GEsp (GEnus species)
-  dummy <- strsplit(d$name," ")
-  dummy2 <- lapply(dummy, function(x){
-    ifelse(length(x)==2 && x[2]=="sp.", 
-           paste(toupper(substr(x[1],1,2)),"sp",sep=""),
-           paste(toupper(substr(x,1,2)),collapse=""))})
-  d$abvr <- row.names(d) <- unlist(dummy2)
-  d
-}
-#######################################
-#' Displays a map of tree inventory data
-#' 
-#' displays tree inventory data
-#'
-#' @param xy data.frame. contains two columns with the X, Y  coordinates of tree centers
-#' @param height vector.  tree heights in meters 
-#' @param diam vector. tree diameters in centimeters
-#' @param species vector. species abreviation as in \code{\link{speciesColor}} for display with corresponding color.
-#' @param add boolean. indicates whether points should be added to an existing plot
-#' @seealso \code{\link{speciesColor}} for a table of species and associated colors
-#' @examples
-#' # load tree inventory data from plot Chablais 3
-#' data("treeinventorychablais3")
-#' 
-#' \dontrun{
-#' # display tree inventory
-#' plotTreeInventory(treeinventorychablais3[,c("x","y")],
-#' treeinventorychablais3$h,
-#' species=as.character(treeinventorychablais3$s))}
-#' @return no return
-#' @export
-#'
-plotTreeInventory <-function(xy, height=NULL, diam=NULL, species=NULL, add=FALSE)
-{
-  # set size of symbol
-  # proportionnal to height if present
-  if(!is.null(height))
-  {
-    size <- height/20
-    # otherwise proportionnal to diameter
-  } else {
-    if (!is.null(diam/40))
-    {
-      size <- diam/20
-    } else { size <- 1}
-  }
-  # apply species-specific colors
-  if (!is.null(species))
-  {
-    # load palette
-    color <- speciesColor()
-    # extract corresponding colors
-    col1 <- color[as.character(species), c("col","abvr")]
-  } else { col1 <- data.frame(col="black")}
-  if (add==FALSE)
-  {
-    #grDevices::dev.new(height=2.5,width=2.5)
-    graphics::par(mar=c(2.5,2.5,0.5,0.5))
-    graphics::plot(xy[,1],xy[,2],cex=size,col=col1$col,xlab="NA",ylab="NA",yaxt="n",xaxt="n",asp=1)
-    graphics::axis(2,mgp=c(0.5,0.5,0))
-    graphics::axis(1,mgp=c(0.5,0.5,0))
-    graphics::mtext(side=2,text="Easting (m)",line=1.3)
-    graphics:: mtext(side=1,text="Northing (m)",line=1.3)
-  } else { 
-    graphics::points(xy[,1],xy[,2],cex=size,col=col1$col)
-  }
-  # add color legend
-  if (!is.null(species))
-  {
-    texte <- sort(unique(col1$abvr))
-    graphics::legend("topleft", texte, col=color[texte,"col"],pch=15,cex=1,y.intersp = 1)
-  }
-}
-##################################
-#' Raster mask by union of buffers around xy positions
-#' 
-#' creates a raster mask by union of circular buffers around xy positions
-#'
-#' @param xy 2 columns matrix or data.frame. xy positions
-#' @param buff vector.  buffers to apply to the xy positions
-#' @param r raster object. target raster
-#' @param binary boolean. should the output mask be boolean (TRUE) or greyscale (FALSE)
-#' @return a raster object
-#' @seealso \code{\link{rasterChullMask}}
-#' @examples
-#' # create raster
-#' r <- raster::raster()
-#' raster::extent(r) <- c(0,40,0,40)
-#' raster::res(r) <- 1
-#' 
-#' # xy positions
-#' xy <- data.frame(c(10,20,31.25,15),
-#'   c(10,20,31.25,25))
-#' # compute mask
-#' mask1 <- rasterXYMask(xy, c(5, 8, 5, 5), r)
-#' mask2 <- rasterXYMask(xy, c(5, 8, 5, 5), r, binary=FALSE)
-#' 
-#' \dontrun{
-#' # display binary raster
-#' raster::plot(mask1)
-#' graphics::points(xy)
-#' 
-#' # display distance raster
-#' raster::plot(mask2)
-#' graphics::points(xy)}
-#' @export
-rasterXYMask <- function(xy, buff, r, binary=TRUE)
-{
-  # compute squared buffers
-  buff2 <- buff^2
-  # compute XY coordinates of cell centers
-  dummyXY <- raster::xyFromCell(r,1:length(r))
-  # create matrix of cell values
-  val <- matrix(NA,nrow=length(r), ncol=nrow(xy))
-  # compute cell value for each position
-  for (i in 1:nrow(xy))
-  {
-    val[,i] <- sqrt(pmax(0,buff2[i] - ((dummyXY[,1]-xy[i,1])^2 + (dummyXY[,2]-xy[i,2])^2)))
-  }
-  # take maximum for each cell
-  val <- apply(val,1,max)
-  # convert to binary if required
-  if (binary) {val <- val >0}
-  raster::values(r) <- val
-  r
-}
-####################################
-#' Raster mask of convex hull
-#' 
-#' creates raster mask corresponding to the convex hull of xy positions
-#'
-#' @param xy 2 columns matrix or data.frame. xy positions
-#' @param r raster object. target raster
-#' @return a raster with 0 or 1
-#' @examples
-#' # create raster
-#' r <- raster::raster()
-#' raster::extent(r) <- c(0,40,0,40)
-#' raster::res(r) <- 1
-#' 
-#' # xy positions
-#' xy <- data.frame(c(10,20,31.25,15),
-#'   c(10,20,31.25,25))