Archive

Posts Tagged ‘map’

Getting and managing occurrence data from gbif.org

The R package dismo by Hijman et al. offers a very powerful function named gbif that allows to get occurence data from the gbif dabase.

For example if we are to map the geographical range of the Mediterranean fruit fly, Ceratitis capitata we simply do:

library(dismo)
distr <- gbif(genus="Ceratitis", species="capitata*", geo=TRUE,removeZeros = TRUE)

Note that “capitata” and “capitata*” won’t lead to similar results because the species name is often associated to author’s name in the database i.e. Ceratitis capitata or Ceratitis capitata (Wiedemann, 1824).
The results can be plotted:

plot(distr$lon,distr$lat, pch=3, cex=1, col="red", asp=1, xlab="longitude", ylab="latitude")
map("world", resolution = 0.5, add=T)

gbif occurrences for Ceratitis capitata

gbif occurrences for Ceratitis capitata

Alas! many occurrences correspond to the same geographical location i.e. there are many duplicated points which might be a source of trouble in subsequent data analyses and distribution modeling. This problem can be easily solved using the R package spatstat by Baddeley and Turner (2005) and its function duplicated.ppp.

We first have to create a spatial point pattern object from the occurrences points using the function ppp and a window created with the ripras function. Then duplicated.ppp will identify duplicated points and we will be able to remove them. The code is:

library(spatstat)
x<-distr$lon ; y<-distr$lat
w <- ripras(x,y)
wp<-ppp(x,y,window=w)
dupv<-duplicated.ppp(wp)
x2<-x[which(dupv==FALSE)] ; y2<-y[which(dupv==FALSE)] # coordinates of points with no duplicates

In the case of Ceratitis capitata gbif returned a total of 1279 occurrences that included not less that 982 duplicated points.


References
A. Baddeley and R. Turner (2005). Spatstat: an R package for analyzing spatial point patterns. Journal of Statistical Software 12 (6), 1-42. ISSN: 1548-7660. URL: http://www.jstatsoft.org

Robert J. Hijmans, Steven Phillips, John Leathwick and Jane Elith (2012). dismo: Species distribution modeling. R package version 0.7-23. http://CRAN.R-project.org/package=dismo



Advertisements
Categories: R Tags: ,

Import worldclim data into GRASS

This post is about how to import worldclim data sets into a GRASS project.

WorldClim – Global Climate Data WorldClim is a set of global climate layers (climate grids) with a spatial resolution of about 1 square kilometer. The data can be used for mapping and spatial modeling in a GIS or with other computer programs.
http://www.worldclim.org/

First you need to create a location using the EPSG code corrresponding to the worldclim data. It is 4326. You can find a list of the EPSG codes at http://spatialreference.org/

Once the location has been created you can launch GRASS with the location and import worldclim data sets. Worldclim offers different resolutions and you have to choose one of them. Then, define the g.region in GRASS using the selected resolution. For instance, we want data at a resolution of 00:02:30.
In GRASS type:

g.region n=90N s=90S e=180E w=180W res=00:02:30

This matches with the extent of the worldclim data and the desired resolution.

The last step is pretty straightforward. Select the directory where the worldclim data have been downloaded.
For instance:

cd /home/rossi/data/worldclim/ 

Then import the data (here we import the Altitude layer)

r.in.bin -s input=alt.bil output=altitude bytes=2 north=90 south=-60 east=180 west=-180 rows=3600 cols=8640 anull=-9999 

Note that the number of columns and rows as well as the east, west, north and south limits are given in the *.bil file available in the download from http://www.worldclim.org.

It is convenient to import automatically a set of layers using a loop :

for i in $(seq 1 19)
do
r.in.bin -s input=bio$i.bil output=bio$i bytes=2 north=90 south=-60 east=180 west=-180 rows=3600 cols=8640 anull=-9999
done

In that example, we have imported the 19 layers of the BIOCLIM variables using a loop and the GRASS function r.in.bin

Now the worldclim data are available within our GRASS location. These raster layers can be used visualisation, analyses and modelling with GRASS but also using R (see the post “Playing with R within a GRASS environment” from this blog).

Have fun!

Playing with R within a GRASS environment

14 August 2012 1 comment

Although there are some powerful GIS utilities in R, I prefer using GRASS to manage my GIS data while I use R to perform scientific computing.

In fact GRASS and R can be very simply interfaced by means of the R package spgrass6.

Under linux operating system (actually Kubuntu Natty Narwhal) I open the console, launch GRASS by typing grass and then select a location (see GRASS manual for details).

console

Once GRASS is running, R can be launched from within the GRASS console…

From now on I am working in R and I can use the library spgrass6 in order e.g. to read or write rasters into the GRASS system.

As an example we will consider the case of the pine shoot beetle Tomicus piniperda (Coleoptera: Curculionidae: Scolytinae). Readers can find the whole story in Horn et al. (2012) available here. Tomicus piniperdaT. piniperda is present throughout Europe. Interestingly, it has long been assumed to be present in North Africa too although molecular studies have recently shown that T. piniperda only rarely occurs in these regions (Horn et al., 2006, 2009).

