Created for VertNet Biodiversity Informatics Training Workshop.

QuickSpeciesRichnessfromRanges.R

#############################################
###VertNet Biodiversity Informatics Training Workshop
###June 2013
#Calculating species richness from vector range maps
#with quick visualization. Introduction to spatial packages
#
###updated Michelle S. Koo, 20130612
###for R.3.0.1
#
##^ Package		Versions
##^ maptools	0.8-23
##^ rgdal 	 	0.8-8
##^ sp 			1.0-9
##^ scales 		0.2.3
##^ rgeos 	0.0-5
##^ raster	 	2.5
#################################################

install.packages(c('maptools','rgeos', 'rgdal', 'scales', 'sp', 'raster'))
library(rgdal);
library(maptools);
gpclibPermit();
library(sp);
library(scales);

#We are going to use some basic shapefiles from GADM (http://www.gadm.org/) for the  African continent,
#and from the IUCN (http://www.iucnredlist.org/) for the range maps of the Family Felidae.

setwd("~/Desktop/BiodivInformatics-Workshop") # set to your working directory
dsn="~/Desktop/BiodivInformatics-Workshop/shp" #for convenience in a few steps
cats_poly<-readShapePoly(file.choose(), proj4string=CRS('+proj=longlat +datum=WGS84'))

#browse for data, choose the file with .shp extension
africa.shp<-readShapePoly(file.choose())

#or if you want a basic world map...
data(wrld_simpl);

#check out the shapefile; what do these functions do?
summary(cats_poly)
attributes(cats_poly@data)
unique(cats_poly$name)


#now to map
plot(wrld_simpl, axes=TRUE, col="light grey")
plot(cats_poly, add=TRUE, col=alpha("orange", 0.4))

##calc richness from polygons without transformation:
library(rgeos)
library(raster)
cats<-readShapeSpatial((system.file("shp/cats.shp")), proj4string=CRS('+proj=longlat +datum=WGS84'))
bb<-gEnvelope(cats_poly) #create a rectangular polygon to sample points from later

#create a set of uniform points
grid<-spsample(x=bb,type='regular', n=1000) #changing n changes the grid resolution

#count number of polygons present for each point {sp}
require(sp)
x=sapply(over(grid, geometry(cats_poly), FUN='count', returnList = TRUE), length)
x=cbind(as.data.frame(grid),x)
cat_rich=rasterFromXYZ(x)

plot(cat_rich)
cellStats(cat_rich, stat='max', na.rm=TRUE)


########
#For felid richness in Africa:
proj4string(africa.shp) = proj4string(cats_poly) #we're using the Africa shapefile that we imported without defining a projection;
#now it's in the same projection as the 'cats_poly' spatial obj
bb1<-gEnvelope(africa.shp)
grid1<-spsample(x=bb1,type='regular',n=1000)
x1=sapply(over(grid1, geometry(cats_poly), FUN='count', returnList = TRUE), length) #'apply' family of functions can be very useful
x1<-cbind(as.data.frame(grid1),x1)
cats_africa=rasterFromXYZ(x1)

plot(cats_africa)
plot(bb1, add=T)
plot(africa.shp, add=T)

summary(cats_africa)

###################
#Other things to think about. We used more than one way to import files, including a command to interactively select a file.
#Which richness calculation is higher resolution, the global richness map or the African one? You know how many pixels are calculated,
#what would be the equivalent spatial resolution and how would you calculate that? How would you export these richness calculations
# into QGIS?
#######################
#Other R packages on R-Forge that gets you spatial data & analysis, webmapping, MODIS, etc
#https://r-forge.r-project.org/softwaremap/trove_list.php?form_cat=355

######################