Area calculation problem

#
# Title: Spatial subsetting reproducable error creation
# Created: June 11th, 2019
# Last Updated: June 11th, 2019
# Author: Brandon Allen
# Objectives: Create a script that reproduces the spatial subsetting error you found in the cure4insect package
# Keywords: 
#

# Initialize the environment
library(cure4insect)

rm(list=ls())
gc()

opar <- set_options(path = "S:/reports", version = "2018")
getOption("cure4insect") # Confirm options

load_common_data() 

# Select species and subset the data to the North Saskatchewan (NSR.results) and the Lower Athabasca (LAR.result)

subset_common_data(id = get_all_id(luf = "North Saskatchewan"), species = "AlderFlycatcher")

# Error is thrown
# Error in base::colSums(x, na.rm = na.rm, dims = dims, ...) : 
#   'x' must be an array of at least two dimensions

NSR.results <- flatten(calculate_results(load_species_data("AlderFlycatcher")))
rownames(NSR.results) <- "NorthSaskatchewan"

clear_subset_data()

subset_common_data(id = get_all_id(luf = "Lower Athabasca"), species = "AlderFlycatcher")

# Error is thrown
# Error in base::colSums(x, na.rm = na.rm, dims = dims, ...) : 
#   'x' must be an array of at least two dimensions

LAR.results <- flatten(calculate_results(load_species_data("AlderFlycatcher")))
rownames(LAR.results) <- "LowerAthabasca"

clear_subset_data()

# Bind results into single data frame for visual inspection
final.results <- rbind(NSR.results, LAR.results)
rm(NSR.results, LAR.results)

final.results

# The values based only on abundances differ between the subset areas as expected.
# The Area_X columns both have values of 0, which in turn affect the Unit_X columns.

# When the kgrid_areas_by_sector R workspace is loaded, it contains a large dgCMatrix called KA_2016
# However, in the subset_common_data() function, it calls either KA_2012 or KA_2014 depending on the taxa.
# This loading happens on lines 123 and 125 in this version of the function below.

subset_common_data <-
        function(id=NULL, species="all")
        {
                if (!is_loaded())
                        stop("common data needed: use load_common_data")

                if (!is.null(dim(species))) # if provided as table, use 1st col
                        species <- as.character(species[,1L])
                vals <- c("all","birds","lichens","mammals","mites","mosses","vplants",
                          "upland", "lowland", "native", "nonnat", "north", "south")
                x <- .c4if$SP
                if (length(species) == 1L && tolower(species)  %in% vals) {
                        if (species %in%  c("all", "birds","lichens","mammals","mites","mosses","vplants"))
                                SPPfull <- get_all_species(taxon=species)
                        if (species %in% c("upland", "lowland"))
                                SPPfull <- get_all_species(habitat=species)
                        if (species %in% c("native", "nonnat"))
                                SPPfull <- get_all_species(status=species)
                        if (species %in% c("north", "south"))
                                SPPfull <- get_all_species(mregion=species)
                } else {
                        SPPfull <- species
                        if (any(SPPfull %ni% rownames(x)))
                                stop("all species must be valid IDs")
                }
                SPPfull <- rownames(x)[rownames(x) %in% SPPfull]
                if (length(SPPfull) <= 0)
                        stop("no species selected")

                if (is.null(id))
                        id <- rownames(.c4if$KT)
                if (inherits(id, "SpatialPolygons"))
                        id <- overlay_polygon(id)
                if (!is.null(dim(id))) # if provided as table, use 1st col
                        id <- as.character(id[,1L])
                if (!is.character(id))
                        id <- as.character(id)
                if (length(id) <= 0)
                        stop("no spatial IDs selected")
                ## QS-to-km mapping
                if (.validate_id(id, type="qs")) {
                        if (.verbose()) {
                                cat("matching quarter sections\n")
                                flush.console()
                        }
                        id <- qs2km(id[id %in% get_all_qsid()])
                }
                ## validating Row_Col IDs
                id <- id[id %in% rownames(.c4if$KT)]
                id <- sort(id)
                if (!.validate_id(id, type="km"))
                        stop("spatial id not valid")
                id10 <- sort(unique(as.character(.c4if$KT[id, "Row10_Col10"])))
                KT <- .c4if$KT
                ## South: >= 0; North: <= 0
                KT$mregion <- as.integer(-1*.select_id("north") + .select_id("south"))

                ## assignment happens only if all checks passed
                if (.verbose()) {
                        cat("arranging subsets\n")
                        flush.console()
                }
                clear_subset_data()
                assign("KTsub", KT[id,,drop=FALSE],
                       envir=.c4is)
                assign("A_2012", Matrix::colSums(.c4if$KA_2012[id,,drop=FALSE]),
                       envir=.c4is)
                assign("A_2014", Matrix::colSums(.c4if$KA_2014[id,,drop=FALSE]),
                       envir=.c4is)
                assign("SPfull", x,
                       envir=.c4is)
                assign("SPsub", x[SPPfull,,drop=FALSE],
                       envir=.c4is)
                invisible(NULL)
        }
ABbiodiversity / cure4insect

Area calculation problem #36
