Updates were made to gDefrag

If you’re a user of the R based package gDefrag then you should update it. Version 0.2 will be available during the following days.

It brings no major features noticeable to the user. The need to develop a new version of the package was caused by changes to other packages, on which gDefrag depends.

Just out: “Structural simplification compromises the potential of common insectivorous bats to provide biocontrol services against the major olive pest Prays oleae”

Check out our new paper on biocontrol services in olive groves.

Here’s the abstract: 

Crop production intensification often leads to the structural simplification of production systems. This structural simplification is expected to have strong impacts on biodiversity and the provisioning of ecosystem services, but information about this topic is scarce. For instance, no information exists for Mediterranean olive (Olea europaea) groves, despite olive farming representing a significant share of the agricultural sector in some European countries. We investigated the impact of in-farm and landscape-level structural simplification on the potential of three common insectivorous bats (i.e., Pipistrellus kuhlii, P. pygmaeus and P. pipistrellus) to provide biocontrol services against one of the most harmful olive pests worldwide, the olive fruit moth Prays oleae. Bats and insect surveys were both carried out in olive groves representing increasing levels of structural simplification and during three sampling seasons (spring, summer and autumn). At grove-level, structural simplification was considered as resulting from reduced planting pattern variability (i.e., tree and row spacing) and tree features (diameter at breast height, height of the trunk and canopy area), while at landscape level was considered as resulting from reduced land-cover types. We found that the Kuhl’s pipistrelle was the most frequently recorded species in all types of olive groves and seasons. Moreover, the activity levels of pipistrelle bats as a whole significantly decreased with the structural simplification of olive groves. The abundance of P. oleae was highest at intermediate levels of structural simplification, irrespective of the season. Forest cover in the surrounding landscape had a significant positive influence on the activity levels of P. kuhlii, and a significant and negative influence on the abundance of P. oleae. Our study demonstrates that structural simplification differentially influences the activity patterns of both insectivorous bats and insect pests within olive groves. Moreover, it suggests that structural simplification may strongly compromise biocontrol services provided by bats on the major olive pest P. oleae.

 

Just accepted “Prioritizing road defragmentation using graph-based tools”

The paper I co-authored with Fernando Ascensão and Márcia Barbosa, entitled “Prioritizing road defragmentation using graph-based tools“, was just accepted for publication in Landscape and Urban Planning.

This has been a great collaboration, which already resulted in the publication of one R package (gDefrag) and one other paper!

Here’s the abstract:

Roads are a main cause of habitat fragmentation but mitigating the full road network is unfeasible. A key goal in the road mitigation planning process is to highlight, at the transportation network level, the most problematic roads, i.e. where mitigation measures are most required, in order to maximize the benefits for biodiversity while keeping implementation costs as low as possible. Grounded on the concepts of habitat amount and accessible habitat, we prioritized roads for mitigation based on dual spatial graphs, where the land polygons delimited by roads are the nodes and the roads themselves are the links. The rationale was to identify those links (roads) that connect the nodes with higher potential biodiversity (as a proxy for quality habitat). We applied this approach to prioritize the defragmentation of the major road network of the Iberian Peninsula, targeting all native mammalian carnivores inhabiting this region. Our goal was to identify those roads that, by dividing areas with the best habitat quality and/or are major potential barriers for connectivity, should be prioritized in the mitigation process. We used two complementary metrics: Area Weighted Metric and the Integral Index of Connectivity. Highlighted roads bisect regions of high potential biodiversity for carnivores in northern Spain and along the Portugal-Spain border. Thirty-five roads were scored as high-priority by both metrics, suggesting that they have particular impact both in the amount of quality habitat and in overall landscape functional connectivity. This approach is completely scalable, allowing a fast assessment from local to continental scales.

Updating MetaLandSim…

A new development version is know available for the package MetaLandSim (v. 1.0.6). It can be downloaded from GitHub.

The user has to run the following code in order to install the package from GitHub:

library(devtools)
install_github(“FMestre1/MetaLandSim”)

These changes were required because of recent improvements made to the package rgrass7 (which allows to work on GRASS from R). MetaLandSim relies on rgrass7 to properly run the function range_raster.

One of the outputs MetaLandSim provides is a simulation of range expansion by a given species (with a given set of characteristics) in a given landscape (with a given connectivity). This function, range_raster, converts the output of the simulation to rasters.

The final result are these two rasters (bellow), but check out the manual, and this paper describing the package to get further details. Also, if you want to see how these results can be used, and combined with Ecological Niche Models, check out this paper.

sample.png

I will update MetaLandSim on CRAN briefly.

Just out: “Bad moon rising? The influence of the lunar cycle on amphibian roadkills”

As a “side-project” I collaborated with people from my previous research group and wrote this paper on the influence on the lunar cycle on amphibian mortality caused by road traffic.

It provides a view on a different, but very relevant, environmental factor, the moon cycle. It might be helpful when designing conservation mesures to mitigate the negative effects of roads on wildlife.

I wish to thank Helena Lopes, Tiago Pinto, Luís Guilherme Sousa, António Mira and Sara Santos for this excellent cooperation.

You can read the full paper here, but here’s the abstract:

Annually, roads, and their associated users, are responsible for millions of roadkills worldwide. Mortality affects multiple taxonomic groups, but amphibians are particularly vulnerable, due to their size and underreporting. In fact, very high mortality frequencies can occur, mostly during short periods of time, when individuals migrate to and from reproduction areas (e.g., ponds). In this study, we assess the influence of the lunar cycle on amphibian roadkills, while accounting for weather conditions. As expected, the main environmental effects explaining roadkill numbers were weather related, with increases in minimum air temperature, average relative air humidity, and cumulative rainfall during the previous 24 h having a positive effect on roadkill numbers for all studied species. However, the lunar cycle also affected roadkills for two of the studied species. Darker nights had higher numbers of roadkills of Pleurodeles waltl, while moonlit nights had higher numbers of Salamandra salamandra. As such, these moon effects are species specific. Animals that are more active in moonlight may be at an advantage if their visual acuity is better than that of their predators. We hypothesize that differences between species in the response to moonlight may be due to differences perceived in predation risk. This information should be considered when designing mitigation measures. Volunteer actions, for instance, can be planned and coordinated keeping in mind the most appropriate weather conditions for the general amphibian community and specific phases of the lunar cycle for particular species.

I Iberian Meeting on Agroecological Research

Ibagreco2018

Last November took place at Évora University the I Iberian Meeting on Agroecological Research. It was an event for science presentation but also for discussion with those that could benefit from this initiative and research topics (scientists, farmers, managers, policy makers, etc).

If the topic interests you, feel free to download the book of abstracts or check out the related paper about the meeting.

A shinny app to compute and visualize roadkill hotspots over linear features

In this post I’m going to explain how to use the package roadHotspots. This is just a development version located in github and I must say that I have no further intentions of continuing to work on this. If someone find this useful and would like to continue where I’ve left, please contact me.

I’ve created this project just to demonstrate that with some (limited) programming knowledge, free open-source tools and approximately 20h of work it’s possible to develop a useful and simple software with a user-friendly interface and a great scientific potential.

To install the package you first need devtools. Then, just run the following code:

library(devtools)
devtools::install_github("bmsasilva/roadHotspots")

This package provides a shiny app that allows to compute and visualize hotspots over linear features (eg. roads) for a set of observations (eg. animal crossings, roadkills). I’ve just implemented two computation algorithms, Gaussian kernel density estimation and Malo’s method (Malo et al. 2004). As I said before this is just a proof of concept and it’s not really my area of expertise. I know that there are more advanced hotspot calculation algorithms but I’ll let that for the experts in this field. To start the app just do:

run_app()

pic1_app

Now that you have the app running you will see a series of buttons and drop down menus. On the top right corner you have two buttons to access the file system (indicated in the image by the gray boxes 1 and 2) to choose the data for the app to work. You must select 2 files: a shapefile with the roads you’re analyzing and a csv file with roadkills/observations. The csv must have three columns: species id, X_coord and Y_coord, by this order. The coordinate system of the shapefile is irrelevant, but the coordinate system used in the csv must be the same.

The drop down menus indicated by the gray boxes 3 and 4 allow to choose some features for the implemented algorithms. Box 3 is for the bandwidth size for the density kernel and box 4 is for the probability of Malo’s method, which is based on a Poisson distribution. The drop down indicated by box 5 is to select the type of map for visualization. Finally, in the space indicated by box 6 the map with the respective density kernel or hotspot segments (I’ve decided to use 500m segments) are indicated.

The packages comes with some sample data for testing. To find the location of these files run:

system.file("extdata", package = "roadHotspots")

After importing the data into the app you will get the algorithms computed and ready for visualization. From the below picture you can see that the water course near the town “Parvoíce” is the one with the greatest density of roadkils/observations.

hot.png

By default the topography map is loaded with both the Malo’s method and the kernel density hotspot algorithms loaded. The kernel is recoded into 3 classes, red, yellow and green (being red the most dense). For Malo’s method the black road segments are the most dense ones. In the check boxes on the top right you can select the method you want to see. Also you can change the map view aswell. Below image shows the topographic view.

topo.png

So, in conclusion: you don’t need ninja programming skills nor absurd funding  to create useful and easy to use softwares. With a bit of motivation and a few hours of work a lot can be achieved nowadays.

