## Warning: multiple methods tables found for 'spTransform'

In this example you will load in some zebra location data and determine how their ranges overlap.

The data originated from MoveBank

But I have cleaned it slightly for this workshop. Read the data file in - it is a shapefile of point data.

require(rgdal) zebras = readOGR("./Data", "zebra")

## OGR data source with driver: ESRI Shapefile ## Source: "./Data", layer: "zebra" ## with 34034 features and 3 fields ## Feature type: wkbPoint with 2 dimensions

```
summary(zebras)
```

## Object of class SpatialPointsDataFrame ## Coordinates: ## min max ## coords.x1 36.6111 37.0311 ## coords.x2 0.2257 0.7166 ## Is projected: FALSE ## proj4string : ## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0] ## Number of points: 34034 ## Data attributes: ## timestamp species individual ## 2007-06-17 06:00:00: 8 Equus burchellii: 5865 GZ4 :11233 ## 2007-06-17 07:00:00: 8 Equus grevyi :28169 GZ2 : 6155 ## 2007-06-17 08:00:00: 8 PZ2 : 5055 ## 2007-06-17 09:00:00: 8 GZ7 : 4198 ## 2007-06-17 10:00:00: 8 GZ3 : 2975 ## 2007-06-17 11:00:00: 8 GZ1 : 1773 ## (Other) :33986 (Other): 2645

```
table(zebras$individual)
```

## ## GZ1 GZ2 GZ3 GZ4 GZ5 GZ6 GZ7 PZ1 PZ2 ## 1773 6155 2975 11233 534 1301 4198 810 5055

This tells us we have records for two species and 9 individuals - 7 Grevy's Zebra, and two Plains (Burchelli) zebra. We also see it is in lat-long projection. Let's change to Google Mercator and get a background image.

zebras = spTransform(zebras, CRS("+init=epsg:3857")) summary(zebras)

## Object of class SpatialPointsDataFrame ## Coordinates: ## min max ## coords.x1 4075530 4122278 ## coords.x2 25125 79776 ## Is projected: TRUE ## proj4string : ## [+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 ## +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null ## +no_defs] ## Number of points: 34034 ## Data attributes: ## timestamp species individual ## 2007-06-17 06:00:00: 8 Equus burchellii: 5865 GZ4 :11233 ## 2007-06-17 07:00:00: 8 Equus grevyi :28169 GZ2 : 6155 ## 2007-06-17 08:00:00: 8 PZ2 : 5055 ## 2007-06-17 09:00:00: 8 GZ7 : 4198 ## 2007-06-17 10:00:00: 8 GZ3 : 2975 ## 2007-06-17 11:00:00: 8 GZ1 : 1773 ## (Other) :33986 (Other): 2645

require(dismo) bg = gmap(extent(zebras))

Now we can plot the zebra locations

plot(bg) plot(zebras, pch = 19, col = "#00000020", add = TRUE)

Now, lets choose nine colours and map them. The ColorBrewer system provides good colours for graphics, so we'll use 9 from one of its categorical colour palettes.

```
require(RColorBrewer)
```

palette(brewer.pal(9, "Set1")) plot(bg) plot(zebras, col = zebras$individual, add = TRUE) legend("left", fill = brewer.pal(9, "Set1"), legend = levels(zebras$individual), bg = "white", inset = 0.1)

First, lets try plotting the bounding box of each zebra.

- Split the SpatialPointsDataFrame into nine SPDFs based on the identifier
- Use a nice palette
- Plot the background map
- Add the zebra point
- Loop over the split list and plot the bounding boxes using the
`extent`

function from the`raster`

package

splitz = split(zebras, zebras$individual) palette(brewer.pal(9, "Set1")) plot(bg) plot(zebras, col = "black", cex = 0.5, add = TRUE) require(raster) for (i in 1:9) plot(extent(splitz[[i]]), add = TRUE, col = i, lwd = 5) title("Bounding box ranges")

Bounding boxes are a poor estimate of an animal's range. Lets try the convex hull. This is the smallest convex polygon that encloses a set of points. For example:

So let's compute the convex hull for each zebra and plot. The procedure is:

- The
`rgeos`

package for its`gConvexHull`

function. - A loop over each of the sets of points
- Save each convex hull in a list
- Loop over the hulls to plot them in colour.

palette(brewer.pal(9, "Set1")) require(rgeos) chulls = list() for (i in 1:9) { chulls[[i]] = gConvexHull(splitz[[i]]) } names(chulls) = names(splitz) plot(bg) plot(zebras, col = "black", cex = 0.5, add = TRUE) for (i in 1:9) { plot(chulls[[i]], border = i, lwd = 4, add = TRUE) }

Now some analysis. I want to see how much of each zebra's range overlaps with
the other zebras. I can compute the area of a spatial object using `gArea`

and
compute intersections of spatial objects using `gIntersection`

