Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(raster)
- set.seed(17) # Seed a reproducible random sequence
- nr <- 30 # Number of rows
- nc <- 50 # Number of columns
- #
- # Create a zone raster for normalizing the probabilities.
- #
- zone <- raster(ncol=nc, nrow=nr)
- zone[] <- 0
- #
- # Create a probability raster (for illustrating the algorithm later).
- #
- p <- raster(ncol=nc, nrow=nr)
- p[] <- (1:(nc*nr) - 1/2) / (nc*nr) + rnorm(nc*nr, sd=0.5)
- p <- abs(focal(p, ngb=5, run=mean))
- z <- zonal(p, zone, stat='sum')
- p <- p / z[[2]] # This normalizes p to sum to unity as required
- #------------------------------------------------------------------------------#
- #
- # The algorithm begins here.
- #
- pvec <- sort(getValues(p), decreasing=TRUE) # The probabilities, sorted
- d <- cumsum(pvec) # Cumulative probabilities
- dpos <- d[d <= 0.95] # Position to stop
- region <- p # Initialize the output
- region[p < pvec[length(dpos)]] <- NA # Exclude the last 5% of the probability
- plot(region) # Display the result
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement