R/addscenario.R

In remi-daigle/BESTMPA:

#' Add a protection scenario
#'
#' @description Adds protection scenarios to the SpatialPolygonsDataFrame grid. Variables in the dataframe of the SpatialPolygonsDataFrame that should otherwise be logical (e.g. scenarios) are converted to numeric since logical variables are incompatible with ESRI shapefiles.
#'
#' @param domain the gridded model domain (SpatialPolygonsDataFrame) as generated by inigrid
#' @param priority A logical vector indicating which grid cells are to be prioritized when assigning MPA protection. The algorithm will select all priority cells regardless of distance/spacing before selecting non-priority cells.
#' @param excluded A logical vector indicating which grid cells are to be excluded from MPA assignment
#' @param included A logical vector indicating which grid cells are to be included in MPA assignment. This accomodates the pre-existing MPAs.
#' @param MPA_coverage Target protection level as a proportion (e.g. 0.2 is 20\% protection)
#' @param replicates A vector of replicates for the main body of the model
#' @param dist The fixed distance for setting MPA distance in m in fixed distance scenarios.
#' @param name The name of the scenario to be added. While any character string will be accepted, using the format "MPA_XX" where XX is any 2 letters which indicates treatment type is suggested to minimize conflict with other functions in the BESTMPA package.
#' @param cell_size cell size of the grid used in initgrid
#' @param cells minimum number of cells that can be considered as an MPA. defaults to 1
#'
#' @return SpatialPolygonsDataFrame
#' @export
#'
#' @examples
#' proj  <- "+proj=utm +zone=20 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
#' p <- initgrid(EEZ=EEZ,cell_size=20000,proj=proj,areas=c(Habitats=Habitats,Breeding=Breeding))
#' # set targeted scenario
#' p <- addscenario(domain=p,included=oldMPA,priority=p$Habitats,excluded=p$Breeding,MPA_coverage=0.1,replicates=1:3,dist=75000,name="MPA_TR",cell_size=20000,cells=3)
addscenario <- function(domain,included=NA,priority=NA,excluded=NA,MPA_coverage=0.1,replicates=25,dist=NA,name="MPA_SQ",cell_size,cells=1){
    if(name=="MPA_SQ"){

        for(i in replicates){
            domain[[paste0(name,"_",i)]] <- as.numeric(included)
        }

        return(domain)

    } else {
        # create neighbours and graphs

        require(igraph)
        require(spdep)

        nb <- poly2nb(domain,snap=100)
        gr <- spdf2graph(domain,nb)

        # dist from m to cells
        if(!is.na(dist)) dist <- round(dist/cell_size)+1

        # numeric to logical
        if(!is.logical(included)) included <- as.logical(included)
        if(!is.logical(priority)) priority <- as.logical(priority)
        if(!is.logical(excluded)) excluded <- as.logical(excluded)


        # loop through replicates

        for(i in replicates){

            print(paste("calculating",name, "replicate",i))

            # set status_quo as MPAs
            domain[[paste0(name,"_",i)]] <- as.numeric(included)

            # generate new MPA sizes
            MPA_sizes <- generatempasizes(cell_size,cells,domain,MPA_coverage,included)

            # set distances
            if(is.na(dist)) dist <- round(sqrt(length(domain))/(sqrt(length(MPA_sizes))-1))


            # plot(domain)
            prioritize <- TRUE
            for(j in MPA_sizes){

                candidates <- NULL
                tol <- 1

                while(length(candidates)<1){
                    candidates <- selectmpa(domain,gr,dist,priority,prioritize,excluded,name,tol,i)
                    tol <- tol*0.9
                }
                # while(length(candidates)<1){
                #
                #     candidates <- unique(unlist(neighborhood(gr,order=dist,node=which(as.logical(domain@data[,names(domain)==paste0(name,"_",i)])))))
                #
                #     # remove non-priority and excluded
                #     if(!any(is.na(priority))&prioritize) candidates <- candidates[priority[candidates]]
                #     if(!any(is.na(excluded))) candidates <- candidates[!excluded[candidates]]
                #
                #     candidates <- candidates[!candidates %in% unlist(neighborhood(gr,order=dist*tol,node=which(as.logical(domain@data[,names(domain)==paste0(name,"_",i)]))))]
                #
                #     tol <- tol*0.9
                # }

                seed <- sample(candidates,1)
                while(length(seed)<j){

                    seed <- unique(c(seed,sample(unlist(nb[seed]))))
                    # # remove non-priority and excluded
                    if(!any(is.na(priority))&prioritize) seed <- seed[priority[seed]]
                    if(!any(is.na(excluded))) seed <- seed[!excluded[seed]]
                    if(length(seed)==1){
                        print(paste("MPA cannot increase in size, selecting different candidate"))
                        candidates <- candidates[!candidates %in% seed]
                        if(length(candidates)==0){
                            tol <- tol*0.9
                            candidates <- selectmpa(domain,gr,dist,priority,prioritize,excluded,name,tol,i)
                        }
                        seed <- sample(candidates,1)
                        count <- 0
                    }
                }

                if(length(seed)>j) seed <- seed[1:j]
                # plot(domain[seed,],col='blue',add=T)
                domain[[paste0(name,"_",i)]][seed] <- 1
                if(sum(domain[[paste0(name,"_",i)]])>=sum(priority,na.rm=T)) prioritize <- FALSE

            }

        }
        return(domain)

    }

}