Fetching all records for a taxonomic group in OBIS

For example, if you want to download the records for all Nudibranchia then in R it is as simple as doing:

    install.packages("robis")

    data <- robis::occurrence("Nudibranchia")

    # get species data
    spdata <- data[!is.na(data$taxonRank) & data$taxonRank == "species",]

    # plot on a map
    robis::leafletmap(spdata)

Alternatively, you can download all data with the new mapper that is being developed:

http://iobis.org/mapper2?scientificname=Nudibranchia&download=true

Error in linbin2D: (list) object cannot be coerced to type 'double'

When calling bkde2D you might encounter the following error:

Error in linbin2D(x, gpoints1, gpoints2) : 
  (list) object cannot be coerced to type 'double'

This issue can manifest itself when the passed in data is a tibble instead of a data.frame.

Easy fix:

data <- as.data.frame(data)

Speeding up mean and standard deviation calculation for a stack of rasters in R

While modelling the distribution of invasive seaweeds during my PhD I wanted to calculate the mean and standard deviation of a set of SDMs that appeared to be equivalent in for the current climate but who gave different results for the future climate. However, the standard naive approach for calculating the mean and standard deviation of a RasterStack was rather slow. I managed to get a 20 fold speed up by calculating the mean and standard deviation by aggregating them in a stream. This rather obvious for calculating the mean, summing all rasters and then dividing by the number of rasters. For the standard deviation, the formula is only slightly more complicated and can be found both on Wikipedia and in code on StackOverflow. If your raster contains very large numbers, you might need to adapt this code to use Welford's method.

# Option 1: shorter but slower

rasterstack_meansd_slow <- function(x) {
mean <- raster::mean(x)
  sd <- raster::calc(x, sd)

  list(mean=mean, sd=sd)

}

# Option 2: faster but more code
rasterstack_meansd_fast <- function(x) {
  s0 <- nlayers(x)
  s1 <- raster(x, layer=1)
  s2 <- s1^2
  for(ri in 2:s0) {
    r <- raster(x, layer=ri)
    s1 <- s1 + r
    s2 <- s2 + r^2
  }
  list(mean=s1/s0, sd=sqrt((s0 * s2 - s1 * s1)/(s0 * (s0 - 1))))
}

As a small example I calculate the mean and standard deviation of the sea surface temperature for different climate scenarios for 2050. The fast version only takes 10 seconds on my machine while the slow version takes 225 seconds to calculate the mean and standard deviation.

library(sdmpredictors)

sstfuture <- load_layers(c('BO2_RCP26_2050_tempmean_ss', 'BO2_RCP45_2050_tempmean_ss', 
                           'BO2_RCP60_2050_tempmean_ss', 'BO2_RCP85_2050_tempmean_ss'))

system.time({ fast <- rasterstack_meansd_fast(sstfuture) }) # 10 seconds
system.time({ slow <- rasterstack_meansd_slow(sstfuture) }) # 225 seconds

Creating a single pixel wide raster along the coastline or another edge

At one point during my PhD I needed a single pixel of raster values along the coastline. In order to get this I used the following code in R:

library(sdmpredictors)
library(raster)

coast_mask <- function(layer) {
  edges <- raster::boundaries(raster(layer, layer=1), type="inner")
  values <- getValues(edges)
  is.na(values) | values == 0
}

l <- load_layers("BO_sstmean", equalarea = FALSE)

mask <- coast_mask(l)
l[mask] <- NA

plot(l, col=rev(heat.colors(255)))

Source code available on GitHub at https://github.com/samuelbosch/blogbits/blob/master/misc/coastalraster.R

Main R Packages used:

• sdmpredictors
• raster

Creating a high resolution png from a raster in R

In this post I'll show you two different ways to create a png from an equal-area bathymetry grid in R.

Version 1: using the raster::plot function