We have a set of localities where T. piniperda was recorded as present or, on the contrary, has been searched and was absent (true absence). GRASS is used to host the data and produce the graphical outputs.

Sites where T. piniperda was either present (black circle) or absent (open circle).

Starting R from the GRASS console and using the dismo package allows us to fit a SDM (Species Distribution Model) and build a raster with the probabilities of presence e.g. using the function predict. This raster layer can be expressed as presence/absence and the resulting data written within the GRASS system using the function writeRAST6 from package spgrass6.
The raster can now be plotted using the GRASS interface and its utilities.

GRASS window showing the occurences of T. piniperda and the associated predicted distribution derived from a SDM run in R


Horn, A., Kerdelhué, C., Lieutier, F., Rossi, J.-P. (2012). Predicting the distribution of the two bark beetles Tomicus destruens and Tomicus piniperda in Europe and the Mediterranean region. Agricultural and Forest Entomology, in press, DOI: 10.1111/j.1461-9563.2012.00576.x.

Horn, A., Roux-Morabito, G., Lieutier, F. & Kerdelhué, C. (2006) Phylogeographic structure and past history of the circum-Mediterranean species Tomicus destruens Woll. (Coleoptera: Scolytinae). Molecular Ecology, 15, 1603–1615.

Horn, A., Stauffer, C., Lieutier, F. & Kerdelhué, C. (2009) Complex postglacial history of the temperate bark beetle Tomicus piniperda L. (Coleoptera, Scolytinae). Heredity, 103, 238–247.


Display pie charts with varying alpha transparencies upon a Google static map

It is sometimes interesting to display pie charts on a map showing geographic details such as Google static maps. This post will illustrate :

  • how to download and manipulate Google static maps using the R package {RgoogleMaps} by M. Loecher
  • how to overlay floating pie charts upon a map using the R package {plotrix} by J. Lemon
  • how to change the alpha value associated to colors {seqinr} by D. Charif, and J.R. Lobry

As an example we will focus on the town of Niamey where my colleagues Madougou Garba and Gauthier Dobigny (Centre Régional AGRHYMET, Département Formation Recherche, Niamey) are currently studying urban rodents communities. Lets say we want to display pie charts corresponding to sampling locations where the rodent community has been sampled using traps.

First we download the static map from Google server using the function GetMap from the package RgoogleMaps.


path<-"/home/rossi/data/blog/data/"
map1.n<-paste(path,"map1",".png",sep="")
zoom<-12 ; lonc<-2.114668 ; latc<-13.520508 ; center<-c(latc,lonc) ; maptype<-"satellite"
niamey <- GetMap(center=center, zoom=zoom, maptype=maptype,destfile=map1.n)

Note that lonc and latc are the coordinates of a points roughtly located at the center of the town. It indicates where the center of the static map should be located.

Once the data are available through the R object niamey we can display it on sreeen using PlotOnStaticMap

PlotOnStaticMap(niamey,destfile=map1.n)

Here is is what you get :

Google static map of Niamey

Now we come to the pie charts. Lets say we have a sample of the urban rodent community collected in a location which coordinates are latitude=13.51204 and longitude=2.09884. Imagine that 20 individuals of R. rattus and 5 individuals of M. natalensis were recorded at that location.

We can plot a point that shows the spatial location of the locality

PlotOnStaticMap(niamey,lat=13.51204, lon=2.09884, pch=3,cex = 3,col="red",verbose=0,destfile=map1.n)

Niamey

If we want a pie chart showing the proportion of each rodent species we use the function floating.pie(xpos,ypos,x,...) where xpos and ypos are the coordintates of the centre of the pie and x is a vector with rodent species counts.
The trick is that xpos and ypos should correspond to the map tile coordinates system and cannot be given in raw degrees.
The package RgoogleMaps provides the function LatLon2XY.centered allowing the conversion:

newcoord<-LatLon2XY.centered(niamey,lat=13.51204,lon=2.09884,zoom=zoom)

We project the pie chart using these new coordinates:

PlotOnStaticMap(niamey,destfile=map1.n)
floating.pie(xpos=newcoord$newX,ypos=newcoord$newY,x=c(20,5),
radius=20,col=c("red","blue"))

radius controls the size of the pie chart, and it can be adjusted to plot a set of pie chart with size proportional to the sample size for example.

When a large number of pie charts are plotted on the same map they may be superimposed. In that case having more or less transparent colors (i.e. changing the alpha) may be interesting.
This is achieved using the function col2alpha from the package seqinr.

col<-c("red","blue")
for (i in 1:(2)) {col[i]<-col2alpha(color=col[i], alpha = 0.50)}
PlotOnStaticMap(niamey,destfile=map1.n)
floating.pie(xpos=newcoord$newX,ypos=newcoord$newY,x=c(20,5),
radius=20,col=col)

Here is the final result, note that the river Niger can now be seen through the pie chart:

Pretty cool he?

Categories: R Tags: ,