Nothing
################################################################################
#### optimiser
################################################################################
#'Function for optimising a landscape genetic analysis based on resistance layers
#'
#'@param landscape resistance layer as a raster object (if not projected dimensions are assumed to be in meters!!)
#'@param nlocations the number of locations
#'@param mindist minimal distance in map units (meter if not specified)
#'@param fixed n fixed locations, provided as a data.frame with dimenstions nx2
#'@param method least cost path algorithm provided by the gdistance package. Options are "leastcost", "SPDistance" and "commute". see function costdistances.
#'@param NN number of next neighboring cells to be used, default is 8. see function costdistances.
#'@param iter number of iterations that should be used to find an optimised design. Try initially the default and if it runs fast (depends on the type of costdistance an d dimenstions of landscape), increase to higher values.
#'@param retry number of retries if optimisation is not possible in a certain iteration (due to mindist and/or fixed locations). Defaults to 10, which should be sufficient in most cases.
#'@param mask a raster object that masks positions which are not available. Areas which are not to be used for locations are coded as NA (not available), all other values are treated as suitable for placing locations.
#'@param plot if true, some plots are presented that show the histogramm, ecdf of the best (and the worst scenario).
#'@return a list object with two slots. The first slot is called opt and contains iter optimisation values (sd.detour and sd.cost) in a iter x 2 dimenstional data.frame. The second slot is called scenario and contains the corrosponding 1:iter scenarios, given as coordinates in a data.frame of dimensions nlocations x 2. Both slots are ordered in decreasing order of sd.detour values. So the best scenario is at 1 and the worst is at position iter.
#'@description opt.landgen places iter times nlocations in a given landscape. For each scenario the pairwise costdistances and Euclidean distances are calculated and the standard deviation of detour (see Gruber et al. in prep) is calculated. This metric evaluates the scenerio in their ability to find a significant effect for the given resistance layer on population structure (based on the causal modelling approach of Cushman et al.). To allow for more realistic designs previously locations, a minimal distance between locations and a mask that indicates "forbidden" areas can be specified.
#'@examples
#'\donttest{
#' library(raster)
#' fric.raster <- readRDS(system.file("extdata","fric.raster.rdata", package="PopGenReport"))
#'opt.landgen(landscape = fric.raster, nlocations = 5, mindist = 3,
#'method = "leastcost", NN = 8, iter = 10)
#'}
#' @export
opt.landgen <- function(landscape, nlocations, mindist=0, fixed = NULL, method, NN=8, iter =100, retry=10, mask=NULL, plot=TRUE)
{
scenario <- list()
opt <- data.frame(sd.cost=NA, sd.detour=NA)#, sd.detour2=NA)
if (!is.null(fixed))
{
specified <- nrow(fixed)
} else specified = 0
nadd <- nlocations - specified
if (is.na(crs(landscape))) crs(landscape) <-"+proj=merc +units=m" #already proojected data ????
for (it in 1:iter)
{
r1 <-landscape
values(r1) <- NA
#mask #NA values are cutted, everything else is left
if (!is.null(mask))
{
r1 <- mask
values(r1)<- ifelse(is.na(values(r1)),1,NA)
}
retryc <- retry
##Define x and y locations
# random no mindist specified
if (mindist==0) {
r1inv <- r1
values(r1inv) <- ifelse(is.na(values(r1inv)),1,NA)
rp<-coordinates(as(r1inv,"SpatialPointsDataFrame") )
choosep <- sample(1:nrow(rp),nadd)
xs <- rp[choosep,1]
ys <- rp[choosep,2]
}
# random with minddist>0
if (mindist>0)
{
if (!is.null(fixed))
{
rd <- rasterize(fixed, r1,1)
rd <- buffer(rd, mindist)
r1 <- sum(r1,rd, na.rm=TRUE)
values(r1) <- ifelse(values(r1)>0,1,NA)
}
left <- sum(is.na(values(r1)))
if (left==0) {message("There is no area left to place the required number of locations after placing the fixed locations. Reduce mindist or the amount of fixed locations.\n"); stop()}
xc <- NA
yc <- NA
i <- 1
rback <- r1
while (i<=nadd)
{
#choose point from left over areas
r1inv <- r1
values(r1inv) <- ifelse(is.na(values(r1inv)),1,NA)
rp<-coordinates(as(r1inv,"SpatialPointsDataFrame") )
choosep <- sample(1:nrow(rp),1)
xs <- rp[choosep,1]
ys <- rp[choosep,2]
xc[i] <- xs
yc[i] <- ys
i <- i+1
#add new point to mask
rd <- rasterize(cbind(xs,ys), r1,1)
rd <- buffer(rd, mindist)
r1 <- sum(r1,rd, na.rm=TRUE)
values(r1) <- ifelse(values(r1)>0,1,NA)
left <- sum(is.na(values(r1)))
if (left==0 & i<=nadd)
{
retryc <- retryc - 1
message(paste("No area left after ",i,"points at iteration",it,". I go back and try again.", retryc, "retries left...\n"))
i <- 1
r1 <- rback
if (retryc<1) {stop("Could not find any good combination, reduce mindist or increase retry or reduce number of fixed locations.\n")}
}
}
xs <- xc
ys <- yc
}
locs <- cbind(xs,ys)
rownames(locs) <- LETTERS[1:nadd]
#put fixed and locs together
locs <- rbind(fixed, locs)
#create a costdistance matrix
cost.mat <- costdistances(landscape, locs, method, NN)
eucl.mat <- as.matrix(dist(locs))
sd.cost <- sd(as.dist(cost.mat))
sd.detour = sd(resid(lm(as.numeric(as.dist(cost.mat))~as.numeric(as.dist( eucl.mat)))))
#sd.detour2 <- sd(as.dist(cost.mat)-as.dist(eucl.mat))
opt[it,] <-c(sd.cost, sd.detour)#, sd.detour2)
scenario[[it]] <- locs
} #end of iter loop
if (plot)
{
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow=c(2,2), mai=c(0.5,0.5,0.2,0.2))
opt.val <- which.max(opt$sd.detour)
locs.opt <- scenario[[opt.val]]
locs <- locs.opt
cost.mat <- costdistances(landscape, locs, method, NN)
eucl.mat <- as.matrix(dist(locs))
detour = resid(lm(as.numeric(as.dist(cost.mat))~as.numeric(as.dist( eucl.mat))))
hist(detour, main="Distrubtion of detour", col="darkgreen")
plot(ecdf(opt$sd.detour), main="Cummulative Density of sd.detour")
abline(v=opt[opt.val,"sd.detour"] ,col="blue")
#best versus worst case...
opt.val <- which.max(opt$sd.detour)
locs.opt <- scenario[[opt.val]]
locs <- locs.opt
plot(landscape, main=paste("best:",round(opt[opt.val,"sd.detour"],2) ))
points(locs[,1], locs[,2], pch=16, cex=1.2, col=c(rep("blue", specified), rep("black", nadd)))
text(locs[,1],locs[,2], row.names(locs), cex=1)
opt.val <- which.min(opt$sd.detour)
locs.opt <- scenario[[opt.val]]
locs <- locs.opt
plot(landscape, main=paste("worst:",round(opt[opt.val,"sd.detour"],2) ))
points(locs[,1], locs[,2], pch=16, cex=1.2, col=c(rep("blue", specified), rep("black", nadd)))
text(locs[,1],locs[,2], row.names(locs), cex=1)
par(oldpar)
}
ord <- order(opt$sd.detour, decreasing = TRUE)
return(list(opt=opt[ord,], scenario = scenario[ord]))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.