Just a cautionary note: all of this was done with free software and I don’t guarantee that it will work in Windows.

Just accepted: “Bad moon rising? The Influence of the lunar cycle on amphibian roadkills”

Does the lunar cycle has some sort of effect on roadkill (wildlife killed on the road)? Is this evident for any amphibian species? Yes! Despite the overwhelming importance of weather-related variables for some species, like Pleurodeles waltl and Salamandra salamandra, the lunar cycle affects roadkill numbers!

This was a great work in cooperation with Helena Lopes, Tiago Pinto, Luís Guilherme Sousa, António Mira and Sara Santos.

Here’s the abstract!

Annually, roads, and their associated users, are responsible for millions of roadkills worldwide. Mortality affects multiple taxonomic groups, but amphibians are particularly vulnerable, due to their size and under-reporting. In fact, very high mortality frequencies can occur, mostly during short periods of time, when individuals migrate to and from reproduction areas (e.g., ponds). In this study we assess the influence of the lunar cycle on amphibian roadkills, while accounting for weather conditions. As expected, the main environmental effects explaining roadkill numbers were weather-related, with increases in minimum air temperature, average relative air humidity and cumulative rainfall during the previous 24 h having a positive effect on roadkill numbers for all studied species. However, the lunar cycle also affected roadkills for two of the studied species. Darker nights had higher numbers of roadkills of Pleurodeles waltl, while moonlit nights had higher numbers of Salamandra salamandra. As such, these moon effects are species-specific. Animals that are more active in moonlight may be at an advantage if their visual acuity is better than that of their predators. We hypothesize that differences between species in the response to moonlight may be due to differences perceived in predation risk. This information should be considered when designing mitigation measures. Volunteer actions, for instance, can be planned and coordinated keeping in mind the most appropriate weather conditions for the general amphibian community and specific phases of the lunar cycle for particular species.

 

 

R as GIS for ecologists

Working with spatial data is a key feature in ecological research. Using R to handle this type of data has the great advantage of keeping both variable extraction and modelling in the same environment, instead of recurring to external GIS softwares to compute some variables and then turning to R for modelling. In this example I’ll use as my base data an altimetry contour line map, then I’ll compute a DEM (digital elevation model) from the altimetry contour lines, derive some maps with variables related to altimetry and finally I’ll measure the variables for sampling sites.

In this example packages “rgdal” and “rgeos” are used to perform vectorial operations, package “raster” to obtain and work with raster maps package “gstat” to interpolation and packages “rgl” and “rasterVis” to visualize 3D plots.

library(raster)
library(rgdal)
library(rgeos)
library(gstat)
library(rgl)
library(rasterVis)

 

The base data can be downloaded here, and is composed of: 1) altimetry, 2) sampling points and 3) study area boundary. As the data is in vector format it can be imported with readOGR().

altimetry <- readOGR(dsn = ".", layer = "altimetry", verbose = F)
study_area <- readOGR(dsn = ".", layer = "study_area", verbose = F)
sampling_points <- readOGR(dsn = ".", layer = "sampling_points", verbose = F)

 

To visualize the data is best to plot all the different information in one plot.

# Plot contour lines
plot(altimetry, col = "dark blue")

# Add study area boundary
plot(study_area, add = TRUE)

# Add sampling points
points(sampling_points, col = "red", cex = 0.5)

fig1-1.png

Spatial data manipulation can be very time consuming. So, it’s wise to clip the altimetry map to the study area to increase speed in further computations. Function intersect() can be used for this end. The new altimetry map is shown below.

# Clip altimetry map to study area
altimetry_clip <- intersect(altimetry, study_area)

# Plot clipped altimetry map
plot(altimetry_clip, col = "dark blue")

# Add sampling points
points(sampling_points, col = "red", cex = 0.5)

fig2-1.png

To create a DEM from altimetry contour lines the following steps are needed: 1) Create a blank raster grid to interpolate the elevation data onto; 2) Convert the contour lines to points so you can interpolate between elevation points; 3) Interpolate the elevation data onto the raster grid.

  1. First we need to create a blank raster grid. The extent and the projection of the raster should be the same as the altimetry map, so we are going to use the information of this shape for our new raster grid. Afterwards the pixel size of the raster is also defined. In this example I’ll use a 5m x 5m pixel.
# Obtain extent
dem_bbox <- bbox(altimetry_clip)

# Create raster
dem_rast <- raster(xmn = dem_bbox[1, 1], 
                   xmx = ceiling(dem_bbox[1, 2]),
                   ymn = dem_bbox[2, 1], 
                   ymx = ceiling(dem_bbox[2, 2]))

# Set projection
projection(dem_rast) <- CRS(projection(altimetry_clip))

# Set resolution
res(dem_rast) <- 5

 

  1. Since interpolation methods were conceived to work with point data, we need to convert the elevation contour lines to elevation points. Essentially we are creating points along each contour line that has as its value the elevation of the associated contour line.
# Convert to elevation points
dem_points <- as(altimetry_clip, "SpatialPointsDataFrame")

 

  1. To perform the interpolation of the point data one two methods are widely used: Nearest Neighbor and Inverse Distance Weighted. The difference between the two methods is that in nearest neighbor all the surrounding points have the same weight. In inverse distance weighted points that are further away get less weight. The function used is the same, gstat(), but for the nearest neighbor methods the argument “idp”, the inverse distance power, must equal zero. For inverse distance some value of idp should be set.
# Compute the interpolation function
dem_interp <- gstat(formula = ALT ~ 1, locations = dem_points,
 set = list(idp = 0), nmax = 5)

# Obtain interpolation values for raster grid
DEM <- interpolate(dem_rast, dem_interp)

 

Now that we have our elevation model ready we can plot it in 2D with some contour lines added for better visualization.

# Subset contour lines to 20m to enhance visualization
contour_plot <- altimetry_clip[(altimetry_clip$ALT) %in% 
                                 seq(min(altimetry_clip$ALT), 
                                     max(altimetry_clip$ALT), 
                                     20), ]                                                                   
# Plot 2D DEM with contour lines
plot(DEM, col = terrain.colors(20))
plot(contour_plot, add = T)

fig3.png

Or make an interactive 3D plot that can be controlled with the mouse.

plot3D(DEM)

 

With the DEM ready other altitude related variables can be derived, such as slope, aspect or roughness, among others. As aspect is a circular variable, i. e., the minimum value (0º) and the maximum value (360º) represent the same thing (north), a better way of using this information is to convert it into two new variables: northness = cos(aspect) and eastness = sin(aspect). Please take note that aspect must be in radians (radian = degree * pi / 180). Northness will take values close to 1 if the aspect is mostly northward, close to -1 if the aspect is southward, and close to 0 if the aspect is either east or west. Eastness behaves similarly, except that values close to 1 represent east-facing slopes. Other approach would be to reclassify aspect into N, S, E and W.

# Obtain DEM derived maps
derived_vars <- terrain(DEM, opt = c('slope', 'roughness', 'aspect'), unit = "degrees")
slope <- derived_vars[["slope"]]
roughness <- derived_vars[["roughness"]]
northness <- cos(derived_vars[["aspect"]] * pi / 180)
eastness <- sin(derived_vars[["aspect"]] * pi / 180)

# Plot maps
par(mfrow = c(2, 2))
plot(slope, col = heat.colors(20), main = "slope", axes = F)
plot(roughness, col = heat.colors(20), main = "roughness", axes = F)
plot(northness, col = heat.colors(20), main = "northness", axes = F)
plot(eastness, col = heat.colors(20), main = "eastness", axes = F)

fig4.png

Having all the maps prepared, last thing to do is to measure the values of these maps in the sampling points and create a data frame that can be used in further modelling. A very simple way would be to just measure the values in the exact location of the points, but a better way is to use a buffer around the sampling points and summarize the values in the buffer, using mean, mode, max, etc.

# Create buffers with 100m radius around sampling points
sp_buff <- gBuffer(sampling_points, width = 100, 
                   quadsegs = 50, byid = TRUE)

# Measure variables in buffer area
sp_alt <- extract(DEM, sp_buff, fun = mean)
sp_slope <- extract(slope, sp_buff, fun = mean) 
sp_rough <- extract(roughness, sp_buff, fun = mean)
sp_north <- extract(northness, sp_buff, fun = mean)
sp_east <- extract(eastness, sp_buff, fun = mean)

# Prepare dataframe
results <- data.frame("id" = sp_buff@data$id,
                      "altitude" = sp_alt,
                      "slope" = sp_slope,
                      "roughness" = sp_rough,
                      "northness" = sp_north,
                      "eastness" = sp_east)
# View results
print(results)
##   id altitude     slope roughness   northness   eastness
## 1  1 487.2113  1.980742 0.4936204 -0.06990387  0.9496363
## 2  2 637.9079 26.458705 6.5107228  0.25682424 -0.9029607
## 3  3 534.3797 28.770571 7.0379747  0.93539489  0.2672350
## 4  4 537.3855 18.411266 4.4562798 -0.24734534  0.6485864
## 5  5 583.4479 27.955300 6.9689737 -0.28556347  0.8485821

 

Using R as GIS to handle spatial data brings not only the advantage of keeping both variable extraction and modelling in the same environment but also gives the possibility of creating a script to extract the variables, with the obvious advantages for future analysis.