The first option for plotting a raster is to use the raster::plot function and write this to a sufficiently large png file. This works reasonably well but the only disadvantage is that both the title and legend are really small. The title is optional but I didn't find a way to disable the legend. The main trick to get the raster::plot function to output a high resolution png consists in setting the maxpixels parameter. Colors where selected on ColorBrewer but with a small addition that the color range was reverted and that the darkest color is repeated to ensure that rare depths (-6000 to -10000) get the same color.

if(!requireNamespace("sdmpredictors")) {
  install.packages("sdmpredictors")
}
library(raster)
x <- sdmpredictors::load_layers("BO_bathymean", equalarea = TRUE)
png("bathymetry_plot1.png", width=ncol(x), height=nrow(x))
col <- rev(c("#f7fbff", "#deebf7", "#c6dbef", "#9ecae1", "#6baed6", "#4292c6","#2171b5","#08519c", rep("#08306b",7)))
plot(x, maxpixels = ncell(x), col = col, colNA = "#818181", 

     main = "Bathymetry", axes = FALSE, ylim=extent(x)[3:4])
dev.off()

Version 2: write to png using leaflet colors

While this version looks and is a bit more complicated, it produces really good looking results. The main gist of this code is that it transforms the values to a range from 0 to 1000 while making sure that extreme values at both extremes will get the same colors as more common values. This is similar to what you general would do when e.g. creating a color scale with QGIS. Once values are mapped to colors they are converted to raw bytes with the right dimensions and written to a png file. Remark that this code was inspired by some of the internal functions in the leaflet package.

if(!requireNamespace("leaflet")) {
  install.packages("leaflet")
}
if(!requireNamespace("sdmpredictors")) {
  install.packages("sdmpredictors")
}
library(sdmpredictors)
library(raster)

# create colors
colors <- leaflet::colorNumeric(rev(c("#f7fbff", "#deebf7", "#c6dbef", "#9ecae1", "#6baed6", "#4292c6","#2171b5","#08519c", "#08306b")),
                                -1:1001, na.color =  "#818181")
cols <- c(colors(-1:1001), colors(NA))
x <- sdmpredictors::load_layers("BO_bathymean", equalarea = TRUE)
# scale values and remove extreme values from the color range 
vals <- values(x)
vals <- scale(vals)
minmax <- quantile(vals, probs=c(0.01, 0.99), na.rm = TRUE)
vals <- round((((vals - minmax[1]) / (minmax[2] - minmax[1])) * 1000))
vals[vals < 0] <- 0
vals[vals > 1000] <- 1000
vals[is.na(vals)] <- 1002

# lookup colors for scaled values, convert to raw and write to file
valcolors <- cols[vals+2] # +2 because -1 and 0 are in cols (value 0 is at index 2 in cols)

rgb_data <- col2rgb(valcolors, alpha = TRUE)
raw_data <- as.raw(rgb_data)
dim(raw_data) <- c(4, ncol(x), nrow(x))

png::writePNG(raw_data, "bathymetry_plot2.png")

All source code is available on GitHub at https://github.com/samuelbosch/blogbits/blob/master/misc/raster2png.R

Main R Packages used:

• sdmpredictors
• raster
• leaflet

dismo pwdSample and ssb for a large number of species distribution records

For my MarineSPEED dataset of species I'm trying to create non-random cross-validation datasets. In order to do this I was looking at the pwdSample and ssb functions from the dismo package. pwdSample calculates a pair-wise distance sampling and ssb calculates the spatial sorting bias. More information about these functions can be found in the excellent paper by Robert Hijmans (2012): "Cross-validation of species distribution models: removing spatial sorting bias and calibration with a null model".

Problem is that both functions calculate a distance matrix between two sets of points which starts to fail when the number of points is (very) large. To fix this I re-implemented the functions replacing the matrix calculation and subsequent row minimum finding with the same functionality but using a divide and conquer strategy by calculating distance matrices on subsets, taking the min distance and then taking the min distances of these min distances to get row based minima.