, both from the `rgeos`

package.

The output is going to be a 9x9 matrix. Each row is the percentage of that zebra's range that overlaps with the zebra in that column. Obviously the diagonal will be 100% since every range is an exact overlap with itself.

# empty matrix to fill zm = matrix(0, 9, 9) # loop over this zebra for (i in 1:9) { # compute area of its range areaI = gArea(chulls[[i]]) # loop over other zebra for (j in 1:9) { # find intersection region overlap = gIntersection(chulls[[i]], chulls[[j]]) # check if regions are distinct if (is.null(overlap)) { zm[i, j] = 0 } else { # compute percentage of this zebra's area zm[i, j] = round(100 * gArea(overlap)/areaI) } } } # make nice labels rownames(zm) = names(splitz) colnames(zm) = names(splitz) zm

## GZ1 GZ2 GZ3 GZ4 GZ5 GZ6 GZ7 PZ1 PZ2 ## GZ1 100 100 0 0 74 35 48 0 18 ## GZ2 28 100 0 0 25 12 15 0 7 ## GZ3 0 0 100 100 0 0 0 81 0 ## GZ4 0 0 20 100 0 0 0 34 0 ## GZ5 81 97 0 0 100 44 47 0 21 ## GZ6 78 99 0 0 91 100 45 0 48 ## GZ7 90 100 0 0 81 37 100 0 2 ## PZ1 0 0 25 52 0 0 0 100 0 ## PZ2 20 27 0 0 22 25 1 0 100

As an exercise, turn that code into a function. Think about the inputs and outputs.

What does this matrix tell us? The first zebra, GZ1, has a range that is enclosed totally within GZ2, but is completely different to GZ3 and GZ4.

The two plains zebras, PZ1 and PZ2, have non-overlapping ranges. We can plot these conclusions to check, and also label the ranges at the region centroid.

plot(bg) plot(chulls[["GZ1"]], add = TRUE, col = "#00FF0080") plot(chulls[["GZ2"]], add = TRUE, col = "#FF00FF80") plot(chulls[["PZ1"]], add = TRUE, col = "#FF000080") plot(chulls[["PZ2"]], add = TRUE, col = "#0000FF80") text(coordinates(chulls[["GZ1"]]), "GZ1") text(coordinates(chulls[["GZ2"]]), "GZ2") text(coordinates(chulls[["PZ2"]]), "PZ2") text(coordinates(chulls[["PZ1"]]), "PZ1")

We can combine regions using the `gUnion`

function. It takes two arguments and
returns a merged object. We have 9 regions to merge. One way to do that would be a
number of nested calls like this:

allRanges = gUnion(chulls[[1]], gUnion(chulls[[2]], gUnion(chulls[[3]], gUnion(chulls[[4]], gUnion(chulls[[5]], gUnion(chulls[[6]], gUnion(chulls[[7]], gUnion(chulls[[8]], chulls[[9]]))))))))

Luckily R has an easier way to do this, using the `Reduce`

function:

allRanges = Reduce(gUnion, chulls) plot(bg) plot(allRanges, add = TRUE)

This shows us that the zebras' ranges fall into two non-overlapping ranges

If we rearrange the matrix we can see the division into two communities more clearly:

zm[c(1, 2, 5, 6, 7, 9, 3, 4, 8), c(1, 2, 5, 6, 7, 9, 3, 4, 8)]

## GZ1 GZ2 GZ5 GZ6 GZ7 PZ2 GZ3 GZ4 PZ1 ## GZ1 100 100 74 35 48 18 0 0 0 ## GZ2 28 100 25 12 15 7 0 0 0 ## GZ5 81 97 100 44 47 21 0 0 0 ## GZ6 78 99 91 100 45 48 0 0 0 ## GZ7 90 100 81 37 100 2 0 0 0 ## PZ2 20 27 22 25 1 100 0 0 0 ## GZ3 0 0 0 0 0 0 100 100 81 ## GZ4 0 0 0 0 0 0 20 100 34 ## PZ1 0 0 0 0 0 0 25 52 100

The `$species`

field of the `zebras`

data frame
gives the species. Make some maps based on species.

You could also use alpha shapes, which are concave, but require a scale parameter, alpha. Also, if the input has coincident points, the algorithm can't cope. To work round that we add a tiny bit of 'jitter' to the coordinates.

Here's some alpha shapes for the first zebra:

require(alphahull) par(mfrow = c(1, 3)) plot(ashape(jitter(coordinates(splitz[[1]])), alpha = 10000), asp = 1) plot(ashape(jitter(coordinates(splitz[[1]])), alpha = 5000), asp = 1) plot(ashape(jitter(coordinates(splitz[[1]])), alpha = 1000), asp = 1)

par(mfrow = c(1, 1))

See what you get for some of the other zebras, and try and produce this map: