Search code examples
rmapslatitude-longitudergdaldismo

Species distribution model to produce psuedo-absence data using the randomPoints() function in the Dismo package in R: Error Message


Aim:

My aim is to build a species distribution prediction model using the function randomPoints() in the Dismo package with the utlimate aim of generating pseudo-absence points and plotting them on a map. These points will be converted into a rastor file in order to extract meta-data (i.e sea surface salinity, chlorophyll levels etc) from MODIS files to act as important ecological predictors to determine how they affect the distribution of blue whales.

The idea is to plug both presence and pseudo absence data with associated meta-data values into general linear mixed (GLM’s), which will ultimately make my models more balanced and accurate.

Problem Outline:

I am attempting to follow this species distribution exercise to generate pseudo-absence points (desired output: see image 1) using the randomPoints() function. However, after running my R-code (see below), I am experiencing this R error message (see below). My r-code also produces a map with the GPS points plotted (image 2).

Error Message:

Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘nlayers’ for signature ‘"standardGeneric"’

I have tried to find the solution. However, I am relatively new in regards to at working with maps in R, and I am feeling really confused as to what the problem is here!

My data frame contains 918 rows and I would like to produce the same number of pseudo-absence points as presence points. Unfortunately, I cannot publish my data but I have provided a mini data frame to use as an example.

If anyone can help me, I would be deeply appreciative!

Many thanks in advance!

R-code:

 ###Open Packages

    library("sp")
    library("rgdal")
    library("raster")
    library("maptools")
    library("rgdal")
    library("dismo")
    library("spatialEco")
    library("ggplot2")
    library("dplyr")
    library("maps")
    library("ggspatial")
    library("GADMTools")
    library("maps")

    ##Mini Dataframe
    Blue.whale_New <- data.frame(longitude = c(80.5, 80.5, 80.5, 80.5, 80.4, 80.4, 80.5, 80.5, 80.4),
                             latitude = c(5.84, 5.82, 5.85, 5.85, 5.89, 5.82, 5.82, 5.84, 5.83))


    ####World Bioclim Data + GADM Object

    ##Creating an array object with just longitude and latitude decimal coordinates
    ##Upload the maps
    ##Plotting the map of Sri Lanka
    ###GADM OBJECT

    dev.new()
    bioclim1.data <- getData('GADM', country='LKA', level=1)

    #####Worldclim raster layers

    bioclim.data <- getData(name = "worldclim",
                            var = "bio",
                            res = 2.5,
                            path = "./")

    ####Get bounding box of Sri Lanka shape file
    bb=bioclim1.data@bbox

    # Determine geographic extent of our data
    max.lat <- ceiling(max(Blue.whale$latitude))
    min.lat <- floor(min(Blue.whale$latitude))
    max.lon <- ceiling(max(Blue.whale$longitude))
    min.lon <- floor(min(Blue.whale$longitude))
    geographic.extent <- extent(x = c(min.lon, max.lon, min.lat, max.lat))

    #####Plot map
    dev.new()

    plot(bioclim1.data, 
         xlim = c(min(c(min.lon,bb[1,1])), max(c(max.lon,bb[1,2]))),
         ylim = c(min(c(min.lat,bb[2,1])), max(c(max.lat,bb[2,2]))),
         axes = TRUE, 
         col = "grey95")

    # Add the points for individual observation
    points(x = Blue.whale$longitude, 
           y = Blue.whale$latitude, 
           col = "olivedrab", 
           pch = 15, 
           cex = 0.50)

    ###Building a model and visualising results

    ##Crop bioclim data to geographic extent of blue whales GADM Map
    bioclim.data.blue.whale_1<-crop(x=bioclim1.data, y=geographic.extent)

    ##Crop bioclim data to geographic extent of blue whales World Clim
    bioclim.data.blue.whale_2<-crop(x = bioclim.data, y = geographic.extent)

    #Build distribution mode using the World Clim)

    bw.model <- bioclim(x = bioclim.data, p = Blue.whale_New)

    # Predict presence from model
    predict.presence <- dismo::predict(object = bw.model, x = bioclim.data, ext = geographic.extent)

    # Plot base map

    dev.new()

    plot(bioclim1.data, 
         xlim = c(min(c(min.lon,bb[1,1])), max(c(max.lon,bb[1,2]))),
         ylim = c(min(c(min.lat,bb[2,1])), max(c(max.lat,bb[2,2]))),
         axes = TRUE, 
         col = "grey95")

    # Add model probabilities
    plot(predict.presence, add = TRUE)

    # Redraw those country borders
    plot(bioclim1.data, add = TRUE, border = "grey5")

    # Add original observations
    points(Blue.whale_New$longitude, 
           Blue.whale_New$latitude, 
           col = "olivedrab", pch = 20, cex = 0.75)

    ##Psuedo Absence Points

    # Use the bioclim data files for sampling resolution
    bil.files <- list.files(path = "data/wc2-5", 
                            pattern = "*.bil$", 
                            full.names = TRUE)

    # Randomly sample points (same number as our observed points)

    ##Mask = provides resolution of sampling points
    ##n = number of random points
    ##ext = Spatially restricts sampling
    ##extf = expands sampling a little bit

    background <- randomPoints(mask = mask,    
                               n = nrow(Blue.whale_New),     
                               ext = geographic.extent, 
                               extf = 1.25)  

#Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘nlayers’ for signature ‘"standardGeneric"’  

Code I plan to run once I have solved my problem

Visualize the pseudo-absence points on a map:

# Plot the base map

dev.new()

plot(bioclim1.data, 
     xlim = c(min(c(min.lon,bb[1,1])), max(c(max.lon,bb[1,2]))),
     ylim = c(min(c(min.lat,bb[2,1])), max(c(max.lat,bb[2,2]))),
     axes = TRUE, 
     col = "grey95")

# Add the background points
points(background, col = "grey30", pch = 1, cex = 0.75)

# Add the observations
points(x = Blue.whale_New$longitude, 
       y = Blue.whale_New$latitude, 
       col = "olivedrab", 
       pch = 20, 
       cex = 0.75)

box()

# Arbitrarily assign group 1 as the testing data group
testing.group <- 1

# Create vector of group memberships
group.presence <- kfold(x = Blue.whale_New, k = 5) # kfold is in dismo package

Image 1 (Desired Output)

enter image description here

Image 2:

enter image description here


Solution

  • I think your problem is that you are not passing a proper mask (e.g. a raster layer) to the randomPoints function, but rather you are passing the mask function itself, which is a standardGeneric, hence the error message.

    You can generate a mask from the predictor raster you want to use and then pass it to the randomPoints function

    mask_r <- bioclim1.data[[1]] #just one raterLayer
    mask_r[!is.na(mask_r)] <- 1 #set all non-NA values to 1
    background <- randomPoints(mask = mask_r, #note the proper mask layer    
                               n = nrow(Blue.whale_New),     
                               ext = geographic.extent, 
                               extf = 1.25)
    

    I haven't checked, but this solution should work.

    Best,

    Emilio