In other words I replaced:

fromd <- apply(distfun(fixed, reference), 1, min)

with:

mindist <- function(distfun, a, b) {
  partition_count <- (1 %/% (100 / NROW(b))) + 1
  parts <- dismo::kfold(x=b, k=partition_count)
  r <- c()
  for(i in 1:partition_count) {
    mind <- apply(distfun(a, b[parts==i,]), 1, min)
    r <- cbind(r, mind)
  }
  apply(r, 1, min)
}
  
fromd <- mindist(distfun, fixed, reference)

The full pwdSample function then becomes:

pwdSample_robust <- function(fixed, sample, reference, tr = 0.33, nearest= TRUE, n=1, lonlat = TRUE, warn = TRUE) {
  ## DIVIDE reference, calculate distance matrix, get min, COMBINE: take min of the min, do rest of the pwdSample logic
  distHaversine <- function(p1, p2) {
    r <- 6378137
    toRad <- pi/180
    p1 <- p1 * toRad
    p2 <- p2 * toRad
    p <- cbind(p1[, 1], p1[, 2], p2[, 1], p2[, 2])
    dLat <- (p[, 4] - p[, 2])
    dLon <- (p[, 3] - p[, 1])
    a <- sin(dLat/2) * sin(dLat/2) + cos(p[, 2]) * cos(p[, 
                                                         4]) * sin(dLon/2) * sin(dLon/2)
    dist <- 2 * atan2(sqrt(a), sqrt(1 - a)) * r
    as.vector(dist)
  }
  distGeo <- function(x, y) {
    n <- nrow(x)
    m <- nrow(y)
    dm <- matrix(ncol = m, nrow = n)
    for (i in 1:n) {
      dm[i, ] <- distHaversine(x[i, , drop = FALSE], y)
    }
    return(dm)
  }
  distPlane <- function(x, y) {
    dfun <- function(x, y) {
      sqrt((x[, 1] - y[, 1])^2 + (x[, 2] - y[, 2])^2)
    }
    n = nrow(x)
    m = nrow(y)
    dm = matrix(ncol = m, nrow = n)
    for (i in 1:n) {
      dm[i, ] = dfun(x[i, , drop = FALSE], y)
    }
    return(dm)
  }
  if (lonlat) {
    distfun <- distGeo
  }
  else {
    distfun <- distPlane
  }
  stopifnot(tr > 0)
  n <- round(n)
  stopifnot(n >= 1)
  if (inherits(fixed, "SpatialPoints")) 
    fixed <- coordinates(fixed)
  if (inherits(sample, "SpatialPoints")) 
    sample <- coordinates(sample)
  if (inherits(reference, "SpatialPoints")) 
    reference <- coordinates(reference)
  fixed <- as.matrix(fixed)[, 1:2]
  sample <- as.matrix(sample)[, 1:2]
  reference <- as.matrix(reference)[, 1:2]
  if (warn) {
    if (nrow(sample) < nrow(fixed)) {
      warning("nrow(sample) < nrow(fixed)")
    }
  }
  
  mindist <- function(distfun, a, b) {
    partition_count <- (1 %/% (100 / NROW(b))) + 1
    parts <- dismo::kfold(x=b, k=partition_count)
    r <- c()
    for(i in 1:partition_count) {
      mind <- apply(distfun(a, b[parts==i,]), 1, min)
      r <- cbind(r, mind)
    }
    apply(r, 1, min)
  }
  
  fromd <- mindist(distfun, fixed, reference) ##apply(distfun(fixed, reference), 1, min)
  tod <- mindist(distfun, sample, reference) ##apply(distfun(sample, reference), 1, min)
  ngb <- matrix(NA, nrow = length(fromd), ncol = n)
  iter <- sample(1:nrow(fixed))
  for (j in 1:n) {
    for (i in iter) {
      d <- abs(tod - fromd[i])
      if (min(d) < (tr * fromd[i])) {
        if (nearest) {
          x <- which.min(d)
        }
        else {
          x <- sample(which(d < (tr * fromd[i])), size = 1)
        }
        ngb[i, j] <- x
        tod[x] <- Inf
      }
    }
  }
  return(ngb)
}

and the ssb function then looks like this:

ssb <- function(p, a, reference, lonlat = TRUE, avg = TRUE) {
  distHaversine <- function(p1, p2) {
    r <- 6378137
    toRad <- pi/180
    p1 <- p1 * toRad
    p2 <- p2 * toRad
    p <- cbind(p1[, 1], p1[, 2], p2[, 1], p2[, 2])
    dLat <- (p[, 4] - p[, 2])
    dLon <- (p[, 3] - p[, 1])
    a <- sin(dLat/2) * sin(dLat/2) + cos(p[, 2]) * cos(p[, 
                                                         4]) * sin(dLon/2) * sin(dLon/2)
    dist <- 2 * atan2(sqrt(a), sqrt(1 - a)) * r
    as.vector(dist)
  }
  distGeo <- function(x, y) {
    n <- nrow(x)
    m <- nrow(y)
    dm <- matrix(ncol = m, nrow = n)
    for (i in 1:n) {
      dm[i, ] <- distHaversine(x[i, , drop = FALSE], y)
    }
    return(dm)
  }
  distPlane <- function(x, y) {
    dfun <- function(x, y) {
      sqrt((x[, 1] - y[, 1])^2 + (x[, 2] - y[, 2])^2)
    }
    n = nrow(x)
    m = nrow(y)
    dm = matrix(ncol = m, nrow = n)
    for (i in 1:n) {
      dm[i, ] = dfun(x[i, , drop = FALSE], y)
    }
    return(dm)
  }
  if (lonlat) {
    distfun <- distGeo
  }
  else {
    distfun <- distPlane
  }
  if (inherits(p, "SpatialPoints")) 
    p <- coordinates(p)
  if (inherits(a, "SpatialPoints")) 
    a <- coordinates(a)
  if (inherits(reference, "SpatialPoints")) 
    reference <- coordinates(reference)
  p <- as.matrix(p)[, 1:2]
  a <- as.matrix(a)[, 1:2]
  reference <- as.matrix(reference)[, 1:2]
  
  mindist <- function(distfun, a, b) {
    partition_count <- (1 %/% (1000 / NROW(b))) + 1
    parts <- dismo::kfold(x=b, k=partition_count)
    r <- c()
    for(i in 1:partition_count) {
      mind <- apply(distfun(a, b[parts==i,]), 1, min)
      r <- cbind(r, mind)
    }
    apply(r, 1, min)
  }
  pd <- mindist(distfun, p, reference) #   pd <- apply(distfun(p, reference), 1, min)
  ad <- mindist(distfun, a, reference) #   ad <- apply(distfun(a, reference), 1, min)
  if (avg) {
    res <- cbind(mean(pd), mean(ad))
    colnames(res) <- c("p", "a")
    return(res)
  }
  else {
    return(list(pd, ad))
  }
}

Spatially partioning lat long data in R with sperrorest partition.kmeans

While trying out sperrorest::partition.kmeans I started wondering whether it did handle the date line (180, -180) correctly. A small clustering test confirmed that sperrorest doesn't take into account the spatial projection of the data and doesn't cluster together points from the same spatial area. To solve this I decided to convert my lat long data into cartesian coordinates and then call the partition.kmeans function.

library(sperrorest)

lonlat_xyz <- function(data, lonlat_cols) {
   rad_coordinates <- data[,lonlat_cols] * (pi/180)
   lon <- rad_coordinates[,lonlat_cols[1]]
   lat <- rad_coordinates[,lonlat_cols[2]]
   r <- cbind(cos(lat)*cos(lon), cos(lat)*sin(lon), sin(lat))
   colnames(r) <- c("x", "y", "z")
   return(r)
}

cartesian <- lonlat_xyz(data, lonlat) 
partitions <- sperrorest::partition.kmeans(cartesian, coords = c("x","y","z"))

Weighted Relative Neighborhood Graph in R based on cccd::rng

The R package cccd contains a nice implementation of the Relative Neighborhood Graph (rng) but in the current version 1.5 it returns a non-weighted igraph. But for one of my experiments I needed the weighted version so I've slightly changed the code to get an igraph with weights.

rng <- function (x = NULL, dx = NULL, r = 1, method = NULL, usedeldir = TRUE,
          open = TRUE, k = NA, algorithm = "cover_tree", weighted = TRUE) {
  if (is.na(k)) {
    if (is.null(dx)) {
      if (is.null(x))
        stop("One of x or dx must be given.")
      dx <- as.matrix(proxy::dist(x, method = method))
    }
    else {
      usedeldir <- FALSE
    }
    n <- nrow(dx)
    A <- matrix(0, nrow = n, ncol = n)
    if (is.vector(x))
      x <- matrix(x, ncol = 1)
    if (usedeldir && ncol(x) == 2) {
      del <- deldir::deldir(x[, 1], x[, 2])
      for (edge in 1:nrow(del$delsgs)) {
        i <- del$delsgs[edge, 5]
        j <- del$delsgs[edge, 6]
        d <- min(apply(cbind(dx[i, -c(i, j)], dx[j, -c(i, j)]), 1, max))
        rd <- r * dx[i, j]
        if ((open && rd < d) || rd <= d) {
          A[i, j] <- A[j, i] <- rd
        }
      }
    } else {
      diag(dx) <- Inf
      for (i in 1:n) {
        for (j in setdiff(1:n, i)) {
          d <- min(apply(cbind(dx[i, -c(i, j)], dx[j, -c(i, j)]), 1, max))
          rd <- r * dx[i, j]
          if ((open && rd < d) || rd <= d) {
            A[i, j] <- A[j, i] <- rd
          }
        }
      }
    }
    diag(A) <- 0
    out <- graph.adjacency(A, mode = "undirected", weighted = weighted)
  } else {
    if (is.null(x))
      stop("x must not be null")
    n <- nrow(x)
    k <- min(k, n - 1)
    dx <- get.knn(x, k = k, algorithm = algorithm)
    edges <- NULL
    weights <- NULL
    for (i in 1:n) {
      i.indices <- dx$nn.index[i, ]
      i.dists <- dx$nn.dist[i, ]
      for (j in 1:k) {
        rd <- r * i.dists[j]/2
        j.indices <- dx$nn.index[i.indices[j], ]
        j.dists <- dx$nn.dist[i.indices[j], ]
        rd <- r * i.dists[j]
        S <- setdiff(intersect(i.indices, j.indices),
                     c(i, i.indices[j]))
        if (length(S) > 0) {
          d <- Inf
          for (si in S) {
            a <- which(i.indices == si)
            b <- which(j.indices == si)
            d <- min(d, max(i.dists[a], j.dists[b]))
          }
          if ((open && rd < d) || rd <= d) {
            edges <- cbind(edges, c(i, i.indices[j]))
            weights <- cbind(weights, rd)
          }
        }
      }
    }
    g <- graph(edges, n = n, directed = FALSE)
    if( weighted ) {
      edge.attributes(g) <- list(weight=weights)
    }
    out <- simplify(g, edge.attr.comb = "first")
  }
  if (!is.null(x)) {
    out$layout <- x
  }
  out$r <- r
  out
}

Workaround object 'couleursIn' not found in R package eVenn version 2.2

When trying to define a venn diagram with custom colors in the eVenn diagram:

set.seed(42)
d <- data.frame(a = sample(rep(c(0,1), 20)), 
                b = sample(rep(c(0,1), 20)), 
                c = sample(rep(c(0,1), 20)))
