La continuité du service gitlab.irstea.fr sera assurée en 2020. Nous envisageons ensuite une évolution vers une forge nationale INRAE encore à construire. Nous vous tiendrons au courant des évolutions futures.

Commit f6b06957 authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu

list of raster replaced by rasterStack in outputs

treeExtraction output is a spatial data.frame
updated raster support inside treeSegmentation
parent eb35fca5
......@@ -7,11 +7,11 @@
#'
#' Compute statistics by aggregating a raster at lower resolution. Aggregation groups are larger cells, new values are computed by applying a user-specified function to original cells contained in the larger cells. Results are provided as a data.frame which also contains the XY coordinates of the larger cells.
#'
#' @param r raster object or data.frame with xy columns
#' @param r raster object or data.frame with xy coordinates in two first columns
#' @param res numeric. Resolution of the aggregation raster, should be a multiple of r resolution if a raster is provided
#' @param fun function. Function to compute metrics in each aggregated cell from the values contained in the initial raster (use x$layer to access raster values) / data.frame (use x$colum_name to access values)
#' @param output string. indicates the class of output object "raster" or "data.frame"
#' @return a data.frame with the XY center coordinates of the aggregated cells, and the values computed with the user-specified function
#' @return a data.frame with the XY center coordinates of the aggregated cells, and the values computed with the user-specified function or a Raster object
#' @examples
#' data(laschablais3)
#'
......@@ -40,15 +40,21 @@
#' @export
rasterMetrics <- function(r, res=20, fun=function(x){data.frame(mean=mean(x$layer), sd=stats::sd(x$layer))}, output="raster")
{
if (is.element(class(r), c("RasterLayer", "RasterBrick")))
if (is.element(class(r), c("RasterLayer", "RasterBrick", "RasterStack")))
{
# convert to data.frame
st <- as.data.frame(raster::rasterToPoints(r))
projinfo <- r@crs
} else {
st <- r
projinfo <- NA
if (class(r) == "SpatialPointsDataFrame")
{
st <- cbind(r@coords, r@data)
projinfo <- r@proj4string
} else {
st <- r
projinfo <- NA
}
}
# compute coordinates of new cell center at metrics resolution
st$X <- round((st[,1]-res/2)/res)*res+res/2
st$Y <- round((st[,2]-res/2)/res)*res+res/2
......
......@@ -63,7 +63,7 @@ createDisk <- function(width=5)
#' # plot image after median and value-dependent Gaussian filters
#' plot(im2$smoothed.image, main="Value-dependent smoothing")}
#' @seealso \code{\link{maximaDetection}}, filters of imager package: \code{\link[imager]{mclosing}}, \code{\link[imager]{medianblur}}, \code{\link[imager]{deriche}}
#' @return A list of two cimg or rasterLayer objects: image after non-linear filter and image after both filters
#' @return A list of two cimg or a RasterStack objects: image after non-linear filter and image after both filters
#' @export
demFiltering <- function(dem, nlFilter="Closing", nlSize=5, sigmap=0.3)
{
......@@ -119,21 +119,23 @@ demFiltering <- function(dem, nlFilter="Closing", nlSize=5, sigmap=0.3)
# convert cimg objects to rasterLayer if necessary
if (class(dem)[1] =="RasterLayer")
{
dem.nl <- cimg2Raster(dem.nl, dem)
dem.gs <- cimg2Raster(dem.gs, dem)
output <- raster::addLayer(cimg2Raster(dem.nl, dem), cimg2Raster(dem.gs, dem))
names(output) <- c("non.linear.image", "smoothed.image")
} else {
output <- list(non.linear.image=dem.nl, smoothed.image=dem.gs)
}
list(non.linear.image=dem.nl, smoothed.image=dem.gs)
output
}
###################################
#' Local maxima extraction on image
#'
#' Variable window size maxima detection is performed on the image to extract local maxima position and calculate the window size where they are global maxima. Gaussian white noise is added to the image to avoid adjacent maxima due to neighbor pixels with identical value.
#'
#' @param dem cimg object (e.g. as created by \code{\link[imager]{cimg}}) or Raster Layer object (e.g. obtained with \code{\link[raster]{raster}})
#' @param dem cimg object (e.g. as created by \code{\link[imager]{cimg}}) or RasterLayer object (e.g. obtained with \code{\link[raster]{raster}})
#' @param dem.res numeric. image resolution, in case \code{dem} is a rasterLayer object, \code{dem.res} is extracted from the object by \code{\link[raster]{res}}
#' @param max.width numeric. maximum kernel width in pixel to check for local maximum
#' @param jitter boolean. indicates if noise should be added to image values to avoid the adjacent maxima due to the adjacent pixels with equal values
#' @return A cimg object or rasterLayer opbject which values are the radius (n) in meter of the square window (width 2n+1) where the center pixel is global maximum
#' @return A cimg object or RasterLayer object which values are the radius (n) in meter of the square window (width 2n+1) where the center pixel is global maximum
#' @examples
#' data(chmchablais3)
#'
......@@ -196,7 +198,7 @@ maximaDetection <- function(dem, dem.res=1, max.width=21, jitter=TRUE)
#' \item values in the maxima image (corresponding to the radius of the neighborhood where they are global maxima) are below a threshold depending on the initial image value.
#' }
#'
#' @param maxi cimg object or rasterLayer object. image with local maxima (typically output from \code{\link{maximaDetection}}, image values correspond to neighborhood radius on which pixels are global maxima in the initial image)
#' @param maxi cimg object or RasterLayer object. image with local maxima (typically output from \code{\link{maximaDetection}}, image values correspond to neighborhood radius on which pixels are global maxima in the initial image)
#' @param dem.nl cimg object. initial image from which maxima were detected
#' @param hmin numeric. minimum value in initial image for a maximum to be selected
#' @param dmin numeric. intercept term for selection of maxima depending on neighborhood radius: \code{maxi >= dmin + dem.nl * dprop}
......@@ -258,7 +260,7 @@ maximaSelection <- function(maxi,dem.nl,hmin=0,dmin=0,dprop=0)
#'
#' @param maxi cimg or rasterLayer object. image with seed points (e.g. from \code{\link{maximaDetection}} or \code{\link{maximaSelection}} )
#' @param dem.nl cimg or rasterLayer object. image for seed propagation (typically initial image used for maxima detection).
#' @return A list of 2 cimg or rasterLayer objects (image of segments id, image of maximum value in segment)
#' @return A list of 2 cimg object or RasterStack object with image of segments id and image of maximum value in segment)
#' @examples
#' data(chmchablais3)
#'
......@@ -320,10 +322,13 @@ segmentation <- function(maxi,dem.nl)
# convert to raster if necessary
if (israster)
{
dem.w <- cimg2Raster(dem.w, dem)
dem.wh <- cimg2Raster(dem.wh, dem)
output <- raster::addLayer(cimg2Raster(dem.w, dem), cimg2Raster(dem.wh, dem))
names(output) <- c("id", "maxvalue")
output
} else {
output <- list(id=dem.w, maxvalue=dem.wh)
}
list(id=dem.w, maxvalue=dem.wh)
output
}
################################
#' Modification of segments based on values
......@@ -422,15 +427,15 @@ segAdjust <- function(dem.w, dem.wh, dem.nl, prop=0.3, min.value=2, min.maxvalue
#' plot(segments$local.maxima, main="Local maxima")
#' #
#' # replace segment with id 0 (not a tree) with NA
#' segments$segments[segments$segments==0] <- NA
#' plot(segments$segments %% 8 , main="Segments")
#' segments$segments.id[segments$segments.id==0] <- NA
#' plot(segments$segments.id %% 8 , main="Segments")
#' #
#' # plot segmentation with other parameters
#' segments2$segments[segments2$segments==0] <- NA
#' plot(segments2$segments %% 8, main="Segments2")}
#' segments2$segments.id[segments2$segments.id==0] <- NA
#' plot(segments2$segments.id %% 8, main="Segments2")}
#'
#' @seealso \code{\link{demFiltering}}, \code{\link{maximaDetection}}, \code{\link{maximaDetection}}, \code{\link{maximaSelection}}, \code{\link{segmentation}}, \code{\link{segAdjust}}, \code{\link{treeExtraction}}
#' @return A list of 4 rasters: selected local maxima (values = distance to higher pixel), segments, non-linear preprocessed dem, smoothed preprocessed dem
#' @return A RasterStack with 4 layers: selected local maxima (values = distance to higher pixel), segments, non-linear preprocessed dem, smoothed preprocessed dem
#' @export
treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0,dprop=0.05,hmin=5,crownProp=0.3,crownMinH=2)
{
......@@ -439,8 +444,6 @@ treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0
if (is.character(dem)) {dem <- raster::raster(dem)}
# convert NAs to 0s
dem[is.na(dem)] <- 0
# convert to cimg
dem.c <- imager::as.cimg(raster::as.matrix(dem))
#
# conversion of Gaussian filter standard deviation from meters to pixels
if (length(sigma==1))
......@@ -450,9 +453,9 @@ treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0
sigma[,2] <- sigma[,2]/raster::res(dem)[1]
}
# dem filtering
temp <- demFiltering(dem.c, nlFilter=nlFilter, nlSize=nlSize, sigmap=sigma)
dem.nl <- temp[[1]]
dem.gs <- temp[[2]]
temp <- demFiltering(dem, nlFilter=nlFilter, nlSize=nlSize, sigmap=sigma)
dem.nl <- temp$non.linear.image
dem.gs <- temp$smoothed.image
#
# maxima detection
maxi <- maximaDetection(dem.gs,raster::res(dem)[1])
......@@ -462,22 +465,17 @@ treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0
#
# segmentation
temp <- segmentation(maxi,dem.nl)
dem.w <- temp[[1]]
dem.wh <- temp[[2]]
dem.w <- temp$id
dem.wh <- temp$maxvalue
#
# segmentation adjustment
dem.w <- segAdjust(dem.w,dem.wh,dem.nl,crownProp,crownMinH, hmin) # removal of trees with top heigth < hmin
# remove trees from r.maxi
maxi[dem.w==0] <- 0
#
# conversion to raster
r.maxi <- raster::raster(as.matrix(maxi))
r.dem.nl <- raster::raster(as.matrix(dem.nl))
r.dem.gs <- raster::raster(as.matrix(dem.gs))
r.dem.w <- raster::raster(as.matrix(dem.w))
raster::extent(r.maxi) <- raster::extent(r.dem.w) <- raster::extent(r.dem.nl) <- raster::extent(r.dem.gs) <- raster::extent(dem)
#
list(local.maxima=r.maxi, segments.id=r.dem.w, filled.dem=r.dem.nl, smoothed.dem=r.dem.gs)
output <- raster::addLayer(maxi, dem.w, dem.nl, dem.gs)
names(output) <- c("local.maxima", "segments.id", "filled.dem", "smoothed.dem")
output
}
################################
#' Tree extraction
......@@ -488,7 +486,7 @@ treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0
#' @param r.maxi raster object. raster positive values at local maxima
#' @param r.dem.w raster object. segmented raster
#' @param r.mask raster object. only segments which maxima are inside the mask are extracted
#' @return A data.frame with tree id, local maximum stats (height, coordinates, dominance radius), surface and volume.
#' @return A spatial data.frame with tree id, local maximum stats (height, dominance radius), segment stats (surface and volume).
#' @examples
#' data(chmchablais3)
#'
......@@ -497,7 +495,7 @@ treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0
#'
#' # tree extraction
#' trees <- treeExtraction(segments$filled.dem, segments$local.maxima, segments$segments.id)
#' summary(trees)
#' trees
#'
#' \dontrun{
#' # plot initial image
......@@ -508,7 +506,7 @@ treeSegmentation <- function(dem, nlFilter="Closing", nlSize=5, sigma=0.3,dmin=0
#' plot(contours, add=T, border="white")
#'
#' # add treetop positions
#' points(trees$x, trees$y, cex=trees$h/20)}
#' plot(trees, cex=trees$h/20, add=T, pch=1)}
#'
#' @seealso \code{\link{treeSegmentation}}
#' @export
......@@ -552,6 +550,8 @@ treeExtraction <- function(r.dem.nl, r.maxi, r.dem.w, r.mask=NULL)
segms <- merge(segms, sp, all.x = T)
segms <- merge(segms, vp, all.x = T)
}
sp::coordinates(segms) <- ~ x+y
segms@proj4string <- r.dem.nl@crs
segms
}
#
......
......@@ -16,7 +16,7 @@ demFiltering(dem, nlFilter = "Closing", nlSize = 5, sigmap = 0.3)
\item{sigmap}{numeric or matrix. if a single number is provided, sigmap is the standard deviation of the Gaussian filter in pixel, 0 corresponds to no smoothing. In case of matrix, the first column corresponds to the standard deviation of the filter, and the second to thresholds for image values (e.g. a filter of standard deviation specified in line \code{i} is applied to pixels in image which values are between thresholds indicated in lines \code{i} and \code{i+1}). Threshold values should be ordered in increasing order.}
}
\value{
A list of two cimg or rasterLayer objects: image after non-linear filter and image after both filters
A list of two cimg or a RasterStack objects: image after non-linear filter and image after both filters
}
\description{
applies two filters to an image:
......
......@@ -7,7 +7,7 @@
maximaDetection(dem, dem.res = 1, max.width = 21, jitter = TRUE)
}
\arguments{
\item{dem}{cimg object (e.g. as created by \code{\link[imager]{cimg}}) or Raster Layer object (e.g. obtained with \code{\link[raster]{raster}})}
\item{dem}{cimg object (e.g. as created by \code{\link[imager]{cimg}}) or RasterLayer object (e.g. obtained with \code{\link[raster]{raster}})}
\item{dem.res}{numeric. image resolution, in case \code{dem} is a rasterLayer object, \code{dem.res} is extracted from the object by \code{\link[raster]{res}}}
......@@ -16,7 +16,7 @@ maximaDetection(dem, dem.res = 1, max.width = 21, jitter = TRUE)
\item{jitter}{boolean. indicates if noise should be added to image values to avoid the adjacent maxima due to the adjacent pixels with equal values}
}
\value{
A cimg object or rasterLayer opbject which values are the radius (n) in meter of the square window (width 2n+1) where the center pixel is global maximum
A cimg object or RasterLayer object which values are the radius (n) in meter of the square window (width 2n+1) where the center pixel is global maximum
}
\description{
Variable window size maxima detection is performed on the image to extract local maxima position and calculate the window size where they are global maxima. Gaussian white noise is added to the image to avoid adjacent maxima due to neighbor pixels with identical value.
......
......@@ -7,7 +7,7 @@
maximaSelection(maxi, dem.nl, hmin = 0, dmin = 0, dprop = 0)
}
\arguments{
\item{maxi}{cimg object or rasterLayer object. image with local maxima (typically output from \code{\link{maximaDetection}}, image values correspond to neighborhood radius on which pixels are global maxima in the initial image)}
\item{maxi}{cimg object or RasterLayer object. image with local maxima (typically output from \code{\link{maximaDetection}}, image values correspond to neighborhood radius on which pixels are global maxima in the initial image)}
\item{dem.nl}{cimg object. initial image from which maxima were detected}
......
......@@ -8,7 +8,7 @@ rasterMetrics(r, res = 20, fun = function(x) { data.frame(mean =
mean(x$layer), sd = stats::sd(x$layer)) }, output = "raster")
}
\arguments{
\item{r}{raster object or data.frame with xy columns}
\item{r}{raster object or data.frame with xy coordinates in two first columns}
\item{res}{numeric. Resolution of the aggregation raster, should be a multiple of r resolution if a raster is provided}
......@@ -17,7 +17,7 @@ rasterMetrics(r, res = 20, fun = function(x) { data.frame(mean =
\item{output}{string. indicates the class of output object "raster" or "data.frame"}
}
\value{
a data.frame with the XY center coordinates of the aggregated cells, and the values computed with the user-specified function
a data.frame with the XY center coordinates of the aggregated cells, and the values computed with the user-specified function or a Raster object
}
\description{
Compute statistics by aggregating a raster at lower resolution. Aggregation groups are larger cells, new values are computed by applying a user-specified function to original cells contained in the larger cells. Results are provided as a data.frame which also contains the XY coordinates of the larger cells.
......
......@@ -12,7 +12,7 @@ segmentation(maxi, dem.nl)
\item{dem.nl}{cimg or rasterLayer object. image for seed propagation (typically initial image used for maxima detection).}
}
\value{
A list of 2 cimg or rasterLayer objects (image of segments id, image of maximum value in segment)
A list of 2 cimg object or RasterStack object with image of segments id and image of maximum value in segment)
}
\description{
performs a seed-based watershed segmentation and outputs two images: one with segment id and one where the segments are filled with the highest value of the segment in the initial image
......
......@@ -16,7 +16,7 @@ treeExtraction(r.dem.nl, r.maxi, r.dem.w, r.mask = NULL)
\item{r.mask}{raster object. only segments which maxima are inside the mask are extracted}
}
\value{
A data.frame with tree id, local maximum stats (height, coordinates, dominance radius), surface and volume.
A spatial data.frame with tree id, local maximum stats (height, dominance radius), segment stats (surface and volume).
}
\description{
creates a dataframe with segment id, height and coordinates of maxima, surface and volume, computed from three images: initial, local maxima and segmented, typically obtained with \code{\link{treeSegmentation}}
......@@ -29,7 +29,7 @@ segments <- treeSegmentation(chmchablais3)
# tree extraction
trees <- treeExtraction(segments$filled.dem, segments$local.maxima, segments$segments.id)
summary(trees)
trees
\dontrun{
# plot initial image
......@@ -40,7 +40,7 @@ contours <- raster::rasterToPolygons(segments$segments.id, dissolve=T)
plot(contours, add=T, border="white")
# add treetop positions
points(trees$x, trees$y, cex=trees$h/20)}
plot(trees, cex=trees$h/20, add=T, pch=1)}
}
\seealso{
......
......@@ -28,7 +28,7 @@ treeSegmentation(dem, nlFilter = "Closing", nlSize = 5, sigma = 0.3,
\item{crownMinH}{numeric. minimum crown height}
}
\value{
A list of 4 rasters: selected local maxima (values = distance to higher pixel), segments, non-linear preprocessed dem, smoothed preprocessed dem
A RasterStack with 4 layers: selected local maxima (values = distance to higher pixel), segments, non-linear preprocessed dem, smoothed preprocessed dem
}
\description{
global function for preprocessing (filtering), maxima detection and selection, segmentation and segmentation adjustment of a raster image.
......@@ -48,12 +48,12 @@ plot(segments$smoothed.dem, main="Filtered image")
plot(segments$local.maxima, main="Local maxima")
#
# replace segment with id 0 (not a tree) with NA
segments$segments[segments$segments==0] <- NA
plot(segments$segments \%\% 8 , main="Segments")
segments$segments.id[segments$segments.id==0] <- NA
plot(segments$segments.id \%\% 8 , main="Segments")
#
# plot segmentation with other parameters
segments2$segments[segments2$segments==0] <- NA
plot(segments2$segments \%\% 8, main="Segments2")}
segments2$segments.id[segments2$segments.id==0] <- NA
plot(segments2$segments.id \%\% 8, main="Segments2")}
}
\references{
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment