\documentclass{article} % \VignetteIndexEntry{adehabitatMA: tools for the analysis of mapped data} \usepackage{graphicx} \usepackage[colorlinks=true,urlcolor=blue]{hyperref} \usepackage{color} \usepackage{url} \usepackage{Sweave} \newcommand{\strong}[1]{{\normalfont\fontseries{b}\selectfont #1}} \let\pkg=\strong \title{Tools for the Analysis of Mapped Data:\\ the {\tt adehabitatMA} Package } \author{Clement Calenge,\\ Office national de la classe et de la faune sauvage\\ Saint Benoist -- 78610 Auffargis -- France.} \date{Mar 2019} \begin{document} \maketitle \tableofcontents <>= oldopt <- options(width=80, warn=-1) .PngNo <- 0 wi <- 480 pt <- 20 @ <>= .PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".png", sep="") png(file=file, width = wi, height = wi, pointsize = pt) @ <>= dev.null <- dev.off() cat("\\includegraphics[height=7cm,width=7cm]{", file, "}\n\n", sep="") @ <>= dev.null <- dev.off() cat("\\includegraphics[height=14cm,width=14cm]{", file, "}\n\n", sep="") @ \section{History of the package adehabitatMA} The package {\tt adehabitatMA} contains functions allowing the analysis of mapped data that were originally available in the package {\tt adehabitat} (Calenge, 2006).\\ I developed the package {\tt adehabitat} during my PhD (Calenge, 2005) to make easier the analysis of habitat selection by animals. The package {\tt adehabitat} was designed to extend the capabilities of the package {\tt ade4} concerning studies of habitat selection by wildlife.\\ Since its first submission to CRAN in September 2004, a lot of work has been done on the management and analysis of spatial data in R, and especially with the release of the package {\tt sp} (Pebesma and Bivand, 2005). The package {\tt sp} provides classes of data that are really useful to deal with spatial data...\\ In addition, with the increase of both the number (more than 250 functions in Oct. 2008) and the diversity of the functions in the package \texttt{adehabitat}, it soon became apparent that a reshaping of the package was needed, to make its content clearer to the users. I decided to ``split'' the package {\tt adehabitat} into four packages:\\ \begin{itemize} \item {\tt adehabitatMA} package provides methods for dealing with maps in R. \item {\tt adehabitatHR} package provides classes and methods for dealing with home range analysis in R. \item {\tt adehabitatHS} package provides classes and methods for dealing with habitat selection analysis in R. \item {\tt adehabitatLT} package provides classes and methods for dealing with animals trajectory analysis in R.\\ \end{itemize} We consider in this document the use of the package {\tt adehabitatMA} to deal with the analysis of mapped data. All the methods available in \texttt{adehabitatMA} are also available in \texttt{adehabitat}, but the classes of data returned by the functions of \texttt{adehabitatMA} are completely different from the classes returned by the same functions in \texttt{adehabitat}. Indeed, the classes of data returned by the functions of \texttt{adehabitatHR} have been designed to be compliant with the classes of data provided by the package \texttt{sp}. Note that functions allowing the conversion from old classes of \texttt{adehabitat} to new classes of \texttt{adehabitatMA} and \texttt{sp} are described on the help page of the function \texttt{kasc2spixdf}. We therefore suppose that the user is familiar with this package.\\ Package {\tt adehabitatMA} is loaded by <>= library(adehabitatMA) @ <>= suppressWarnings(RNGversion("3.5.0")) set.seed(13431) adeoptions(shortprint=TRUE) @ \section{The aim of the package} \subsection{sp and adehabitatMA} The package \texttt{sp} is really useful for the analysis of mapped data, and for this reason, I encourage the user to become familiar with this package. However, the package \texttt{adehabitat} originally contained functions which were also useful in spatial analysis though not available in \texttt{sp}, and especially in the field of analysis of animal space use. For this reason, I included these functions in the package \texttt{adehabitatMA}. I demonstrate the use of these functions in this document.\\ I will demonstrate the use of this package using the example dataset \texttt{lynxjura}. <<>>= data(lynxjura) @ This dataset contains the results of the monitoring of the lynx in the Jura mountains (France) by the French lynx network of the \textit{Office national de la chasse et de la faune sauvage}. This dataset has two components: \texttt{\$locs} is a \texttt{SpatialPointsDataFrame} containing the locations of presence indices of the lynx, along with the date of collection and the type of indices (attacks on livestocks, captured animal, tracks, etc.). Look at the first rows of this object: <<>>= head(lynxjura$locs) @ The component \texttt{\$map} of this dataset is a \texttt{SpatialPixelsDataFrame} contains the maps of four environmental variables measured on the study area: <<>>= lynxjura$map @ The whole content of this object can be displayed using the function \texttt{mimage}: <>= mimage(lynxjura$map) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \section{Working with raster maps only} \subsection{Exploring a raster map interactively} \texttt{adehabitatMA} provides one function allowing the exploration of objects of class \texttt{SpatialPixelsDataFrame}, the function \texttt{explore}. This allows to interactively explore a map, to find the value of the variables at a given place, to measure the distance between two points, to zoom on a given part. Note that this function requires the package \texttt{tkrplot}: <>= explore(lynxjura$map) @ \includegraphics[keepaspectratio]{screenshot1.png} This function also has an argument named \texttt{panel.last}, which allows to pass an expression to be evaluated after plotting has taken place. For example, we can use it to explore the distribution of lynx indices in relationship with the maps. Copy and paste the following expression: <>= explore(lynxjura$map, panel.last=function() points(lynxjura$locs, pch=3)) @ \subsection{Labelling connected features on a map} It is sometimes useful to identify automatically separated components on a raster map. For example, consider the maps in the object \texttt{lynxjura} loaded previously: <<>>= map <- lynxjura$map @ Look at the values of the mapped variables: <>= hist(map) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Consider the variable ``forets'', which represents the percentage of forested area in each pixel of the map. We will consider areas with at least 95\% of forests: <>= forest <- map[,1] forest[[1]][forest[[1]]<95] <- NA image(forest, col="green") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We may wonder how many forested areas are identified using this criterium. in the present case, we can count them ourselves: 9 areas. However, the function \texttt{labcon} provides an efficient way to do it: <>= lab <- labcon(forest) image(lab) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The resulting object is a \texttt{SpatialPixelsDataFrame} mapping one integer vector taking a unique value per connected component. Therefore, we can count the number of components in the map. <<>>= lab max(lab[[1]]) @ It is easy to measure the area of each component, because we know the area covered by a pixel: <<>>= gridparameters(lab) @ Therefore, the area covered by each component can be computed by: <<>>= table(lab[[1]])*500*500 @ The result is here returned in squared meters as the original units of the maps were in meters.\\ We can also measure the average density of rivers (second variable, names ``hydro'' in the object \texttt{map}) in the component 1. We first have to transform the objects \texttt{lab} and \texttt{map} as full grid to allow the comparison: <<>>= fullgrid(lab) <- TRUE fullgrid(map) <- TRUE @ Then we compute the mean of the density of rivers in the component 1: <<>>= mean(map[[2]][lab[[1]]==1], na.rm=TRUE) @ Or getting the map of ``hydro'' for the first connected component in a separate map: <>= comp1 <- map[2] comp1[[1]][map[[1]]<95] <- NA comp1[[1]][lab[[1]]!=1] <- NA image(comp1) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \subsection{Computing the contour of a mapped area on a raster map} \label{sec:contour} Consider again the map of forested area built in the previous section: <>= image(forest, col="red") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} It may sometimes be useful to compute the coordinates of the vertices of the contour polygon of the forested area: <>= con <- getcontour(forest) @ \begin{Schunk} \begin{Sinput} > con <- getcontour(forest) \end{Sinput} \begin{Soutput} Warning message: In getcontour(forest) : At least 3 pixels are required to compute a contour. 3 components erased \end{Soutput} \end{Schunk} The resulting object is of the class \texttt{SpatialPolygons}. It can therefore be handled with other functions of the package \texttt{sp}, and can be exported toward a GIS using the functions of the package \texttt{sf} after conversion with \texttt{st\_as\sf} from this package.\\ However, note the warning here. To understand it, look at the resulting object: <>= plot(con, col="green") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} There are 6 polygons plotted here. However, we saw in the previous section that there are 9 forested areas in the studied regions. Three of the forested areas are built by less than three pixels, so that no contour can be calculated for these ``small forests''. \subsection{Mathematical morphology} In the last case, we may be annoyed to ``loose'' these three small forests. It is easy to understand that at least three points are required to build a polygon, and that the function cannot find the contour polygon of only one point (one pixel).\\ There is however a way to circumvent this drawback. We can use the function \texttt{morphology}, which performs operations of mathematical morphology: it allows to erode or to dilate the object (see the help page of the function. In our case, we may dilate each forested area by one pixel: <<>>= for1 <- morphology(forest, "dilate", nt=1) for1 @ The resulting object is a \texttt{SpatialPixelsDataFrame}. Now compare the area covered by the components of \texttt{for1} with the area covered by the components of \texttt{forest}: <>= image(for1, col="blue") image(forest, col="yellow", add=TRUE) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Each forested area has been expanded by exactly one pixel on all sides. Note that the components on the boundary of the map cannot be expanded. Now if we use the function \texttt{getcontour} again: <>= plot(getcontour(for1), col="green") @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} All forested areas are now present in this vectorized object. This operation of course assume that the area increase caused by the use of the function can be considered as negligible.\\ Note that the inverse operation (erosion) is also possible with the function \texttt{morphology}. \subsection{Changing the resolution of a map} The package \texttt{adehabitatMA} provides a function allowing to change the resolution of the map. For example, consider the map of the study area used previously: <>= map <- lynxjura$map mimage(map) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The parameters of this map are: <<>>= gridparameters(map) @ We may want to change the resolution of this map so that pixels cover 5 $\times$ 5 km instead of 500 $\times$ 500 m. That is, we want to merge together 10 pixels into one pixel. The function \texttt{lowres} allows to perform this operation. We first define \texttt{map} as a \texttt{SpatialPixelsDataFrame}: <>= mimage(lowres(map, 10)) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The function computes the average value in each pixel of the new map. However, note that in some case, it is not sensible to compute an average to summarize the values observed in a given pixel. For example, imagine that we transform the variable ``forets'' in the following way: <>= map[[1]] <- as.numeric(cut(map[[1]],3)) image(map, 1) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The map is now a factor with three classes. In such a case, it does not make sense to compute a mean to summarize the values of the pixels. In such a case, the function \texttt{lowres} assigns to the pixel the most frequent level. When several levels are equally represented in the pixel, the function randomly samples one of these levels. We have to indicate to the function which variable(s) in the object are factors thanks to the argument \texttt{which.fac}: <>= image(lowres(map, 10, which.fac=1)) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \subsection{Subsetting a map} Consider again the map of the forests: <>= image(forest, col="green") box() @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The forested areas are located in the eastern part of the area. We can use the function \texttt{subsetmap} to get that part of the map: <>= for2 <- subsetmap(forest) @ Then click on the map to indicate the lower-left and upper-right boundaries of the new map: <>= for2 <- subsetmap(forest, xlim=c(850254.2, 878990.2), ylim=c(2128744, 2172175)) @ <>= image(for2, col="green") box() @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \section{Working with points and maps} We now consider the case where we work with points and raster maps. Remember that the functions of \texttt{adehabitatMA} were originally available in \texttt{adehabitat} to study issues related to habitat selection by animals. Issues involving both points and maps are extremely frequent in this field. The points may be the relocations of animals collected using radio-tracking, or presence indices of animals on an area (as is the case for the \texttt{lynxjura} data). The maps therefore represent the environment of the animals.\\ We programmed two commonly used functions to deal with points and maps:\\ \begin{itemize} \item \texttt{count.points} allows to count the number of points in each pixel of the raster map; \item \texttt{join} allows to identify the values of mapped variables at specified locations.\\ \end{itemize} Although the function \texttt{overlap} of the package \texttt{sp} could easily be used to perform these operations, we prefered to provide specific functions, as these operations are very common in the field of habitat selection studies. \subsection{Counting the number of points in pixels} Consider the data on the lynx in the Jura mountains: <<>>= data(lynxjura) map <- lynxjura$map class(map) locs <- lynxjura$loc class(locs) @ Display the image: <>= image(map, 1) points(locs, pch=3) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We can count the number of points in each pixel of the map: <>= cp <- count.points(locs, map) @ \begin{Schunk} \begin{Sinput} > cp <- count.points(locs, map) \end{Sinput} \begin{Soutput} Warning message: In count.points(locs, map) : several columns in the SpatialPointsDataFrame, no id considered \end{Soutput} \end{Schunk} Note the warning message. We will talk about it later. Now, display the result: <>= image(cp) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} The warning returned by the functions is explained by the object containing the points. Indeed, this object should be, according to the help page, an object of class \texttt{SpatialPoints}, or \texttt{SpatialPointsDataFrame} \textit{with one column}. In the latter case, the column is considered as a factor giving, for each point, the membership of the point to a set (e.g. the identity of the animal). In this case, the \texttt{SpatialPointsDataFrame} contains more than one column: <<>>= head(locs) @ The function does not know which column should be considered as the variable defining the membership of points to a set. Therefore, it treats the object as a \texttt{SpatialPoints} object (which is the desired behaviour). However, we could have wanted to compute several maps, with one map per type of index, indicating the number of indices of each type in each pixel of the map: <<>>= cpr <- count.points(locs[,"Type"], map) cpr @ The results can be displayed: <>= mimage(cpr) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} There is one map per type of points. \subsection{Finding the value of mapped variables at specified locations} The other possible approach consists to perform a spatial join: in other words, to find the value of mapped variables at specified locations For example: <<>>= df <- join(locs, map) @ The result is a data frame containing the value of the variables at the specified locations: <<>>= head(df) @ This data frame can then be used in further analyses of habitat selection. \subsection{Generating a raster map from points data} In some cases it may be useful to generate a grid of pixels from a set of points. The package \texttt{sp} contains a function named \texttt{makegrid}, which generates a regular grid of points. However, I also included in \texttt{adehabitatMA} a function which computes directly a \texttt{SpatialPixelsDataFrame} from a set of points. This function allows to specify either the size of the pixel or the number of rows/columns of the grid. For example, to generate a \texttt{SpatialPixelsDataFrame} object from the set of locations of presence indices of the lynx, with a pixel size of 5000 m $\times$ 5000 m: <<>>= asc <- ascgen(locs, cellsize=5000) asc @ The resulting object is an object of class \texttt{SpatialPixelsDataFrame}, with one column containing the number of points in the pixels of the map. <>= image(asc) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \section{Computing buffer regions} \texttt{adehabitatMA} provides a function allowing the computation of buffer region. Buffer regions can be computed from points, lines or polygons. For example, consider the presence indices of the lynx of type ``sighting'' (labelled "O" in the data) <<>>= po <- locs[locs[["Type"]]=="O",] @ We may want to identify all areas located within 3000 m from a lynx presence index: <>= image(buffer(po, map, 3000)) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Now consider again the vector map of forested areas that we built in the section \ref{sec:contour} (map \texttt{con}): <>= plot(con) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} We may want to identify all the areas located within 3000 metres from one of these forested areas: <>= image(buffer(con, map, 3000)) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} Alternatively, we may have wanted to identify the ``ecotone'' (boundary habitat) between the forest and the open areas... for example, to identify all the areas located within 500 metres from the boundary: <>= sl <- as(con, "SpatialLines") image(buffer(sl, map, 500)) @ <>= <> @ \begin{center} <>= <> <> <> @ \end{center} \section{Conclusion} I included in the package \texttt{adehabitatMA} several functions adding to the existing set of tools available in R to perform spatial analyses. Originally, these functions were included in the package \texttt{adehabitat} to allow the study of very specific ecological issues related to habitat selection, but as most users who sent me questions/suggestions concerning these functions were using them to study issues in other fields (e.g. biogeography, landscape ecology, etc.), I decided to include them in a separate package much more compliant with existing packages for spatial analysis within R.\\ However, readers interested in the study of animal space use and habitat selection should also have a look at the brother packages. All the packages \texttt{adehabitat**} contain a vignette similar to this one, which explains not only the functions, but also in some cases the philosophy underlying the analysis of animal space use. \section*{References} \begin{description} \item Calenge, C. 2005. Des outils statistiques pour l'analyse des semis de points dans l'espace ecologique. Universite Claude Bernard Lyon 1. \item Calenge, C. 2006. The package adehabitat for the R software: a tool for the analysis of space and habitat use by animals. Ecological modelling, 197, 516--519. \item Pebesma, E. and Bivand, R.S. 2005. Classes and Methods for Spatial data in R. R News, 5, 9--13. \end{description} <>= options(oldopt) @ \end{document} % vim:syntax=tex