evenn(matLists=as.matrix(d), display = TRUE, couleurs = c("#90BA6E","#956EAD","#9F5845"))

I encountered the following error:

Error: object 'couleursIn' not found

This error only occurs when you pass a vector with multiple colors to couleurs and you've set Solid = TRUE, which is the default value. First thing to do after getting the error is calling dev.off() as the error will have left a connection to the png file open. To fix the error you have to define a variable couleursIn with colors with alpha values if you want to create similar venn diagrams as the default ones but with different colors:

couleursIn <- c("#90BA6E80","#956EAD80","#9F584580")
evenn(matLists=as.matrix(d), display = TRUE, couleurs = c("#90BA6E","#956EAD","#9F5845"))

Note the 80 at the end of each color which stands for an alpha value of 128 which is 50%. This can be checked with col2rgb("#90BA6E80", alpha = TRUE).

Workaround missing support for as.formula in R package maxlike version 0.1-5

Today I noticed that the predict function for maxlike didn't support models build with formulas from "as.formula". The error I got was:

Error in predict.maxlikeFit(model, data.frame(p), ...) : 
  at least 1 covariate in the formula is not in rasters.

I send in a pull request but in the mean time you can workaround this problem with the following wrapper for the maxlike function:

## workaround missing support for as.formula in maxlike
maxlike <- function(formula, ...) {
  m <- maxlike::maxlike(formula, ...)
  m$call$formula <- formula
  return(m)
}

Workaround "n.trees" is missing in R package gbm version 2.1.1

Just short note to anyone encountering following error in the gbm package version 2.1.1:

 Error in paste("Using", n.trees, "trees...\n") : 
  argument "n.trees" is missing, with no default 

This issue will be fixed in the next version but in the mean time I'll post my workaround. If you add following function to your code then everything should work:

## work around bug in gbm 2.1.1
predict.gbm <- function (object, newdata, n.trees, type = "link", single.tree = FALSE, ...) {
  if (missing(n.trees)) {
    if (object$train.fraction < 1) {
      n.trees <- gbm.perf(object, method = "test", plot.it = FALSE)
    }
    else if (!is.null(object$cv.error)) {
      n.trees <- gbm.perf(object, method = "cv", plot.it = FALSE)
    }
    else {
      n.trees <- length(object$train.error)
    }
    cat(paste("Using", n.trees, "trees...\n"))
    gbm::predict.gbm(object, newdata, n.trees, type, single.tree, ...)
  }
}

Benchmarking reading binary values from a file with F#, Python, Julia, R, Go and OCaml

In one of my recent posts I showed some F# code for storing and reading integers as binary data. Since then I've created different versions of the reading part in Python, Julia, Go, OCaml and R to be able to compare the performance of the different languages on this simple task and to get a feeling on what it's like to program in Julia, Go and OCaml as I hadn't created any programs in these languages yet.

Below table shows the benchmark results for reading 10 times 10000 values and 10000 times 10 values with links to the source code files in the different programming languages:

Language	10 times 10000 values	10000 times 10 values
F#	6 seconds	20 seconds
Python	26 seconds	40 seconds
Julia	45 seconds	72 seconds
Go	8 seconds	25 seconds
OCaml	2.5 seconds	48 seconds
R	110 seconds	NA

The overall fastest version is the one written in F# but note that it's also the version I have tweaked the most. As I'm not very experienced in most of the languages so any performance tips are welcome. Note that I tried using memory mapped files in .NET and Python this improved performance when querying lots of values from the same file but also made it worse in other cases.

The implementation of the functionality is most of the times rather similar in the different languages. Some notable differences where:

• Julia apparently doesn't have a null value so I refrained from checking whether the read integer value was equal to the int32 minimum value (-2147483648).
• In Go converting the bytes to integers was faster with a custom function.
• I didn't find a function in the OCaml Core library to convert bytes to a 32-bit integer, but luckily I found one on Stack Overflow.

F#:

open System
open System.IO

let readValue (reader:BinaryReader) cellIndex = 
    reader.BaseStream.Seek(int64 (cellIndex*4), SeekOrigin.Begin) |> ignore
    match reader.ReadInt32() with
    | Int32.MinValue -> None
    | v -> Some(v)
        
let readValues indices fileName = 
    use reader = new BinaryReader(File.Open(fileName, FileMode.Open, FileAccess.Read, FileShare.Read))
    let values = Array.map (readValue reader) indices
    values

Python:

def read_values(filename, indices):
    # indices are sorted and unique
    values = []
    with open(filename, 'rb') as f:
        for index in indices:
            f.seek(index*4L, os.SEEK_SET)
            b = f.read(4)
            v = struct.unpack("@i", b)[0]
            if v == -2147483648:
                v = None
            values.append(v)
    return values

Julia:

function readvalue(stream, position)
    seek(stream, position)
    return read(stream, Int32)
end

function readvalues(filename::String, indices)
    stream = open(filename, "r")
    try
        return Int32[readvalue(stream, index*4) for index in indices]
    finally
        close(stream)
    end
end

Go:

import ("os")
func bytes2int(b []byte) int32 {
    v := int32(0)
    for  i := 0; i < 4; i++ {
        v = v | (int32(b[i]) << (uint(8*i)))
    }
    return v
}

func readValues(indices []int, filename string) []int32 {
    results := make([]int32, len(indices))
    b := make([]byte, 4)
    f,_ := os.Open(filename)
    for i, cellIndex := range indices {
        f.Seek(int64(cellIndex*4), os.SEEK_SET)
        f.Read(b)
        value := bytes2int(b) // around 10-20% faster then binary.Read
        if value != -2147483648 {
            results[i] = value
        } else {
            results[i] = 99999
        }
    }
    return results
}

OCaml:

let input_le_int32 inchannel = (* http://stackoverflow.com/a/6031286/477367 *)
  let res = ref 0l in
    for i = 0 to 3 do
      let byte = input_byte inchannel in
        res := Int32.logor !res (Int32.shift_left (Int32.of_int byte) (8*i))
    done;

    match !res with
      | -2147483648l -> None
      | v -> Some(v)

let readvalue inchannel index =
  seek_in inchannel (index*4);
  input_le_int32 inchannel

let readvalues (indices:int array) filename =
  let inchannel = open_in_bin filename in
    try
      let result = Array.map (readvalue inchannel) indices in
        close_in inchannel;
        result
    with e ->
      close_in_noerr inchannel;
      raise e

R:

read.values <- function(filename, indices) {
  conn <- file(filename, "rb")
  read.value <- function(index) {
    seek(conn, where=index*4)
    readBin(conn, integer(), size = 4, n = 1, endian = "little")
  }
  r <- sapply(indices,read.value)
  close(conn)
  r[r==-2147483648] <- NA
  r
}

Any suggestions for improving the implementation in one of the above programming languages ? Which language would you like to compare my results with ? Which other language(s) do you expect to be faster than my benchmark results ? Can you help me out with a version in your favorite language or in C, fortran, Common Lisp, Scheme, Clojure, Java or J ?

BIOMOD2: invalid term in model formula

While trying out BIOMOD2 for species distribution modeling I encountered the following error when trying to create a random forest (RF) model:

Error in terms.formula(formula, data = data) : 
  invalid term in model formula
Error in predict(model.bm, Data[, expl_var_names, drop = FALSE], on_0_1000 = TRUE) : 
  error in evaluating the argument 'object' in selecting a method for function 'predict': Error: object 'model.bm' not found

After lots of digging in the source code of BIOMOD I discovered what happened and how to fix it.

The line where the error was occurring looked like this:

model.sp <- try(randomForest(formula = makeFormula(resp_name, 
    head(Data), "simple", 0), data = Data[calibLines, 
    ], ntree = Options@RF$ntree, importance = FALSE, 
    norm.votes = TRUE, strata = factor(c(0, 1)), 
    nodesize = Options@RF$nodesize, maxnodes = Options@RF$maxnodes))

As you can see a formula is created from our data. The generated formula looked like this:

654987634 ~ 1 + nitrate_mean + par_mean + sst_mean + sst_range

The number 654987634 is my response variable name that I set when calling BIOMOD_FormatingData which apparently has been added to the Data data.frame as the first column. 

Solution: adding a letter prefix to my response variable name fixed my issue.

Lessons learned:

1. It's incredibly useful that you can take a look at external source code in an RSession by typing its name, for iternals you add the package name and 3 double dots e.g. biomod2:::.Biomod.Models. Generics functions can be fetched with getMethod e.g. getMethod(biomod2::predict, signature="RF_biomod2_model")
2. Read the docs very carefully: this is the description for resp.name: response variable name (character). The species name.
3. An assertion in BIOMOD_FormatingData would have saved me an hour or two. But note that a naive check like is.character(resp.name) wouldn't work in this case, because once the formula is made it isn't a character anymore.
4. Adding more diagnostic information in your error messages saves people time. If the generated formula would have been printed as well then I would probably have noticed the issue directly.
5. Libraries make you win time, but you lose some as well.

If you liked this content then don't forget to subscribe on the top left.

Creating a kernel density estimate map in R

In this post I'm going to create a kernel density estimate map in R from a file with latitude/longitude coordinates. WARNING: depending on your application the following gives incorrect results because a non-spherical kernel density estimator is used with spherical data (big thanks too Brian Rowlingson for pointing that out). I didn't find any spherical kernel density estimation implementation in R but if you know any please let me know. While exploring this issue I also wrote a similar implementation as the below code for the generation of a kernel density estimate with the spatial statistics package splancs instead of KernSmooth (splancs_kernel_density.R).

To create this map we'll use the bkde2D function in the KernSmooth package to generate the kernel density estimates and the raster package for outputting a raster file. If you don't have the KernSmooth and raster packages then they can as usually be installed with:

install.packages("KernSmooth")
install.packages("raster")

In order to create the 2D binned kernel density map I first loaded the KernSmooth and raster package.

library("KernSmooth")
library("raster")

Then we read the input file with the coordinates of the points that we want to generate the kernel density estimate off. This input file is an export from the Ocean Biographic Information System (OBIS) and represents distribution records from the seaweed Alaria esculenta.

input <- "Alaria_esculenta.csv"
output <- "Alaria_esculenta_density.asc"

# get the coordinates
records <- read.csv(input)
coordinates <- records[,2:3]

Next we compute the 2D binned kernel density estimate. In order to limit the size of the generated output files I set all values smaller then 0.00001 to 0.

# compute the 2D binned kernel density estimate
est <- bkde2D(coordinates, 
              bandwidth=c(3,3), 
              gridsize=c(4320,2160),
              range.x=list(c(-180,180),c(-90,90)))
est$fhat[est$fhat<0.00001] <- 0 ## ignore very small values

Finally we create a raster file from the output of bkde2D, inspect it visually and export it as an ascii file.

# create raster
est.raster = raster(list(x=est$x1,y=est$x2,z=est$fhat))
projection(est.raster) <- CRS("+init=epsg:4326")
xmin(est.raster) <- -180
xmax(est.raster) <- 180
ymin(est.raster) <- -90
ymax(est.raster) <- 90
# visually inspect the raster output
plot(est.raster)
# write an ascii file
writeRaster(est.raster,output,"ascii")

And the output looks like this. As you can see most records are located near 50°N and 0° which corresponds to the distribution records map generated by OBIS.

The output ascii raster file can be used as an input file for OccurrenceThinner which uses the kernel density estimate grid for filtering out records from densely sampled regions in order to avoid overfitting while building species distribution models (SDM).

Code and data for this post are in my blog github repository.